Arc Length of Bézier Curves

See also: answers with code on GameDev.SE

How can I find out the arc length of a Bézier curve? For instance, the arc length of a linear Bézier curve is simply:

$$s = \sqrt{(x_1 - x_0)^2 + (y_1 - y_0)^2}$$

But what of quadratic, cubic, or nth-degree Bézier curves?

$$\mathbf{B}(t) = \sum_{i=0}^n {n\choose i}(1-t)^{n-i}t^i\mathbf{P}_i$$

Solution 1:

I showed in the comments how to get Mathematica to generate the arclength function for a quadratic Bézier curve; for this answer I'll give an explicit derivation.

Consider the parametrization

$$\begin{align*}x&=(1-u)^2 x_1+2u(1-u) x_2+u^2 x_3 \\ y&=(1-u)^2 y_1+2u(1-u) y_2+u^2 y_3\end{align*}$$

Letting $\Delta x_i=x_{i+1}-x_i$ and similarly for $\Delta y_i$, the arclength integral corresponding to this parametrization is

$\displaystyle \scriptsize 2\int\sqrt{\Delta x_1^2+\Delta y_1^2+2\left(\Delta x_1\left(\Delta x_2-\Delta x_1\right)+\Delta y_1\left(\Delta y_2-\Delta y_1\right)\right)u+\left(\left(\Delta x_2-\Delta x_1\right)^2+\left(\Delta y_2-\Delta y_1\right)^2\right)u^2}\mathrm du$

We let $c=\Delta x_1^2+\Delta y_1^2$, $b=\Delta x_1\left(\Delta x_2-\Delta x_1\right)+\Delta y_1\left(\Delta y_2-\Delta y_1\right)$, and $a=\left(\Delta x_2-\Delta x_1\right)^2+\left(\Delta y_2-\Delta y_1\right)^2$ to further simplify things.

Consider now the integral

$$2\int \sqrt{c+2bu+au^2}\mathrm du$$

Completing the square yields

$$2\sqrt{a}\int \sqrt{\left(u+\frac{b}{a}\right)^2+\frac{ac-b^2}{a^2}}\mathrm du$$

Skipping the details (but see here for how one might derive the answer), the integral evaluates to

$$\sqrt{a}\left(\left(u+\frac{b}{a}\right)\sqrt{\left(u+\frac{b}{a}\right)^2+\frac{ac-b^2}{a^2}}+\left(\frac{ac-b^2}{a^2}\right)\mathrm{arsinh}\left(\frac{au+b}{\sqrt{ac-b^2}}\right)\right)$$

or

$$\left(u+\frac{b}{a}\right)\sqrt{c+2bu+au^2}+\left(\frac{ac-b^2}{a^{3/2}}\right)\mathrm{arsinh}\left(\frac{au+b}{\sqrt{ac-b^2}}\right)$$

As I've mentioned in the comments, the closed form for the quadratic case is quite complicated (even more so for the cubic case), and you're better off with using numerical quadrature to compute the arclength.

Solution 2:

See these papers:

• Approximate Arc Length Parametrization, in SIBGRAPI 1996.

• Adaptive sampling of parametric curves, in Graphics Gems V, 1995.

• Computing the arc length of parametric curves, IEEE Computer Graphics and Applications, 1990.

• Point-based methods for estimating the length of a parametric curve, Journal of Computational and Applied Mathematics, 2006.

Solution 3:

See also "Adaptive subdivision and the length and energy of Bézier curves" by Jens Gravesen in Computational Geometry Volume 8, Issue 1, June 1997, Pages 13–31.

The idea is that the arc length of Bezier curve lies between chord-length (distance from first to last control point) and polygon-length (distance between each successive pair of control points). Repeated subdivision of the curve shrinks the arc length interval, up to arbitrary close precision.

Solution 4:

I suspect the answer you seek differs from the question you asked. Let me consider a few 3-D scenarios (0ne post each).

SCENARIO #1: You load up a 3D program like Blender, and add a Bezier curve to the scene. There are control handles, and the curve is approximated by segments. Increase the resolution of the curve, and the approximation gets better and better. Zoom far away from the curve, and at some point the resolution of your curve will exceed the resolution displayed. It would be nice to cull unnecessary points.

If this is what you mean, you are in luck:
The drawn curve is already a linear approximation. Increase and decreasing the resolution has no effect on the order of the curve. Let us build such a curve.

Let the order be 3, and the given points be A,B,C,D. and let the 3D curve be given by
$(1a)\;\;\;\;\; {\bf x}(t) = A + 3 P t + 3Q t^2 + R t^3, \;\;\;\;\; 0\le t \le 1$
${\small \hspace{60pt} P = (B\!-\!A),\;\; Q = (C\!-\!2B\!+\!A),\;\; R= (D\!-\!3C\!+\!3B\!-\!A)}$
or, immediately in the given points
$(1b)\;\;\;\;\; {\bf x}(t) = A(1\!-\!t)^3 + 3 Bt (1\!-\!t)^2 + 3C t^2 (1\!-\!t) + D t^3 \;\;\;\;\; 0\le t \le 1$

