Fast maximum curvature calculation for a 3d bezier curve - Graphics

"sam" <asimnaseer@gmail.com> wrote in message
news:1157597215.061285.10380@p79g2000cwp.googlegroups.com...

> I have been looking at efficient ways of calculating the maximum
> curvature of a bezier curve,

What is the degree of your curves?

--
Dave Eberly
http://www.geometrictools.com

Dave Eberly wrote:
> "sam" <asimnaseer@gmail.com> wrote in message
> news:1157597215.061285.10380@p79g2000cwp.googlegroups.com...
>
> > I have been looking at efficient ways of calculating the maximum
> > curvature of a bezier curve,
>
> What is the degree of your curves?
>
> --
> Dave Eberly
> http://www.geometrictools.com

Its a cubic curve. Sorry, should have mentioned.

Sam.

"sam" <asimnaseer@gmail.com> wrote in message
news:1157608649.309785.238610@i42g2000cwa.googlegroups.com...
>
> Its a cubic curve. Sorry, should have mentioned.

Thanks for the clarification.

Sometimes people use the term "curvature" to mean
magnitude of the second-derivative vector. They are
in fact different, so I am assuming you really mean
"curvature" in the differential geometric sense.

=====
About the subdivision of the cubic curve using a test
for flatness. The typical fast subdivision involves
a stopping criterion involving the magnitude of the
second-derivative vector and the length of the current
subinterval.
void Subdivide (
PointList L,
float t0, float t1,
Point x0, Point x1,
Point sd0, Point sd1)
{
sdmid = 0.5*(sd0 + sd1);
dt = t1 - t0;
nonlinearity = dt*dt*sdmid;
if (SquaredLength(nonlinearity) > epsilon)
{
tmid = 0.5*(t0 + t1);
xmid = 0.5*(x0 + x1 - nonlinearity);
insert xmid in L between x0 and x1;
Subdivide(t0,tmid,x0,xmid,sd0,sdmid);
Subdivide(tmid,t1,xmid,x1,smid,sd1);
}
}
where t is the curve parameter, x is position, and
sd is second derivative of position. The epsilon
is a user-specified tolerance The subdivision is
initiated with
x0 = P(0); x1 = P(1); sd0 = P"(0); sd1 = P"(1);
L = {x0,x1};
Subdivide(L,0,1,x0,x1,sd0,sd1);
The following pathological example (in 2D) shows how
the subdivision approach can grossly underestimate
the maximum curvature. Naturally, your data might
not lead to such problems.

The control points are (0,0), (0,2), (e,2), and (e,1)
for a small positive e. The Bezier curve is
p"(t) = (x(t),y(t))
with
x(t) = e*(3*(1-t)*t^2 + t^3) = e*u(t)
y(t) = 6*(1-t)^2*t + 6*(1-t)*t^2 + t^3
for 0 <= t <= 1. The last equality in the x(t) equation
defines the polynomial u(t). The curvature is
k(t) = [x'(t)*y"(t) - x"(t)*y'(t)] / s(t)^3
where
s(t) = sqrt(x'(t)^2 + y'(t)^2)
Notice that y'(r) = 0 for r = sqrt(2)/(sqrt(2)+1). It
is easily shown that y"(r) and u'(r) are both not zero.
The curvature at this parameter is
k(r) = |y"(r)|/(e*u'(r))^2
As you make e smaller, k(r) becomes very large. Now at
any other value of t, y'(t) is not zero. As you make
e smaller, the curvature k(t) becomes smaller. When
you look at the graph of the curve, it is essentially
two very flat subcurves connected by a tiny curve with
very large curvature.

The second derivative of the curve is
p"(t) = 6*(e*(1-2t),t-2)
For small e, the x-component is nearly zero, so the
stopping criterion for the subdivision will effectively
be controlled by the length of the subintervals. When
you evaluate the curvature at the subdivision points,
you can always choose e small enough so that those
curvatures are nearly zero, in which case you would
conclude that the curve is nearly flat everywhere. The
reality is that the curve has a tiny region of very
large curvature.

=====
Knowing now that such a subdivision scheme can incorrectly
report the maximum curvature, I will suggest another
subdivision scheme

