Signed distance to quadratic bezier
category: code [glöplog]
Having limited math skills I cant get my head around calculating the distance from a 2d point to a 2d quadratic curve analytically.
But I have a gut feeling this might prove to the fastest solution,
even so its not a direct solution.
http://research.microsoft.com/en-us/um/people/cloop/loopblinn05.pdf
The antialiasing method give a signed distance...
So is the solution simply to get the baricentric coordinates in relation to the 3 control point and evaluating the signed distance as described in this pixel shader?
http://http.developer.nvidia.com/GPUGems3/gpugems3_ch25.html
float d = SignedDistance(float px, float py, float sx[3], float sy[3]) {
....
}
My goal is to draw thick & soft edge curve without geometry.
But I have a gut feeling this might prove to the fastest solution,
even so its not a direct solution.
http://research.microsoft.com/en-us/um/people/cloop/loopblinn05.pdf
The antialiasing method give a signed distance...
So is the solution simply to get the baricentric coordinates in relation to the 3 control point and evaluating the signed distance as described in this pixel shader?
http://http.developer.nvidia.com/GPUGems3/gpugems3_ch25.html
float d = SignedDistance(float px, float py, float sx[3], float sy[3]) {
....
}
My goal is to draw thick & soft edge curve without geometry.
Just as a suggestion, you could find the closest point on the curve to your shading point and subtract your stroke radius from the euclidean distance to that point. This has the benefit of giving a meaningful distance to the stroke's outline.
unlimited bezier curves geometry edges.
If your quadratic bezier is given parametrically as (at^2+bt+c, dt^2+et+f), then the tangent vector at the point for a particular t is the derivative of the parameterization, i.e. (2at+b, 2dt+e). At the point on the curve closest to some arbitrary point (x, y), the line from (x, y) to the point on the curve is perpendicular to the tangent, so (at^2+bt+c - x)*(2at+b) + (dt^2+et+f - y)*(2dt+e) = 0.
This is a cubic equation in t. It can have 1, 2 or 3 solutions, some of which might not correspond to the minimum distance. Wikipedia gives some suggestions on how to find all solutions.
This is a cubic equation in t. It can have 1, 2 or 3 solutions, some of which might not correspond to the minimum distance. Wikipedia gives some suggestions on how to find all solutions.
i wonder if T21 got what you said ;)
i didnt, but i didnt even try to understand what you said yet. (i would if i would want to achieve that effect maybe ;) )
why i am posting is: i guess using some sin/cos-ed(position) unlimited_in_height_Cylinder does the same job, with lesser brainPower used and also less amount of code. yielding in less amount of bytes in the final prod due to your crinkler being so good at compressing ascii ;)
i didnt, but i didnt even try to understand what you said yet. (i would if i would want to achieve that effect maybe ;) )
why i am posting is: i guess using some sin/cos-ed(position) unlimited_in_height_Cylinder does the same job, with lesser brainPower used and also less amount of code. yielding in less amount of bytes in the final prod due to your crinkler being so good at compressing ascii ;)
fizzer: Thats the goal. I have this working for a few primitives. Just not quatradic curve.
blue: I seen a few solution, but they all look massive and not compute friendly. Like this one. the thirdDegreeEquation evaluation look painfull
http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html
I just dont have the skill to even state : "This cant be simplified"
If this is the most optical method to analytically find a distance to a quadratic curve, then I cant use that and need to come up with a hack.
Hence trying to derive from the trick of rendering a curve via interpolated UV of a bounding triangle of the bezier curve.
hardy: What I find wild is how complex the math is for solving this, yet if you use barricentric coordinate of the bounding triangle it seem trivial.
Or am I missing something in my assumption.
I have to wonder, because I dont need the actual point on the curve, just the minimum distance. Could this be leveraged for some simplifications?
blue: I seen a few solution, but they all look massive and not compute friendly. Like this one. the thirdDegreeEquation evaluation look painfull
http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html
I just dont have the skill to even state : "This cant be simplified"
If this is the most optical method to analytically find a distance to a quadratic curve, then I cant use that and need to come up with a hack.
Hence trying to derive from the trick of rendering a curve via interpolated UV of a bounding triangle of the bezier curve.
hardy: What I find wild is how complex the math is for solving this, yet if you use barricentric coordinate of the bounding triangle it seem trivial.
Or am I missing something in my assumption.
I have to wonder, because I dont need the actual point on the curve, just the minimum distance. Could this be leveraged for some simplifications?
What Blueberry said. But if you want something very simple then you could treat the curve as a series of line segments and take the minimum of the distances to those segments. The suitability of this approach depends on your application.
What Blueberry outlined would work, but its just to slow if it cant be optimized further then the .as code I linked to.
But before doing an approximation using subdivision.
Wouldn't leveraging this be simpler then what blueberry outlined ? http://http.developer.nvidia.com/GPUGems3/gpugems3_ch25.html
Because I dont have geometry, I would generate the UV over the bounding triangle, made from the 3 control point, using a barricentric formula then just do d = (u*u-v) scale by the gradient to get the signed distance ?
What I'm not sure is what happen for UV outside 0-1 (point out of the bounding triangle) would this still generated correct distances?
But before doing an approximation using subdivision.
Wouldn't leveraging this be simpler then what blueberry outlined ? http://http.developer.nvidia.com/GPUGems3/gpugems3_ch25.html
Because I dont have geometry, I would generate the UV over the bounding triangle, made from the 3 control point, using a barricentric formula then just do d = (u*u-v) scale by the gradient to get the signed distance ?
What I'm not sure is what happen for UV outside 0-1 (point out of the bounding triangle) would this still generated correct distances?
as Blueberry said, it's a cubic thing:
Your 3d curve is f(t) = at² +bt + c (a, b and c are 3d vectors, no need to work with coordinates). Your point is x (a vector too, no need for x, y and z). The distance squared |f(t) - x|² is
|f(t) - x|² = |at² +bt + c - x|² = a²t^4 + 2abt^3 + 2a(c-x)t^2 + b²t² + 2b(c-x)t + |c-x|²
(you can square and multiply vectors, it's called "dot product"!)
your local minima is when
2a²t^3 + 3abt² + (2a(c-x)+b²)t + b(c-x) = 0
which is a cubic as Blueberry pointed. You can solve the cubic with this piece of code (i've been using it for a while now):
which needs you to normalize the cubic so the coefficient of the cubic term is 1, ie:
meaning,
Of course, make sure your root is between 0 and 1
disclaimer - i wrote these maths directly here in the post box, and it probably contains erros
Your 3d curve is f(t) = at² +bt + c (a, b and c are 3d vectors, no need to work with coordinates). Your point is x (a vector too, no need for x, y and z). The distance squared |f(t) - x|² is
|f(t) - x|² = |at² +bt + c - x|² = a²t^4 + 2abt^3 + 2a(c-x)t^2 + b²t² + 2b(c-x)t + |c-x|²
(you can square and multiply vectors, it's called "dot product"!)
your local minima is when
2a²t^3 + 3abt² + (2a(c-x)+b²)t + b(c-x) = 0
which is a cubic as Blueberry pointed. You can solve the cubic with this piece of code (i've been using it for a while now):
Code:
#include <math.h>
static __inline float cuberoot( float x )
{
if( x<0.0f ) return -powf(-x,1.0f/3.0f);
return powf(x,1.0f/3.0f);
}
// Calculates the 3 roots of the cubic. Returns the number of real roots (1 or 3)
// 1 real root:
// res[0] = real root
// res[2] = Re{root 2}
// res[3] = Im{root 2}
// res[4] = Re{root 3}
// res[5] = Im{root 3}
// 3 real roots
// res[0] = real root 1
// res[2] = real root 2
// res[4] = real root 3
int solveCubic( float b, float c, float d, float *res )
{
// x^3 + a2·x^2 + a1·x + a0 = 0
float p = (3.0f*c-b*b)/3.0f;
float q = (9.0f*b*c - 27.0f*d - 2.0f*b*b*b)/27.0f;
// x^3 + px + q = 0
// x = w-(p/3w)
// (w^3)^2 - q*(w^3) - (p^3)/27 = 0
float h = q*q*.25f + p*p*p/27.0f;
if( h>=0.0f )
{
// one single real solution!
h = sqrtf(h);
float r = (q*.5f) + h;
float t = (q*.5f) - h;
float s = cuberoot( r );
float u = cuberoot( t );
res[2*0+0] = (s + u) - (b/3.0f);
res[2*0+0] = 0.0f;
float re = -(s + u)*.5f - (b/3.0f);
float im = (s-u)*(sqrtf(3.0f)*.5f);
res[2*1+0] = re;
res[2*1+1] = im;
res[2*2+0] = re;
res[2*2+1] = -im;
return 1;
}
else
{
// three real solutions!
const float sqrt3 = sqrtf(3.0f);
const float ovsqrt27 = 1.0f/sqrtf(27.0f);
const float ov3 = 1.0f/3.0f;
float i = p*sqrtf( -p )*ovsqrt27;
float j = cuberoot( i );
float k = ov3*acosf( (q/(2.0f*i)) );
float m = cosf(k);
float n = sinf(k)*sqrt3;
float s = -b*ov3;
res[2*0+0] = 2.0f*j*m + s;
res[2*1+0] = -j*(m + n) + s;
res[2*2+0] = -j*(m - n) + s;
return 3;
}
return 1;
}
which needs you to normalize the cubic so the coefficient of the cubic term is 1, ie:
meaning,
Code:
int nroots = solveCubic( 3ab/2a², a(c-x)+b²)/a², b(c-x)2a²/, results );
Of course, make sure your root is between 0 and 1
disclaimer - i wrote these maths directly here in the post box, and it probably contains erros
fixes:
and
Code:
res[2*0+1] = 0.0f;and
Code:
int nroots = solveCubic( 3ab/2a², (a(c-x)+b²)/a², b(c-x)/2a², results );Thats the cleanest code I seen for solving the roots.
So doing this is the full solution ??
(I think I'm missing something big:)
So doing this is the full solution ??
(I think I'm missing something big:)
Code:
x = point
a = curve start point
b = curve control control point
c = curve end point
float DistanceSquared(vector x, vector a, vector b, vector c) {
float results[6];
int nroots = solveCubic(
dot(a*3,b) / dot(a*2,a*2),
(dot(a, (c-x)) + dot(b,b)) / dot(a,a),
dot(b, (c-x)) / (dot(a*2,a*2)),
results);
float minD = 10e10;
for(int i=0; i<nroots; i++) {
float t = results[2*i+0];
if(t>=0 && t<=1) {
vector s = Eval(a, b, c, t);
vector d = s-x;
minD = dot(d,d);
}
}
return minD;
}
no. a, b and c are not the curve control points. they are a combination od the control points. you'll have to compute them. for example, if your control points are P0 as curve start, P2 as curve end, and P1-P0 as a tangent at t=0, you can do
f(0) = a
f(1) = a+b+c = P2
f'(0) = b = P1-P0
and then you can solve a, b and c, and proceed from there. You can also set a different set of condutions, like P0 is curve start, P1 is curve at t=0.5 and P2 is curve end.
You can invent the conditions tha suit your problem best. For quadratic Bezier curve, I don't know the conditions my heart.
f(0) = a
f(1) = a+b+c = P2
f'(0) = b = P1-P0
and then you can solve a, b and c, and proceed from there. You can also set a different set of condutions, like P0 is curve start, P1 is curve at t=0.5 and P2 is curve end.
You can invent the conditions tha suit your problem best. For quadratic Bezier curve, I don't know the conditions my heart.
ah, here we go, indeed
f(t) = (1-t)²·P0 + 2(1-t)t·P1 + t²·P2
where indeed f(0) = P0 and f(1)=P2. now, if you expand this you get (and again, i'm doing this here, so you better double check it)
a = P0 - 2·P1 + P2
b = 2·(P0 + P1)
c = P0
f(t) = (1-t)²·P0 + 2(1-t)t·P1 + t²·P2
where indeed f(0) = P0 and f(1)=P2. now, if you expand this you get (and again, i'm doing this here, so you better double check it)
a = P0 - 2·P1 + P2
b = 2·(P0 + P1)
c = P0
so, something like this???
Code:
float distSquaredToQuadratocBezier( in vec3 P0, in vec3 P1, in vec3 P2, in vec3 x )
{
// f(t) = (1-t)²·P0 + 2t(1-t)·P1 + t²·P2
vec3 a = P0 - 2.0f*P1 + P2
vec3 b = 2.0f*(P0 + P1)
vec3 c = P0
float k3 = 2.0f*dot(a,a)
float k2 = 3.0f*dot(a,b)
float k1 = dot(a,c-x) + dot(b,b)
float k0 = dot(b,c-x)
float res[6];
int n = solveCubic( k2/k3, k1/k3, k0/k3, res );
float minDis = 1e20f;
for( int i=0; i<n; i++ )
{
float t = res[2*i+0];
if( t<0.0f ) t=0.0f; else if( t>1.0f ) t=1.0f;
vec3 p = (1.0f-t)*(1.0f-t)*P0 + 2.0f*t*(1.0f-t)*P1 + t*t*P2;
float dis = dot(p-x,p-x);
if( dis<minDis ) minDis = dis;
}
return minDis;
}why bother even calculating the complex roots? if it's not real, it's not in [0,1] and thus not interesting as a solution to this problem.
The code as is seem to always return the distance to P0,
and I cant verify the correctness because the math is beyond me.
I assume the solveCubic() is a time tested function (I cant imagine anyone typing this freestyle in a BBS edit box :), so I focus on the generation of the parameters.
Only one thing got my attention so far
is not identical mathematically to
and I cant verify the correctness because the math is beyond me.
I assume the solveCubic() is a time tested function (I cant imagine anyone typing this freestyle in a BBS edit box :), so I focus on the generation of the parameters.
Only one thing got my attention so far
Code:
int nroots = solveCubic( 3ab/2a², (a(c-x)+b²)/>>a²<<, b(c-x)/2a², results );
is not identical mathematically to
Code:
float k3 = 2.0f*dot(a,a)
int n = solveCubic( k2/k3, k1/>>k3<<, k0/k3, res );
side note, I tried to derive a distance from the anti aliasing trick in the papper I linked. And I got this (draw 4 curves to build a circle)
a) Get UV coordinates of triangle formed by the 3 control point
b) Use pixel shader code to get distance
Its not a complete failure, but the distance is incorrect :)
a) Get UV coordinates of triangle formed by the 3 control point
b) Use pixel shader code to get distance
Its not a complete failure, but the distance is incorrect :)
I think this is the bezier to quadratic
And if you have those terms, you can compute p like this
The result is still not correct... but now I'm at a lose as to whats is going on now.
Code:
// f(t) = (1-t)²·P0 + 2t(1-t)·P1 + t²·P2
vec3 a = P0
vec3 b = (P1 - P0) * 2.0f (and not (P0-P1)*2)
vec3 c = P0 - 2.0f*P1 + P2
And if you have those terms, you can compute p like this
Code:
vec3 p = a + (b + c * t) * t;
vs
vec3 p = (1.0f-t)*(1.0f-t)*P0 + 2.0f*t*(1.0f-t)*P1 + t*t*P2;
The result is still not correct... but now I'm at a lose as to whats is going on now.
I couldn't find the typo in the posted formula, ever after quadruple checking against posted solution like .
http://www.1728.org/cubic2.htm
But I found some other sample code
http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html
And that worked for me.
And now the fun part begin (optimizations)
Mostly take the constant out and make them part of the quadratic curve structure.
But first on my list, find a fast cuberoot() approximation... anyone ? :)
Also here is the solver
http://www.1728.org/cubic2.htm
But I found some other sample code
http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html
And that worked for me.
And now the fun part begin (optimizations)
Mostly take the constant out and make them part of the quadratic curve structure.
But first on my list, find a fast cuberoot() approximation... anyone ? :)
Also here is the solver
Code:
int thirdDegreeEquation(float a, float b, float c, float* results)
{
float p = b - a*a / 3;
float q = a * (2*a*a - 9*b) / 27 + c;
float p3 = p*p*p;
float D = q*q + 4*p3 / 27;
float offset = -a / 3;
if(D >= 0) {
float z = sqrtf(D);
float u = ( -q + z) / 2;
float v = ( -q - z) / 2;
u = cuberoot(u);
v = cuberoot(v);
results[0] = u + v + offset;
return 1;
} else {
float u = 2 * sqrtf( -p / 3);
float v = acosf( -sqrtf( -27 / p3) * q / 2) / 3;
results[0] = u * cosf(v) + offset;
results[1] = u * cosf(v + 2 * 3.1416 / 3) + offset;
results[2] = u * cosf(v + 4 * 3.1416 / 3) + offset;
return 3;
}
}
I'm happy, I got my thick & soft quadratic bezier curves done now via SDF
note: the control point in the pic are also SDF primitives.
note: the control point in the pic are also SDF primitives.
yeah, my code had a few bugs (pouet coding!), among them the one you found for b = 2.0f*(P1-P0) and the other one being k1 = dot(a,c-x) + dot(b,b) ^_^
i mean, float k1 = 2.0f*dot(a,c-x) + dot(b,b)
sdf saved just another soul ;D
Thx iq, here is the final working code for reference.
(Drew a few True type font character outlines to test)
BTW, this code, look so much shorter then anything I came across.
(Drew a few True type font character outlines to test)
BTW, this code, look so much shorter then anything I came across.
Code:
int solveCubic(float a, float b, float c, float* r)
{
float p = b - a*a / 3;
float q = a * (2*a*a - 9*b) / 27 + c;
float p3 = p*p*p;
float d = q*q + 4*p3 / 27;
float offset = -a / 3;
if(d >= 0) { // Single solution
float z = sqrtf(d);
float u = (-q + z) / 2;
float v = (-q - z) / 2;
u = cuberoot(u);
v = cuberoot(v);
r[0] = offset + u + v;
return 1;
}
float u = sqrtf(-p / 3);
float v = acos(-sqrtf( -27 / p3) * q / 2) / 3;
float m = cos(v), n = sin(v)*1.732050808f;
r[0] = offset + u * (m + m);
r[1] = offset - u * (n + m);
r[2] = offset + u * (n - m);
return 3;
}
Code:
float SignedDistanceSquared(point2D p, QBSpline2D s)
{
float minDis = 1e20f;
float res[3];
point2D D = s.p0 - p;
float a = s.A;
float b = s.B;
float c = s.C + dot(D, s.c);
float d = dot(D, s.d);
int n = solveCubic(b*a, c*a, d*a, res);
for(int j=0; j<n; j++) {
float t = Clamp(res[j]);
point2D d = s.p0 + (s.b + s.c*t)*t - p;
minDis = min(minDis, dot(d,d));
}
return minDis;
}
Unless anyone wants to talk optimization etc.. I can close this thread with this image showing the stepping used to render outlines.
The green show when and only when the distance function is evaluated. Gray is just no.op
It is not much different then sphere tracing, but in 2D.
x += max(1, d-edge); // Skip empty space
instead of the naive
x += 1
The green show when and only when the distance function is evaluated. Gray is just no.op
It is not much different then sphere tracing, but in 2D.
x += max(1, d-edge); // Skip empty space
instead of the naive
x += 1