The points A and D are the endpoints of the curve, and if you join segments AB, CD, you will find that these are the familiar handles. Now consider two tasks

1. Change the resolution. How do we draw such a curve? Say we want n segments. Then we will need n+1 t-values:
$\;\;\;\;\;\{t_k\} = 0, \tfrac{1}{n}, \tfrac{2}{n}, ..., \tfrac{n}{n} = 1,\;\;$ or
$\;\;\;\;\;\{t_k\} = {k/n}, \;\;\;k = 0, \ldots, n $
Then evaluate (1a) or (1b) at the n t-values. Join them with straight lines. This curve is your approximation. The distance between points is the straight line distance.

Now let's think this through. Even if this is not what you mean, are you sure you want the arc length? I have two points in space, A, and B. I will know if two points overlap on my screen by checking their ordinary Euclidean distance. The distance between them along the curve is not only much harder to calculate, it doesn't answer the question.

2. Increase the number of control points/handles. High-order polynomials are unruly monsters. Rather than increase the order of the curve, I will concatenate 3rd order curves (let them share common endpoints, and maybe tangents). This has the familiar look: a point on the curve has a handle that extends on both sides. As a general rule, we will use 3rd order curves everywhere, so that between any two control points on the spline, there is a single, 3rd order curve, and the handles give the direction, and strength, of the tangent at that point,
like this.

In both cases, I avoid raising the order, All curves are order 3.

Let me take this implementation seriously. I am not convinced that I need to change the resolution, except by some very general rules (say, based on camera distance and bounding box). Here is why.

Consider representation (1b). You posted in a Game forum, so I'll assume you are willing to use a graphics card or SSE2 to do some vector math. I choose an arbitrary resolution: say, 32 points (31 segments). Let me agree not to change this number, but instead use it everywhere, for all splines displayed on the screen.

Consider the kth point of my spline. I will write it as a vector product:

$\left[\begin{array}{} A & B & C&D \end{array}\right] \left[\begin{array}{} (1\!-\!t_k)^3 \\(1\!-\!t_k)^2 t_k \\ (1\!-\!t_k){t_k}^2\\{t_k}^3 \end{array}\right]\;\;\;$ or, $\left[\begin{array}{} A & B & C&D \end{array}\right] \cdot T_k$,
where Tk is a column vector. Then the 32 poitns on my spline are $\left[\begin{array}{} A & B & C&D \end{array}\right] \left[ \left[\begin{array}{} \;\\ \!T_0 \!\\\; \end{array}\right] \left[\begin{array}{} \;\\ \!T_1\! \\\; \end{array}\right] \left[\begin{array}{} \;\\ \!\ldots \! \\\; \end{array}\right] \left[\begin{array}{} \;\\ \!T_{31}\! \\\; \end{array}\right] \right] $

The 4x32 matrix on the right, containing all the powers of tk, is constant, and will be used for all splines on my screen. The whole problem can be passed to the graphics card (or SSE) as the (partitioned) matrix multiplication of [Pts] $\cdot$ [Constants], and then drawing the straight lines joining points on the same curve.

Being the thing to be done
Onward. (TO DO: EDIT.)

Solution 5:

I worked out the closed form expression of length for a 3 point (quadratic) Bezier (below). I've not attempted to work out a closed form for 4+ points. This would most likely be difficult or complicated to represent and handle. However, a numerical approximation technique such as a Runge-Kutta integration algorithm would work quite well by integrating using the arc length formula.

Here is some Java code for the arc length of a 3 point Bezier, with points a,b, and c.

    v.x = 2*(b.x - a.x);
    v.y = 2*(b.y - a.y);
    w.x = c.x - 2*b.x + a.x;
    w.y = c.y - 2*b.y + a.y;

    uu = 4*(w.x*w.x + w.y*w.y);

    if(uu < 0.00001)
    {
        return (float) Math.sqrt((c.x - a.x)*(c.x - a.x) + (c.y - a.y)*(c.y - a.y));
    }

    vv = 4*(v.x*w.x + v.y*w.y);
    ww = v.x*v.x + v.y*v.y;

    t1 = (float) (2*Math.sqrt(uu*(uu + vv + ww)));
    t2 = 2*uu+vv;
    t3 = vv*vv - 4*uu*ww;
    t4 = (float) (2*Math.sqrt(uu*ww));

    return (float) ((t1*t2 - t3*Math.log(t2+t1) -(vv*t4 - t3*Math.log(vv+t4))) / (8*Math.pow(uu, 1.5)));