Menu
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=)
311
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:
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:
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 ...
114
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
