Is it possible to express "t" variable from Cubic Bezier Curve equation?
I want draw bezier curve only by fragment shader to connect nodes in my editor. I know all 4 points that define Bezier Curve. And Fragment Shader is called for every pixel, so i can just check: if "t" for gl_Coord.x is between 0 and 1 then set frag_color to Red for example. I want to avoid loops in shader that's inefficient. Best way, i think, is check for points that lay on the curve. But how to do it for Bezier Curves?
Is it possible to express "t" variable from cubic bezier equation?
x = ((1-t)^3 * p0.x) + (3 * (1-t)^2 * t * p1.x) + (3 * (1 - t) * t^2 * p2.x) + (t^3 * p3.x);
t = ?
Website Wolfram Aplha give me that formula(in GetBezierT function). But formula give me wrong "t" values and i have half of parabola instead of curve:
#version 150
.....
layout (origin_upper_left, pixel_center_integer) in vec4 gl_FragCoord;
out vec4 frag_color;
.....
vec4 BackgroundColor = vec4(0.15, 0.15, 0.15, 1.0);
vec2 p0 = vec2(61.0f,87.0f);
vec2 p1 = vec2(181.0f, 39.0f);
vec2 p2 = vec2(283.0f, 178.0f);
vec2 p3 = vec2(416.0f, 132.0f);
float getBezierT(float x, float a, float b, float c, float d)
{
return float(sqrt(3) *
sqrt(-4 * b * d + 4 * b * x + 3 * c * c + 2 * c * d - 8 * c * x - d * d + 4 * d * x)
+ 6 * b - 9 * c + 3 * d)
/ (6 * (b - 2 * c + d));
}
void main() {
.....
frag_color = BackgroundColor;
.....
float tx = getBezierT(gl_FragCoord.x, p0.x, p1.x, p2.x, p3.x);
float ty = getBezierT(gl_FragCoord.y, p0.y, p1.y, p2.y, p3.y);
if (tx >= 0.0f && tx <= 1.0f && ty >= 0.0f && ty <= 1.0f)
{
if(abs(tx-ty) < 0.01f) // simple check is that one point with little bias
frag_color = vec4(1.0f, 0.0f, 0.0f, 1.0f);
}
}
UPDATE
Made a mistake. I thought there was no point in looking for t. I thought I would put up with it. But after the answer given by Salix alba and Stratubas, I realized that if tX is equal to tY, this means that this point will lie on the curve, because in the formula for each point one value of t is substituted for both x and y. Maybe there are cases when different tX and tY can also give a point on this curve, but we can just ignore that. The algorithm for constructing a bezier curve implies that we linearly increase t and substitute it into the formula and it does not matter how much the curve is twisted, the algorithm returns the coordinates of each next point sequentially along the curve.
Therefore, first of all, I again open this question: how to express the variable t from a cubic bezier equation?
Tried to express t, but it is insanely difficult for me. It is necessary to evaluate the effectiveness of this approach for "scientific purposes" =). Before asking a question here, I searched a lot, but never found that someone would try to use this method. I need to understand why.
UPDATE 2
You have done an excellent job! I did not expect to receive such detailed answers. Exactly what i needed. Give me time to check everything=)
UPDATE 3
Conclusions: Accurate expression of t from the Cubic Bezier equation. Time-consuming task, but approximate values don't have practical use. To solve this problem, it is necessary to analyze the equation data, find patterns and develop new formula for constructing bezier curves. With a new relations of variables among themselves, then it will become possible to express t in a different way. If we represent the Cubic Bezier formula in the form of the sum of the products of the x coordinates of the control points by four coefficients ( v0 -v3) generated by the functions in the four parts of the equation depending on the value of t. This gives the formula x = a.x * v0 + b.x * v1 + c.x * v2 + d.x * v3. And if you look at the table below, you can get the idea that the expression for the variable t is an equation with four unknowns. Because both the values and the relations of some of the V coefficients between themselves change in an unpredictable way from iteration to iteration. Finding that new abstract formula is beyond the scope of this question and my competence.
Many thanks to all for your work, especially Spektre for the unique development and efforts made to optimize the rendering algorithm. Your approach is the best choice for me=)
Solution 1:
What you need is to search your cubic path and remember closest point. This can be done recursively with increasing precisions here small C++ GL example:
//---------------------------------------------------------------------------
double pnt[]= // cubic curve control points
{
-0.9,-0.8,0.0,
-0.6,+0.8,0.0,
+0.6,+0.8,0.0,
+0.9,-0.8,0.0,
};
const int pnts3=sizeof(pnt)/sizeof(pnt[0]);
const int pnts=pnts3/3;
//---------------------------------------------------------------------------
double cubic_a[4][3]; // cubic coefficients
void cubic_init(double *pnt) // compute cubic coefficients
{
int i;
double *p0=pnt,*p1=p0+3,*p2=p1+3,*p3=p2+3;
for (i=0;i<3;i++) // cubic BEZIER coefficients
{
cubic_a[0][i]= ( p0[i]);
cubic_a[1][i]= (3.0*p1[i])-(3.0*p0[i]);
cubic_a[2][i]= (3.0*p2[i])-(6.0*p1[i])+(3.0*p0[i]);
cubic_a[3][i]=( p3[i])-(3.0*p2[i])+(3.0*p1[i])-( p0[i]);
}
}
//---------------------------------------------------------------------------
double* cubic(double t) // return point on cubic from parameter
{
int i;
static double p[3];
double tt=t*t,ttt=tt*t;
for (i=0;i<3;i++)
p[i]=cubic_a[0][i]
+(cubic_a[1][i]*t)
+(cubic_a[2][i]*tt)
+(cubic_a[3][i]*ttt);
return p;
}
//---------------------------------------------------------------------------
double cubic_d(double *p) // return closest distance from point to cubic
{
int i,j;
double t,tt,t0,t1,dt,
l,ll,a,*q;
tt=-1.0; ll=-1.0; t0=0.0; t1=1.001; dt=0.05;
for (j=0;j<3;j++)
{
for (t=t0;t<=t1;t+=dt)
{
q=cubic(t);
for (l=0.0,i=0;i<3;i++) l+=(p[i]-q[i])*(p[i]-q[i]);
if ((ll<0.0)||(ll>l)){ ll=l; tt=t; }
}
t0=tt-dt; if (t0<0.0) t0=0.0;
t1=tt+dt; if (t1>1.0) t1=1.0;
dt*=0.2;
}
return sqrt(ll);
}
//---------------------------------------------------------------------------
void gl_draw()
{
int i;
double t,p[3],dp;
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glEnable(GL_CULL_FACE);
// GL render
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glDisable(GL_DEPTH_TEST);
glColor3f(0.2,0.2,0.2); glBegin(GL_LINE_STRIP); for (i=0;i<pnts3;i+=3) glVertex3dv(pnt+i); glEnd();
glPointSize(5); glColor3f(0.0,0.0,0.7); glBegin(GL_POINTS); for (i=0;i<pnts3;i+=3) glVertex3dv(pnt+i); glEnd(); glPointSize(1);
cubic_init(pnt);glColor3f(0.2,0.7,0.7); glBegin(GL_LINE_STRIP); for (t=0.0;t<1.001;t+=0.025) glVertex3dv(cubic(t)); glEnd();
glColor3f(0.0,0.7,0.0); glBegin(GL_POINTS);
p[2]=0.0; dp=0.01;
for (p[0]=-1.0;p[0]<1.001;p[0]+=dp)
for (p[1]=-1.0;p[1]<1.001;p[1]+=dp)
if (cubic_d(p)<0.05)
glVertex3dv(p);
glEnd();
glFlush();
SwapBuffers(hdc);
}
//---------------------------------------------------------------------------
so first you call cubic_init once to compute the coefficients and then to obtain point on curve as function of parameter use:
double pnt[3] = cubic(double t);
Now the reverse (I return closest distance ll but you can easily change it to return the tt)
double dist = cubic_d(double pnt[3]);
Now you just port this to shader and determine if fragment is close enough to curve to render it (hence the distance instead of t also for speed you can get rid of the last sqrt and use powered values latter).
The gl_draw function renders control points(blue) / lines(gray) the bezier curve (aqua) with GL and then emulates fragment shader to render curve with thickness 2*0.05 in (green)...
Preview:
Now its just a matter of porting that into GLSL. In order to use GLSL native way of passing vertexes you need to enlarge the area a bit like in here:
- Draw Quadratic Curve on GPU
But you need to change the geometry a bit to account for 4 control points instead of just 3. That stuff should be in geometry shader ...
So in geometry shader you should do the cubic_init, and in fragment shader discard if distance cubic_d is bigger than thickness.
The search is based on:
- How approximation search works
which I develop for problems like this. The search loop itself can be tweaked a bit to improve performance/precision ... but beware the initial search should sample the curve to at least 4-5 chunks otherwise it might stop working properly for some shapes.
[Edit1] after some thinking here the GLSL version
Vertex
// Vertex
#version 400 core
layout(location = 0) in vec2 pos; // control points (QUADS)
layout(location = 3) in vec3 col; // color
out vec2 vpos;
out vec3 vcol;
void main()
{
vpos=pos;
vcol=col;
gl_Position=vec4(pos,0.0,1.0);
}
Geometry:
//------------------------------------------------------------------------------
// Geometry
//------------------------------------------------------------------------------
#version 400 core
layout(lines_adjacency) in;
layout(triangle_strip, max_vertices = 4) out;
uniform float d=0.05; // half thickness
in vec2 vpos[];
in vec3 vcol[];
out vec2 a0,a1,a2,a3; // cubic coefficients
out vec3 fcol; // color
out vec2 fpos; // position
//------------------------------------------------------------------------------
void main()
{
vec4 p0,p1,p2,p3,a,b;
p0=gl_in[0].gl_Position;
p1=gl_in[1].gl_Position;
p2=gl_in[2].gl_Position;
p3=gl_in[3].gl_Position;
// compute BEZIER coefficients
a0.x= ( p0.x);
a1.x= (3.0*p1.x)-(3.0*p0.x);
a2.x= (3.0*p2.x)-(6.0*p1.x)+(3.0*p0.x);
a3.x=(p3.x)-(3.0*p2.x)+(3.0*p1.x)-( p0.x);
a0.y= ( p0.y);
a1.y= (3.0*p1.y)-(3.0*p0.y);
a2.y= (3.0*p2.y)-(6.0*p1.y)+(3.0*p0.y);
a3.y=(p3.y)-(3.0*p2.y)+(3.0*p1.y)-( p0.y);
// compute BBOX
a=p0; b=p0;
if (a.x > p1.x) a.x=p1.x; if (b.x < p1.x) b.x=p1.x;
if (a.x > p2.x) a.x=p2.x; if (b.x < p2.x) b.x=p2.x;
if (a.x > p3.x) a.x=p3.x; if (b.x < p3.x) b.x=p3.x;
if (a.y > p1.y) a.y=p1.y; if (b.y < p1.y) b.y=p1.y;
if (a.y > p2.y) a.y=p2.y; if (b.y < p2.y) b.y=p2.y;
if (a.y > p3.y) a.y=p3.y; if (b.y < p3.y) b.y=p3.y;
// enlarge by d
a.x-=d; a.y-=d;
b.x+=d; b.y+=d;
// pass it as QUAD
fcol=vcol[0];
fpos=vec2(a.x,a.y); gl_Position=vec4(a.x,a.y,0.0,1.0); EmitVertex();
fpos=vec2(a.x,b.y); gl_Position=vec4(a.x,b.y,0.0,1.0); EmitVertex();
fpos=vec2(b.x,a.y); gl_Position=vec4(b.x,a.y,0.0,1.0); EmitVertex();
fpos=vec2(b.x,b.y); gl_Position=vec4(b.x,b.y,0.0,1.0); EmitVertex();
EndPrimitive();
}
//------------------------------------------------------------------------------
Fragment:
// Fragment
#version 400 core
uniform float d=0.05; // half thickness
in vec2 fpos; // fragment position
in vec3 fcol; // fragment color
in vec2 a0,a1,a2,a3; // cubic coefficients
out vec4 col;
vec2 cubic(float t) // return point on cubic from parameter
{
float tt=t*t,ttt=tt*t;
return a0+(a1*t)+(a2*tt)+(a3*ttt);
}
void main()
{
vec2 p;
int i;
float t,tt,t0,t1,dt,l,ll;
tt=-1.0; ll=-1.0; dt=0.05; t0=0.0; t1=1.0; l=0.0;
for (i=0;i<3;i++)
{
for (t=t0;t<=t1;t+=dt)
{
p=cubic(t)-fpos;
l=length(p);
if ((ll<0.0)||(ll>l)){ ll=l; tt=t; }
}
t0=tt-dt; if (t0<0.0) t0=0.0;
t1=tt+dt; if (t1>1.0) t1=1.0;
dt*=0.2;
}
if (ll>d) discard;
col=vec4(fcol,1.0); // ll,tt can be used for coloring or texturing
}
It expect 4 BEZIER control points per CUBIC in form of GL_LINES_ADJACENCY since GL_QUADS are no more :( When I use it like this (inside gl_draw):
glUseProgram(prog_id); // use our shaders
i=glGetUniformLocation(prog_id,"d"); // set line half thickness
glUniform1f(i,0.02);
glColor3f(0.2,0.7,0.2); // color
glBegin(GL_LINES_ADJACENCY);
for (i=0;i<pnts3;i+=3)
glVertex3dv(pnt+i);
glEnd();
glUseProgram(0);
The result looks like this:
and of coarse is a lot of faster than the old api dotted shader emulation :). I know old api and new style GLSL shaders should not be mixed so you should create VAO/VBO instead of using glBegin/glEnd... I am too lazy to do that just for the purpose of this answer ...
Here the non function (more y per single x) example (compared with the CPU side dots):
double pnt[]= // cubic curve control points
{
+0.9,-0.8,0.0,
-2.5,+0.8,0.0,
+2.5,+0.8,0.0,
-0.9,-0.8,0.0,
};
As you can see both approaches matches the shape (dots used bigger thickness). In order this to work the search coefficients (dt) must be set properly to not miss a solution...
PS solving the cubic your way leads to 2 set of these:
Which I strongly doubt can be much computed faster than simple search.
[Edit2] further improvements
I simply changed the geometry shader so that it sample the curve into 10 segments and emit BBOX for each separatelly eliminating a lot of empty space that needed to be processed before. I changed the color layout and rendering order a bit.
This is new result (identical to previous one but several times faster due lower empty space ratio):
This is how the coverage looks now:
Before the coverage was BBOX of control points + enlargement by d which in this case was much bigger then curve itself (2 control points are outside view).
Here updated Geometry shader:
//------------------------------------------------------------------------------
// Geometry
//------------------------------------------------------------------------------
#version 400 core
layout(lines_adjacency) in;
layout(triangle_strip, max_vertices = 40) out; // 4*n <= 60
uniform float d=0.05; // half thickness
in vec2 vpos[];
in vec3 vcol[];
out vec2 a0,a1,a2,a3; // cubic coefficients
out vec3 fcol; // color
out vec2 fpos; // position
//------------------------------------------------------------------------------
vec2 cubic(float t) // return point on cubic from parameter
{
float tt=t*t,ttt=tt*t;
return a0+(a1*t)+(a2*tt)+(a3*ttt);
}
//------------------------------------------------------------------------------
void main()
{
float t,dt=1.0/10.0; // 1/n
vec2 p0,p1,p2,p3,a,b;
p0=gl_in[0].gl_Position.xy;
p1=gl_in[1].gl_Position.xy;
p2=gl_in[2].gl_Position.xy;
p3=gl_in[3].gl_Position.xy;
// compute BEZIER coefficients
a0.x= ( p0.x);
a1.x= (3.0*p1.x)-(3.0*p0.x);
a2.x= (3.0*p2.x)-(6.0*p1.x)+(3.0*p0.x);
a3.x=(p3.x)-(3.0*p2.x)+(3.0*p1.x)-( p0.x);
a0.y= ( p0.y);
a1.y= (3.0*p1.y)-(3.0*p0.y);
a2.y= (3.0*p2.y)-(6.0*p1.y)+(3.0*p0.y);
a3.y=(p3.y)-(3.0*p2.y)+(3.0*p1.y)-( p0.y);
p1=cubic(0.0);
for (t=dt;t < 1.001;t+=dt)
{
p0=p1; p1=cubic(t);
// compute BBOX
a=p0; b=p0;
if (a.x > p1.x) a.x=p1.x; if (b.x < p1.x) b.x=p1.x;
if (a.y > p1.y) a.y=p1.y; if (b.y < p1.y) b.y=p1.y;
// enlarge by d
a.x-=d; a.y-=d;
b.x+=d; b.y+=d;
// pass it as QUAD
fcol=vcol[0];
fpos=vec2(a.x,a.y); gl_Position=vec4(a.x,a.y,0.0,1.0); EmitVertex();
fpos=vec2(a.x,b.y); gl_Position=vec4(a.x,b.y,0.0,1.0); EmitVertex();
fpos=vec2(b.x,a.y); gl_Position=vec4(b.x,a.y,0.0,1.0); EmitVertex();
fpos=vec2(b.x,b.y); gl_Position=vec4(b.x,b.y,0.0,1.0); EmitVertex();
EndPrimitive();
}
}
//------------------------------------------------------------------------------
My gfx card has 60 vertex limit so as I output triangle strips emulating QUADs the limit on segments is 60/4 = 15 I used n=10 just to be sure it runs on lower HW. In order to change the number of segments see the 2 lines with comment containing n
[Edit3] even better coverage useful/empty space ratio
I changed the AABB BBOX coverage to ~OOB BBOX without overlaps. This also allows to pass actual range of t into fragment speeding up the search ~10 times. Updated shaders:
Vertex:
// Vertex
#version 400 core
layout(location = 0) in vec2 pos; // control points (QUADS)
layout(location = 3) in vec3 col; // color
out vec2 vpos;
out vec3 vcol;
void main()
{
vpos=pos;
vcol=col;
gl_Position=vec4(pos,0.0,1.0);
}
Geometry:
//------------------------------------------------------------------------------
// Geometry
//------------------------------------------------------------------------------
#version 400 core
layout(lines_adjacency) in;
layout(triangle_strip, max_vertices = 40) out; // 4*n <= 60
uniform float d=0.05; // half thickness
in vec2 vpos[];
in vec3 vcol[];
out vec2 a0,a1,a2,a3; // cubic coefficients
out vec3 fcol; // color
out vec2 fpos; // position
out vec2 trange; // t range of chunk
//------------------------------------------------------------------------------
vec2 cubic(float t) // return point on cubic from parameter
{
float tt=t*t,ttt=tt*t;
return a0+(a1*t)+(a2*tt)+(a3*ttt);
}
//------------------------------------------------------------------------------
void main()
{
int i,j,n=10,m=10; // n,m
float t,dd,d0,d1,dt=1.0/10.0; // 1/n
float tt,dtt=1.0/100.0; // 1/(n*m)
vec2 p0,p1,p2,p3,u,v;
vec2 q0,q1,q2,q3;
p0=gl_in[0].gl_Position.xy;
p1=gl_in[1].gl_Position.xy;
p2=gl_in[2].gl_Position.xy;
p3=gl_in[3].gl_Position.xy;
// compute BEZIER coefficients
a0.x= ( p0.x);
a1.x= (3.0*p1.x)-(3.0*p0.x);
a2.x= (3.0*p2.x)-(6.0*p1.x)+(3.0*p0.x);
a3.x=(p3.x)-(3.0*p2.x)+(3.0*p1.x)-( p0.x);
a0.y= ( p0.y);
a1.y= (3.0*p1.y)-(3.0*p0.y);
a2.y= (3.0*p2.y)-(6.0*p1.y)+(3.0*p0.y);
a3.y=(p3.y)-(3.0*p2.y)+(3.0*p1.y)-( p0.y);
q2=vec2(0.0,0.0);
q3=vec2(0.0,0.0);
// sample curve by chunks
for (p1=cubic(0.0),i=0,t=dt;i<n;i++,t+=dt)
{
// sample point
p0=p1; p1=cubic(t); q0=q2; q1=q3;
// compute ~OBB enlarged by D
u=normalize(p1-p0);
v=vec2(u.y,-u.x);
// resample chunk to compute enlargement
for (d0=0.0,d1=0.0,tt=t-dtt,j=2;j<m;j++,tt-=dtt)
{
dd=dot(cubic(tt)-p0,v);
d0=max(-dd,d0);
d1=max(+dd,d1);
}
d0+=d; d1+=d; u*=d;
d0*=1.25; d1*=1.25; // just to be sure
// enlarge radial
q2=p1+(v*d1);
q3=p1-(v*d0);
// enlarge axial
if (i==0)
{
q0=p0+(v*d1)-u;
q1=p0-(v*d0)-u;
}
if (i==n-1)
{
q2+=u;
q3+=u;
}
// pass it as QUAD
fcol=vcol[0]; trange=vec2(t-dt,t);
fpos=q0; gl_Position=vec4(q0,0.0,1.0); EmitVertex();
fpos=q1; gl_Position=vec4(q1,0.0,1.0); EmitVertex();
fpos=q2; gl_Position=vec4(q2,0.0,1.0); EmitVertex();
fpos=q3; gl_Position=vec4(q3,0.0,1.0); EmitVertex();
EndPrimitive();
}
}
//------------------------------------------------------------------------------*
Fragment:
// Fragment
#version 400 core
//#define show_coverage
uniform float d=0.05; // half thickness
in vec2 fpos; // fragment position
in vec3 fcol; // fragment color
in vec2 a0,a1,a2,a3; // cubic coefficients
in vec2 trange; // t range of chunk
out vec4 col;
vec2 cubic(float t) // return point on cubic from parameter
{
float tt=t*t,ttt=tt*t;
return a0+(a1*t)+(a2*tt)+(a3*ttt);
}
void main()
{
vec2 p;
int i,n;
float t,tt,t0,t1,dt,l,ll;
tt=-1.0; ll=-1.0; l=0.0;
#ifdef show_coverage
t0=0.0; t1=1.0; dt=0.05; n=3;
#else
t0=trange.x; n=2;
t1=trange.y;
dt=(t1-t0)*0.1;
#endif
for (i=0;i<n;i++)
{
for (t=t0;t<=t1;t+=dt)
{
p=cubic(t)-fpos;
l=length(p);
if ((ll<0.0)||(ll>l)){ ll=l; tt=t; }
}
t0=tt-dt; if (t0<0.0) t0=0.0;
t1=tt+dt; if (t1>1.0) t1=1.0;
dt*=0.2;
}
#ifdef show_coverage
if (ll>d) col=vec4(0.1,0.1,0.1,1.0); else
#else
if (ll>d) discard;
#endif
col=vec4(fcol,1.0);
}
And preview (curve + coverage):
And just curve:
as You can see the seam at the crossing wit hcoverage is due to coverage rendering without blending. The curve itself is OK.
The d0,d1 parameters are the max perpendicular distances to th actual chunk OBB axial axis (u) enlarged by d and scaled up by 25% just to be sure. Looks like it fits very good. I doubt there is much to be gained by further optimizations as this result is pretty close to perfect fit of the coverage...
the #define show_coverage just enables to view what geometry is passed to fragment shader ...
Solution 2:
See this tricky bezier curve:
There is no one solution for t, there are (up to) 3 solutions.
(edit1: As stated in Salix alba's answer, it doesn't mean you can't find them. When you thought that there was only one tx and one ty, you checked if they're (almost) equal. Going to 3 solutions, you could find the tx's and ty's and check if there is an (almost) common real value, but I think it should be sufficient (and faster) to check if bezierY(tx) is (almost) equal to glFragCoord.y for any tx, without calculating any ty. Also since tx's are the same for every pixel that has the same x, see if you can calculate them only once for each unique x.)
I haven't worked with bezier curves much, and never with glsl, so here's an idea that might be bad:
Every time your control points change, do a t loop to generate a list of {x,y} points, and possibly store them in some kind of unordered map. Then, in your shader, for every pixel, if that pixel exists in that map, apply the desired effect.
You can add nearby points too, and store the distance from the curve as the value in the map, so you can do some kind of anti-aliasing if you want to.
The step size in the t loop will have to be small enough, so that no points will be missed, but large enough, so that it will be fast. You can implement a dynamic t step, by checking how close the next point is to the previous point. If it's too close, increase the step. If it's too far, decrease the step.
You can also try using a 2d array instead of a map, something like 512x512 booleans. Initialize every element with false, and change the values to true for the points generated in your t loop. Meanwhile, store a list of the array indices that are currently true, so you can only initialize the 2d array once, and when your curve changes, flip every true back to false, empty your list of indices, and repeat the t loop etc.
(edit2, after your update)
Instead of searching "how to express the variable t from a cubic bezier equation", you can search for a "cubic equation solution" generally. If I'm not mistaken, the bezier equations (of x or y) can be written as
(-a + 3b - 3c + d) t^3 + (3a - 6b + 3c) t^2 + (-3a + 3b) t + (a - x) = 0
where a, b, c and d are the x (or y) components of the control points, and x is the x (or y) component of the curve, so they're just cubic equations. See that x appears only in the last coefficient, which might make things easier when you need to solve lots of them and their only difference is the value of x.
There should be simpler solutions, but if you have access to complex arithmetic (or are willing to write it yourself using vec2, see Spektre's answer an "How to compute Discrete Fourier Transform"), you can try these 3 solutions for t I got from Mathematica (I is the imaginary unit):
(-2*(a - 2*b + c) + (2*2^(1/3)*(b^2 + c^2 + a*(-c + d) - b*(c + d)))/(-2*b^3 + 3*a*b*c + 3*b^2*c - 6*a*c^2 + 3*b*c^2 - 2*c^3 - a^2*d + 3*a*b*d - 6*b^2*d + 3*a*c*d + 3*b*c*d - a*d^2 + a^2*x - 6*a*b*x + 9*b^2*x + 6*a*c*x - 18*b*c*x + 9*c^2*x - 2*a*d*x + 6*b*d*x - 6*c*d*x + d^2*x + Sqrt[(a - 3*b + 3*c - d)^2*(4*b^3*(d - x) + a^2*(d - x)^2 + x*(-4*c^3 + 9*c^2*x - 6*c*d*x + d^2*x) - 3*b^2*(c^2 - 2*c*x + (4*d - 3*x)*x) + 2*a*(2*c^3 - 6*c^2*x + 3*c*x*(d + x) - d*x*(d + x)) + 6*b*(a*(c - x)*(-d + x) + x*(c^2 + c*d - 3*c*x + d*x)))])^(1/3) + 2^(2/3)*(-2*b^3 + 3*a*b*c + 3*b^2*c - 6*a*c^2 + 3*b*c^2 - 2*c^3 - a^2*d + 3*a*b*d - 6*b^2*d + 3*a*c*d + 3*b*c*d - a*d^2 + a^2*x - 6*a*b*x + 9*b^2*x + 6*a*c*x - 18*b*c*x + 9*c^2*x - 2*a*d*x + 6*b*d*x - 6*c*d*x + d^2*x + Sqrt[(a - 3*b + 3*c - d)^2*(4*b^3*(d - x) + a^2*(d - x)^2 + x*(-4*c^3 + 9*c^2*x - 6*c*d*x + d^2*x) - 3*b^2*(c^2 - 2*c*x + (4*d - 3*x)*x) + 2*a*(2*c^3 - 6*c^2*x + 3*c*x*(d + x) - d*x*(d + x)) + 6*b*(a*(c - x)*(-d + x) + x*(c^2 + c*d - 3*c*x + d*x)))])^(1/3))/(2*(-a + 3*b - 3*c + d))
(-36*(a - 2*b + c) - ((18*I)*2^(1/3)*(-I + Sqrt[3])*(b^2 + c^2 + a*(-c + d) - b*(c + d)))/(-2*b^3 + 3*a*b*c + 3*b^2*c - 6*a*c^2 + 3*b*c^2 - 2*c^3 - a^2*d + 3*a*b*d - 6*b^2*d + 3*a*c*d + 3*b*c*d - a*d^2 + a^2*x - 6*a*b*x + 9*b^2*x + 6*a*c*x - 18*b*c*x + 9*c^2*x - 2*a*d*x + 6*b*d*x - 6*c*d*x + d^2*x + Sqrt[(a - 3*b + 3*c - d)^2*(4*b^3*(d - x) + a^2*(d - x)^2 + x*(-4*c^3 + 9*c^2*x - 6*c*d*x + d^2*x) - 3*b^2*(c^2 - 2*c*x + (4*d - 3*x)*x) + 2*a*(2*c^3 - 6*c^2*x + 3*c*x*(d + x) - d*x*(d + x)) + 6*b*(a*(c - x)*(-d + x) + x*(c^2 + c*d - 3*c*x + d*x)))])^(1/3) + (9*I)*2^(2/3)*(I + Sqrt[3])*(-2*b^3 + 3*a*b*c + 3*b^2*c - 6*a*c^2 + 3*b*c^2 - 2*c^3 - a^2*d + 3*a*b*d - 6*b^2*d + 3*a*c*d + 3*b*c*d - a*d^2 + a^2*x - 6*a*b*x + 9*b^2*x + 6*a*c*x - 18*b*c*x + 9*c^2*x - 2*a*d*x + 6*b*d*x - 6*c*d*x + d^2*x + Sqrt[(a - 3*b + 3*c - d)^2*(4*b^3*(d - x) + a^2*(d - x)^2 + x*(-4*c^3 + 9*c^2*x - 6*c*d*x + d^2*x) - 3*b^2*(c^2 - 2*c*x + (4*d - 3*x)*x) + 2*a*(2*c^3 - 6*c^2*x + 3*c*x*(d + x) - d*x*(d + x)) + 6*b*(a*(c - x)*(-d + x) + x*(c^2 + c*d - 3*c*x + d*x)))])^(1/3))/(36*(-a + 3*b - 3*c + d))
(-36*(a - 2*b + c) + ((18*I)*2^(1/3)*(I + Sqrt[3])*(b^2 + c^2 + a*(-c + d) - b*(c + d)))/(-2*b^3 + 3*a*b*c + 3*b^2*c - 6*a*c^2 + 3*b*c^2 - 2*c^3 - a^2*d + 3*a*b*d - 6*b^2*d + 3*a*c*d + 3*b*c*d - a*d^2 + a^2*x - 6*a*b*x + 9*b^2*x + 6*a*c*x - 18*b*c*x + 9*c^2*x - 2*a*d*x + 6*b*d*x - 6*c*d*x + d^2*x + Sqrt[(a - 3*b + 3*c - d)^2*(4*b^3*(d - x) + a^2*(d - x)^2 + x*(-4*c^3 + 9*c^2*x - 6*c*d*x + d^2*x) - 3*b^2*(c^2 - 2*c*x + (4*d - 3*x)*x) + 2*a*(2*c^3 - 6*c^2*x + 3*c*x*(d + x) - d*x*(d + x)) + 6*b*(a*(c - x)*(-d + x) + x*(c^2 + c*d - 3*c*x + d*x)))])^(1/3) - 9*2^(2/3)*(1 + I*Sqrt[3])*(-2*b^3 + 3*a*b*c + 3*b^2*c - 6*a*c^2 + 3*b*c^2 - 2*c^3 - a^2*d + 3*a*b*d - 6*b^2*d + 3*a*c*d + 3*b*c*d - a*d^2 + a^2*x - 6*a*b*x + 9*b^2*x + 6*a*c*x - 18*b*c*x + 9*c^2*x - 2*a*d*x + 6*b*d*x - 6*c*d*x + d^2*x + Sqrt[(a - 3*b + 3*c - d)^2*(4*b^3*(d - x) + a^2*(d - x)^2 + x*(-4*c^3 + 9*c^2*x - 6*c*d*x + d^2*x) - 3*b^2*(c^2 - 2*c*x + (4*d - 3*x)*x) + 2*a*(2*c^3 - 6*c^2*x + 3*c*x*(d + x) - d*x*(d + x)) + 6*b*(a*(c - x)*(-d + x) + x*(c^2 + c*d - 3*c*x + d*x)))])^(1/3))/(36*(-a + 3*b - 3*c + d))
They are large, but they contain many common sub-expressions (like (a - 2*b + c)) which you can evaluate once and reuse, to improve performance (if all this works at all).
For the tricky bezier I posted, here are the 3 solutions:
red = (6 + (4*2^(1/3))/(-9 + 49*x + 7*Sqrt[1 + x*(-18 + 49*x)])^(1/3) + 2^(2/3)*(-9 + 49*x + 7*Sqrt[1 + x*(-18 + 49*x)])^(1/3))/14
green = (12 - ((4*I)*2^(1/3)*(-I + Sqrt[3]))/(-9 + 49*x + 7*Sqrt[1 - 18*x + 49*x^2])^(1/3) + I*2^(2/3)*(I + Sqrt[3])*(-9 + 49*x + 7*Sqrt[1 - 18*x + 49*x^2])^(1/3))/28
blue = (12 + ((4*I)*2^(1/3)*(I + Sqrt[3]))/(-9 + 49*x + 7*Sqrt[1 - 18*x + 49*x^2])^(1/3) - 2^(2/3)*(1 + I*Sqrt[3])*(-9 + 49*x + 7*Sqrt[1 - 18*x + 49*x^2])^(1/3))/28
(edit3) Following Spektre's suggestion, using the coefficients of a cubic directly
x = a*t^3 + b*t^2 + c*t + d
(instead of using the control points' coordinates) gives cleaner expressions:
1st(red) = (-2*b + (2*2^(1/3)*(b^2 - 3*a*c))/(-2*b^3 + 9*a*b*c - 27*a^2*d + Sqrt[-4*(b^2 - 3*a*c)^3 + (2*b^3 - 9*a*b*c + 27*a^2*(d - x))^2] + 27*a^2*x)^(1/3) + 2^(2/3)*(-2*b^3 + 9*a*b*c - 27*a^2*d + Sqrt[-4*(b^2 - 3*a*c)^3 + (2*b^3 - 9*a*b*c + 27*a^2*(d - x))^2] + 27*a^2*x)^(1/3))/(6*a)
2nd(green) = (-4*b + (2*2^(1/3)*(1 + I*Sqrt[3])*(-b^2 + 3*a*c))/(-2*b^3 + 9*a*b*c - 27*a^2*d + Sqrt[-4*(b^2 - 3*a*c)^3 + (2*b^3 - 9*a*b*c + 27*a^2*(d - x))^2] + 27*a^2*x)^(1/3) + I*2^(2/3)*(I + Sqrt[3])*(-2*b^3 + 9*a*b*c - 27*a^2*d + Sqrt[-4*(b^2 - 3*a*c)^3 + (2*b^3 - 9*a*b*c + 27*a^2*(d - x))^2] + 27*a^2*x)^(1/3))/(12*a)
3rd(blue) = -(4*b - ((2*I)*2^(1/3)*(I + Sqrt[3])*(b^2 - 3*a*c))/(-2*b^3 + 9*a*b*c - 27*a^2*d + Sqrt[-4*(b^2 - 3*a*c)^3 + (2*b^3 - 9*a*b*c + 27*a^2*(d - x))^2] + 27*a^2*x)^(1/3) + 2^(2/3)*(1 + I*Sqrt[3])*(-2*b^3 + 9*a*b*c - 27*a^2*d + Sqrt[-4*(b^2 - 3*a*c)^3 + (2*b^3 - 9*a*b*c + 27*a^2*(d - x))^2] + 27*a^2*x)^(1/3))/(12*a)
You can easily transform your control point's coordinates into these coordinates:
direct a = control (-a + 3 b - 3 c + d)
direct b = control (3 a - 6 b + 3 c)
direct c = control (-3 a + 3 b)
direct d = control a