of Cubic Bezier Curves

Some time ago I asked:

and promised to summarize the responses. Herewith the summary.Given a cubic bezier defined by four points P0 to P3, and a parameter value t, what is the expression which gives the length of the curve from its origin at P0 to the parameter value t?

The length of a segment of a parametric curve P(t) = [x(t) y(t) z(t)] is:

/t1 s = \ |P'(t)| dt / where P'(t) is the curve tangent vector at t t0 dx(t) dy(t) dz(t) P'(t) = [x' y' z'] = [----- ----- -----] dt dt dtUsing equivalent alternative notations:

/t1 s = \ sqrt(P'(t) .dot. P'(t)) dt / t0 /t1 = \ sqrt(x'*x' + y'*y' + z'*z') dt / t0There is no closed-form solution, in general, to this integral for cubic polynomial curves. It is common to use Gauss-Legendre quadrature to efficiently evaluate the integral.

Once you have a function to evaluate curve length, you can constrain the length of a parametric curve defined by control points using an iterative numerical procedure. If you are using Bezier or Hermite curves, it is straightforward, for example, to fix the start and end control points and tangent vector directions, while adjusting the tangent vector magnitudes to converge to the desired curve length constraint (ie., for cubic Bezier curves, slide the two intermediate control points along the lines defined by the tangents at the start and end control points). The minimum curve length, of course, is the straight line distance between the start and end control points. The curve length is a monotonically increasing function of start and end control point tangent vector magnitudes, so convergence should not be a problem for any reasonable iterative numerical method.

Here's a simple approximation algorithm based on Simpson:

#define sqr(x) (x * x) #define _ABS(x) (x < 0 ? -x : x) const double TOLERANCE = 0.0000001; // Application specific tolerance extern double sqrt(double); struct point2d { double x, y; }; double q1, q2, q3, q4, q5; // These belong to balf() //--------------------------------------------------------------------------- double balf(double t) // Bezier Arc Length Function { double result = q5 + t*(q4 + t*(q3 + t*(q2 + t*q1))); result = sqrt(result); return result; } //--------------------------------------------------------------------------- // NOTES: TOLERANCE is a maximum error ratio // if n_limit isn't a power of 2 it will be act like the next higher // power of two. double Simpson ( double (*f)(double), double a, double b, int n_limit, double TOLERANCE) { int n = 1; double multiplier = (b - a)/6.0; double endsum = f(a) + f(b); double interval = (b - a)/2.0; double asum = 0; double bsum = f(a + interval); double est1 = multiplier * (endsum + 2 * asum + 4 * bsum); double est0 = 2 * est1; while(n < n_limit && (_ABS(est1) > 0 && _ABS((est1 - est0) / est1) > TOLERANCE)) { n *= 2; multiplier /= 2; interval /= 2; asum += bsum; bsum = 0; est0 = est1; double interval_div_2n = interval / (2.0 * n); for (int i = 1; i < 2 * n; i += 2) { double t = a + i * interval_div_2n; bsum += f(t); } est1 = multiplier*(endsum + 2*asum + 4*bsum); } return est1; } // //--------------------------------------------------------------------------- // double BezierArcLength(point2d p1, point2d p2, point2d p3, point2d p4) { point2d k1, k2, k3, k4; k1 = -p1 + 3*(p2 - p3) + p4; k2 = 3*(p1 + p3) - 6*p2; k3 = 3*(p2 - p1); k4 = p1; q1 = 9.0*(sqr(k1.x) + sqr(k1.y)); q2 = 12.0*(k1.x*k2.x + k1.y*k2.y); q3 = 3.0*(k1.x*k3.x + k1.y*k3.y) + 4.0*(sqr(k2.x) + sqr(k2.y)); q4 = 4.0*(k2.x*k3.x + k2.y*k3.y); q5 = sqr(k3.x) + sqr(k3.y); double result = Simpson(balf, 0, 1, 1024, 0.001); return result; }