For a 3D curve p(t), the curvature is
k(t) = |Cross(p'(t),p"(t))|/|p'(t)|^3
The squared curvature is a ratio of two polynomials in t.
An ****ytical approach to computing the maximum curvature
is to compute the critical points for k(t)^2. From calculus,
These are the t-values where the derivative is zero, and
at the endpoints t = 0 and t = 1. The numerator of the
derivative of k(t)^2 is
n(t) = p11*(p11*p23 - p12*p13) - 3*p12*(p11*p22-p12^2)
where p1(t) = p'(t), p2(t) = p"(t), and pij = Dot(pi,pj).
If you just count the degrees of the various polynomial
terms in n(t), you get degree 9. It turns out that the
degree is actually 7. So to compute the critical points,
you must compute the roots of a degree 7 polynomial. This
sounds heavy handed, but you could take advantage of root
bounding schemes to get _approximations_ to the roots and
use these to estimate the maximum curvature. You could
also compute the Sturm sequence of polynomials for n(t)
and determine if n(t) has any roots on [0,1]. If it does
not, then the maximum curvature must occur at t = 0 or
t = 1.

Consider the same ****ysis for a quadratic Bezier curve.
The numerator n(t) of the derivative of k(t)^2 happens
to be of degree 1. This is easy enough to solve for t,
so you have only three critical points to consider, the
root for n(t) = 0 and t-values 0 and 1. The short of
it is that it is very easy to compute the maximum
curvature for a quadratic Bezier curve.

Rather than subdivide your cubic Bezier curve into a
sequence of line segments (linear Bezier curves as it
were), subdivide it into a sequence of quadratic Bezier
curves. Each quadratic curve will hopefully be a good
estimate of the cubic subcurve it is associated with.

A simple way to subdivide is to half the parameter
intervals (as in the linear case). Use a least-squares
method to fit the cubic curve
p(t) = (1-t)^3*p0 + 3*(1-t)^2*t*p1 +
3*(1-t)*t^2*p2 + t^3*p3
with a quadratic curve
q(t) = (1-t)^2*q0 +2*(1-t)*t*q1 + t^2*q2
with q0 = p0 and q2 = p3. Use an integral norm to
minimize the integral of squared distances between
the two curves. That is,
E(q1) = integral[0,1] |q(t) - p(t)|^2 dt
When you work through all the math, the unknown control
point q1 must be
q1 = (-p0 + 3*p1 - 3*p2 + p3)/4
Whereas your stopping criterion for the line segments
relied on magnitude of second derivatives, the stopping
criterion for the quadratic curves can rely on the
value of E(q1).

--
Dave Eberly
http://www.geometrictools.com

Thank you for your reply Dave. I am going to try that out today with
some test cases, and will let you know how it goes.

Sam.

Dave Eberly wrote:
> "sam" <asimnaseer@gmail.com> wrote in message
> news:1157608649.309785.238610@i42g2000cwa.googlegroups.com...
> >
> > Its a cubic curve. Sorry, should have mentioned.
>
> Thanks for the clarification.
>
> Sometimes people use the term "curvature" to mean
> magnitude of the second-derivative vector. They are
> in fact different, so I am assuming you really mean
> "curvature" in the differential geometric sense.
>
> =====
> About the subdivision of the cubic curve using a test
> for flatness. The typical fast subdivision involves
> a stopping criterion involving the magnitude of the
> second-derivative vector and the length of the current
> subinterval.
> void Subdivide (
> PointList L,
> float t0, float t1,
> Point x0, Point x1,
> Point sd0, Point sd1)
> {
> sdmid = 0.5*(sd0 + sd1);
> dt = t1 - t0;
> nonlinearity = dt*dt*sdmid;
> if (SquaredLength(nonlinearity) > epsilon)
> {
> tmid = 0.5*(t0 + t1);
> xmid = 0.5*(x0 + x1 - nonlinearity);
> insert xmid in L between x0 and x1;
> Subdivide(t0,tmid,x0,xmid,sd0,sdmid);
> Subdivide(tmid,t1,xmid,x1,smid,sd1);
> }
> }
> where t is the curve parameter, x is position, and
> sd is second derivative of position. The epsilon
> is a user-specified tolerance The subdivision is
> initiated with
> x0 = P(0); x1 = P(1); sd0 = P"(0); sd1 = P"(1);
> L = {x0,x1};
> Subdivide(L,0,1,x0,x1,sd0,sd1);
> The following pathological example (in 2D) shows how
> the subdivision approach can grossly underestimate
> the maximum curvature. Naturally, your data might
> not lead to such problems.
>
> The control points are (0,0), (0,2), (e,2), and (e,1)
> for a small positive e. The Bezier curve is
> p"(t) = (x(t),y(t))
> with
> x(t) = e*(3*(1-t)*t^2 + t^3) = e*u(t)
> y(t) = 6*(1-t)^2*t + 6*(1-t)*t^2 + t^3
> for 0 <= t <= 1. The last equality in the x(t) equation
> defines the polynomial u(t). The curvature is
> k(t) = [x'(t)*y"(t) - x"(t)*y'(t)] / s(t)^3
> where
> s(t) = sqrt(x'(t)^2 + y'(t)^2)
> Notice that y'(r) = 0 for r = sqrt(2)/(sqrt(2)+1). It
> is easily shown that y"(r) and u'(r) are both not zero.
> The curvature at this parameter is
> k(r) = |y"(r)|/(e*u'(r))^2
> As you make e smaller, k(r) becomes very large. Now at
> any other value of t, y'(t) is not zero. As you make
> e smaller, the curvature k(t) becomes smaller. When
> you look at the graph of the curve, it is essentially
> two very flat subcurves connected by a tiny curve with
> very large curvature.
>
> The second derivative of the curve is
> p"(t) = 6*(e*(1-2t),t-2)
> For small e, the x-component is nearly zero, so the
> stopping criterion for the subdivision will effectively
> be controlled by the length of the subintervals. When
> you evaluate the curvature at the subdivision points,
> you can always choose e small enough so that those
> curvatures are nearly zero, in which case you would
> conclude that the curve is nearly flat everywhere. The
> reality is that the curve has a tiny region of very
> large curvature.
>
> =====
> Knowing now that such a subdivision scheme can incorrectly
> report the maximum curvature, I will suggest another
> subdivision scheme
>
> For a 3D curve p(t), the curvature is
> k(t) = |Cross(p'(t),p"(t))|/|p'(t)|^3
> The squared curvature is a ratio of two polynomials in t.
> An ****ytical approach to computing the maximum curvature
> is to compute the critical points for k(t)^2. From calculus,
> These are the t-values where the derivative is zero, and
> at the endpoints t = 0 and t = 1. The numerator of the
> derivative of k(t)^2 is
> n(t) = p11*(p11*p23 - p12*p13) - 3*p12*(p11*p22-p12^2)
> where p1(t) = p'(t), p2(t) = p"(t), and pij = Dot(pi,pj).
> If you just count the degrees of the various polynomial
> terms in n(t), you get degree 9. It turns out that the
> degree is actually 7. So to compute the critical points,
> you must compute the roots of a degree 7 polynomial. This
> sounds heavy handed, but you could take advantage of root
> bounding schemes to get _approximations_ to the roots and
> use these to estimate the maximum curvature. You could
> also compute the Sturm sequence of polynomials for n(t)
> and determine if n(t) has any roots on [0,1]. If it does
> not, then the maximum curvature must occur at t = 0 or
> t = 1.
>
> Consider the same ****ysis for a quadratic Bezier curve.
> The numerator n(t) of the derivative of k(t)^2 happens
> to be of degree 1. This is easy enough to solve for t,
> so you have only three critical points to consider, the
> root for n(t) = 0 and t-values 0 and 1. The short of
> it is that it is very easy to compute the maximum
> curvature for a quadratic Bezier curve.
>
> Rather than subdivide your cubic Bezier curve into a
> sequence of line segments (linear Bezier curves as it
> were), subdivide it into a sequence of quadratic Bezier
> curves. Each quadratic curve will hopefully be a good
> estimate of the cubic subcurve it is associated with.
>
> A simple way to subdivide is to half the parameter
> intervals (as in the linear case). Use a least-squares
> method to fit the cubic curve
> p(t) = (1-t)^3*p0 + 3*(1-t)^2*t*p1 +
> 3*(1-t)*t^2*p2 + t^3*p3
> with a quadratic curve
> q(t) = (1-t)^2*q0 +2*(1-t)*t*q1 + t^2*q2
> with q0 = p0 and q2 = p3. Use an integral norm to
> minimize the integral of squared distances between
> the two curves. That is,
> E(q1) = integral[0,1] |q(t) - p(t)|^2 dt
> When you work through all the math, the unknown control
> point q1 must be
> q1 = (-p0 + 3*p1 - 3*p2 + p3)/4
> Whereas your stopping criterion for the line segments
> relied on magnitude of second derivatives, the stopping
> criterion for the quadratic curves can rely on the
> value of E(q1).
>
> --
> Dave Eberly
> http://www.geometrictools.com

sam wrote:
> Thank you for your reply Dave. I am going to try that out today with
> some test cases, and will let you know how it goes.
>
> Sam.
>
> Dave Eberly wrote:
> > "sam" <asimnaseer@gmail.com> wrote in message
> > news:1157608649.309785.238610@i42g2000cwa.googlegroups.com...
> > >
> > > Its a cubic curve. Sorry, should have mentioned.
> >
> > Thanks for the clarification.
> >
> > Sometimes people use the term "curvature" to mean
> > magnitude of the second-derivative vector. They are
> > in fact different, so I am assuming you really mean
> > "curvature" in the differential geometric sense.
> >
> > =====
> > About the subdivision of the cubic curve using a test
> > for flatness. The typical fast subdivision involves
> > a stopping criterion involving the magnitude of the
> > second-derivative vector and the length of the current
> > subinterval.
> > void Subdivide (
> > PointList L,
> > float t0, float t1,
> > Point x0, Point x1,
> > Point sd0, Point sd1)
> > {
> > sdmid = 0.5*(sd0 + sd1);
> > dt = t1 - t0;
> > nonlinearity = dt*dt*sdmid;
> > if (SquaredLength(nonlinearity) > epsilon)
> > {
> > tmid = 0.5*(t0 + t1);
> > xmid = 0.5*(x0 + x1 - nonlinearity);
> > insert xmid in L between x0 and x1;
> > Subdivide(t0,tmid,x0,xmid,sd0,sdmid);
> > Subdivide(tmid,t1,xmid,x1,smid,sd1);
> > }
> > }
> > where t is the curve parameter, x is position, and
> > sd is second derivative of position. The epsilon
> > is a user-specified tolerance The subdivision is
> > initiated with
> > x0 = P(0); x1 = P(1); sd0 = P"(0); sd1 = P"(1);
> > L = {x0,x1};
> > Subdivide(L,0,1,x0,x1,sd0,sd1);
> > The following pathological example (in 2D) shows how
> > the subdivision approach can grossly underestimate
> > the maximum curvature. Naturally, your data might
> > not lead to such problems.
> >
> > The control points are (0,0), (0,2), (e,2), and (e,1)
> > for a small positive e. The Bezier curve is
> > p"(t) = (x(t),y(t))
> > with
> > x(t) = e*(3*(1-t)*t^2 + t^3) = e*u(t)
> > y(t) = 6*(1-t)^2*t + 6*(1-t)*t^2 + t^3
> > for 0 <= t <= 1. The last equality in the x(t) equation
> > defines the polynomial u(t). The curvature is
> > k(t) = [x'(t)*y"(t) - x"(t)*y'(t)] / s(t)^3
> > where
> > s(t) = sqrt(x'(t)^2 + y'(t)^2)
> > Notice that y'(r) = 0 for r = sqrt(2)/(sqrt(2)+1). It
> > is easily shown that y"(r) and u'(r) are both not zero.
> > The curvature at this parameter is
> > k(r) = |y"(r)|/(e*u'(r))^2
> > As you make e smaller, k(r) becomes very large. Now at
> > any other value of t, y'(t) is not zero. As you make
> > e smaller, the curvature k(t) becomes smaller. When
> > you look at the graph of the curve, it is essentially
> > two very flat subcurves connected by a tiny curve with
> > very large curvature.
> >
> > The second derivative of the curve is
> > p"(t) = 6*(e*(1-2t),t-2)
> > For small e, the x-component is nearly zero, so the
> > stopping criterion for the subdivision will effectively
> > be controlled by the length of the subintervals. When
> > you evaluate the curvature at the subdivision points,
> > you can always choose e small enough so that those
> > curvatures are nearly zero, in which case you would
> > conclude that the curve is nearly flat everywhere. The
> > reality is that the curve has a tiny region of very
> > large curvature.
> >
> > =====
> > Knowing now that such a subdivision scheme can incorrectly
> > report the maximum curvature, I will suggest another
> > subdivision scheme
> >
> > For a 3D curve p(t), the curvature is
> > k(t) = |Cross(p'(t),p"(t))|/|p'(t)|^3
> > The squared curvature is a ratio of two polynomials in t.
> > An ****ytical approach to computing the maximum curvature
> > is to compute the critical points for k(t)^2. From calculus,
> > These are the t-values where the derivative is zero, and
> > at the endpoints t = 0 and t = 1. The numerator of the
> > derivative of k(t)^2 is
> > n(t) = p11*(p11*p23 - p12*p13) - 3*p12*(p11*p22-p12^2)
> > where p1(t) = p'(t), p2(t) = p"(t), and pij = Dot(pi,pj).
> > If you just count the degrees of the various polynomial
> > terms in n(t), you get degree 9. It turns out that the
> > degree is actually 7. So to compute the critical points,
> > you must compute the roots of a degree 7 polynomial. This
> > sounds heavy handed, but you could take advantage of root
> > bounding schemes to get _approximations_ to the roots and
> > use these to estimate the maximum curvature. You could
> > also compute the Sturm sequence of polynomials for n(t)
> > and determine if n(t) has any roots on [0,1]. If it does
> > not, then the maximum curvature must occur at t = 0 or
> > t = 1.
> >
> > Consider the same ****ysis for a quadratic Bezier curve.
> > The numerator n(t) of the derivative of k(t)^2 happens
> > to be of degree 1. This is easy enough to solve for t,
> > so you have only three critical points to consider, the
> > root for n(t) = 0 and t-values 0 and 1. The short of
> > it is that it is very easy to compute the maximum
> > curvature for a quadratic Bezier curve.
> >
> > Rather than subdivide your cubic Bezier curve into a
> > sequence of line segments (linear Bezier curves as it
> > were), subdivide it into a sequence of quadratic Bezier
> > curves. Each quadratic curve will hopefully be a good
> > estimate of the cubic subcurve it is associated with.
> >
> > A simple way to subdivide is to half the parameter
> > intervals (as in the linear case). Use a least-squares
> > method to fit the cubic curve
> > p(t) = (1-t)^3*p0 + 3*(1-t)^2*t*p1 +
> > 3*(1-t)*t^2*p2 + t^3*p3
> > with a quadratic curve
> > q(t) = (1-t)^2*q0 +2*(1-t)*t*q1 + t^2*q2
> > with q0 = p0 and q2 = p3. Use an integral norm to
> > minimize the integral of squared distances between
> > the two curves. That is,
> > E(q1) = integral[0,1] |q(t) - p(t)|^2 dt
> > When you work through all the math, the unknown control
> > point q1 must be
> > q1 = (-p0 + 3*p1 - 3*p2 + p3)/4
> > Whereas your stopping criterion for the line segments
> > relied on magnitude of second derivatives, the stopping
> > criterion for the quadratic curves can rely on the
> > value of E(q1).
> >
> > --
> > Dave Eberly
> > http://www.geometrictools.com

Actually your expression for calculating a best fit quadratic seesm to
be incorrect. How did you calculate this ? Graphically the quadratics
seem to be way off.

I am going to try to convert to quadratics using the following.

http://www.timotheegroleau.com/Flash...r_in_flash.htm

Sam.

"sam" <asimnaseer@gmail.com> wrote in message
news:1157851165.104066.113390@i42g2000cwa.googlegroups.com...

I wrote:
>> When you work through all the math, the unknown control
>> point q1 must be
>> q1 = (-p0 + 3*p1 - 3*p2 + p3)/4

> Actually your expression for calculating a best fit quadratic seesm to
> be incorrect. How did you calculate this ? Graphically the quadratics
> seem to be way off.

Typographical sign errors. This should be
q1 = (-p0 + 3*p1 + 3*p2 - p3)/4

The Mathematica code just to verify, where 'x' is q1:
p = (1-t)^3*p0 + 3*(1-t)^2*t*p1 + 3*(1-t)*t^2*p2 + t^3*p3
q = (1-t)^2*p0 + 2*(1-t)*t*x + t^2*p3
err = Together[Expand[Integrate[(p-q)^2,{t,0,1}]]]
derr = Simplify[D[err,x]]
Solve[derr == 0, x]

--
Dave Eberly
http://www.geometrictools.com

That did it. Thank you.

Sam.

Dave Eberly wrote:
> "sam" <asimnaseer@gmail.com> wrote in message
> news:1157851165.104066.113390@i42g2000cwa.googlegroups.com...
>
> I wrote:
> >> When you work through all the math, the unknown control
> >> point q1 must be
> >> q1 = (-p0 + 3*p1 - 3*p2 + p3)/4
>
> > Actually your expression for calculating a best fit quadratic seesm to
> > be incorrect. How did you calculate this ? Graphically the quadratics
> > seem to be way off.
>
> Typographical sign errors. This should be
> q1 = (-p0 + 3*p1 + 3*p2 - p3)/4
>
> The Mathematica code just to verify, where 'x' is q1:
> p = (1-t)^3*p0 + 3*(1-t)^2*t*p1 + 3*(1-t)*t^2*p2 + t^3*p3
> q = (1-t)^2*p0 + 2*(1-t)*t*x + t^2*p3
> err = Together[Expand[Integrate[(p-q)^2,{t,0,1}]]]
> derr = Simplify[D[err,x]]
> Solve[derr == 0, x]
>
> --
> Dave Eberly
> http://www.geometrictools.com