Re: Fast maximum curvature calculation for a 3d bezier curve
- From: "Dave Eberly" <dNOSPAMeberly@xxxxxxxxxxxxxxx>
- Date: Thu, 07 Sep 2006 13:58:58 GMT
"sam" <asimnaseer@xxxxxxxxx> wrote in message
news:1157608649.309785.238610@xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
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 analytical 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 analysis 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
.
- Follow-Ups:
- References:
- Prev by Date: Evaluating an object recognition algorithm
- Next by Date: Re: Stochastic Positioning of A Point in A 3D Gaussian Distribution
- Previous by thread: Re: Fast maximum curvature calculation for a 3d bezier curve
- Next by thread: Re: Fast maximum curvature calculation for a 3d bezier curve
- Index(es):
Relevant Pages
|