From: Jens Gravesen

By subdividing the curve at parameter value t you only have to find the length of a full Bezier curve.

If you denote the length of the control polygon by L1 i.e.:

L1 = |P0 P1| +|P1 P2| +|P2 P3|and the length of the cord by L0 i.e.:

L0 = |P0 P3|then

L = 1/2*L0 + 1/2*L1is a good approximation to the length of the curve, and the difference

ERR = L1-L0is a measure of the error. If the error is to large, then you just subdivide curve at parameter value 1/2, and find the length of each half.

If m is the number of subdivisions then the error goes to zero as 2^-4m.

If you dont have a cubic curve but a curve of degree n then you put

L = (2*L0 + (n-1)*L1)/(n+1)A reference is:

Jens Gravesen: "Adaptive subdivision and the length of Bezier curves" mat-report no. 1992-10, Mathematical Institute, The Technical University of Denmark.

/************************ split a cubic bezier in two ***********************/ static void bezsplit(V, Left, Right) point *V; /* Control pts */ point *Left; /* RETURN left half ctl pts */ point *Right; /* RETURN right half ctl pts */ { int i, j; /* Index variables */ point Vtemp[4][4]; /* Triangle Matrix */ /* Copy control points */ for (j =0; j <= 3; j++) Vtemp[0][j] = V[j]; /* Triangle computation */ for (i = 1; i <= 3; i++) { for (j =0 ; j <= 3 - i; j++) { Vtemp[i][j].x = 0.5 * Vtemp[i-1][j].x + 0.5 * Vtemp[i-1][j+1].x; Vtemp[i][j].y = 0.5 * Vtemp[i-1][j].y + 0.5 * Vtemp[i-1][j+1].y; } /* end for i */ } /* end for j */ for (j = 0; j <= 3; j++) Left[j] = Vtemp[j][0]; for (j = 0; j <= 3; j++) Right[j] = Vtemp[3-j][j]; } /* end splitbez */ /********************** add polyline length if close enuf *******************/ static void addifclose(V,length, error) point *V; double *length; double error; { point left[4], right[4]; /* bez poly splits */ double len = 0.0; /* arc length */ double chord; /* chord length */ int index; /* misc counter */ for (index = 0; index <= 2; index++) len = len + V2DistanceBetween2Points(&V[index],&V[index+1]); chord = V2DistanceBetween2Points(&V[0],&V[3]); if((len-chord) > error) { bezsplit(V,left,right); /* split in two */ addifclose(left, length, error); /* try left side */ addifclose(right, length, error); /* try right side */ return; } *length = *length + len; return; } /* end addifclose */ /************************* arc length of a 2d bezier ************************/ double arclen(V, error) point *V; double error; { double length; /* length of curve */ addifclose(V, &length, error); /* kick off recursion */ return(length); /* that's it! */ } /* end arclen */

The same code can be used for a quick Bezier plot routine by replacing "addifclose" with a "drawifclose" that draws the polygon instead of adding the length to the arc length. The value of "error" can be based on your screen parameters. My thanks to Jens for replying, and I hope his method at least gets into some FAQ because it deserves to be better known. Again, my code is a quick hack and *seems* to work; if it sends your process into hyperspace please email me the corrections so I can fix my version :-)

There were a couple of other responses. Several suggested plotting the curve as a polyline and adding up the lengths, which is what I was doing and which is *slow* compared to the Gravesen method. A couple of people also recommended the book "Curves and Surfaces for Computer Aided Geometric Design" by Gerald Farin. I have one on order and will post a review when I get it. The recommendations have been such that it, too, seems a candidate for a FAQ.

Brian Guenter noted that he had

... published a fast numerical algorithm for computing the arc length of space curves. It appeared in IEEE Computer Graphics and Applications in 1990. I really should be able to give you a better reference than this but I've just moved to Seattle and I don't have all my papers with me. The algorithm is simple and fast and does just what you want it to.

I was unable to track this down, but it should be mentioned.