arc length

Posted on by svark

In this post we look at an algorithm to compute arc lengths of a differentiable curve. First lets recall the definitions:

Let \mathbf{r} be a the curve parametrized by t \in [0,1].
\displaystyle \gamma(t) = \int_{[0,t]}\sqrt{\mathbf{r}'(t) \cdot \mathbf{r}'(t) } \,dt

The thing that is nice about arc lengths is that they are integrals and behave well even when our curve fails to be differentiable at places. The thing that is not nice (computationally) is that for general polynomial curves this integral may not have closed form expression. Try this in wolfram alpha for example:

Integrate[Sqrt[1 + x + 2 x^2 + x^5], x]

We thus have to resort to numerical techniques, such as the Gaussian quadrature to get an answer. Gaussian quadrature is a remarkably effective tool too(esp for smooth functions) and whats more the convergence is fast for all practical purposes( error ~ O(f^{(2n)}(\eta)) for n term evaluation where f is what is being integrated). This does not mean we should stop playing with this and see where it can take us. The complexity comes from the square root right, lets get rid of it and start playing.

Lets write
\displaystyle \Gamma(t) = \int_{[0,t]} {\mathbf{r}'(t) \cdot \mathbf{r}'(t) } \,dt

This is a tractable problem and definitely can be evaluated for polynomial curves but it is way off the result. Lets us see the connection between the two

\displaystyle \Gamma'(t) = (\gamma'(t))^2

This depends purely on the derivatives and there appears no direct way to retrieve the arc length from this. Lets try to bringing in a friend of of \gamma to see a better connection. Define \displaystyle \eta(t) = \gamma(t)/\gamma'(t) and differentiate.

\displaystyle \eta'(t) = 1 - \frac{\gamma(t)\gamma''(t)}{\gamma'(t)^2}
or

\displaystyle \eta'(t) = 1 - \eta(t)\frac{\Gamma''(t)}{2 \Gamma'(t)}

Since we have come this far, let take this as far as it gets us. Rewrite the above with \eta on lhs

\displaystyle \eta(t)\Gamma''(t) = 2 { \Gamma'(t)} (1 - \eta'(t))
or
\displaystyle \eta(t)\Gamma''(t) + \eta'(t)\Gamma'(t) = 2\Gamma'(t) - \Gamma'(t) \eta'(t)
or
\displaystyle (\eta(t) \Gamma'(t))' = (2 - \eta'(t))\Gamma'(t)

Note that lhs is ready for integration but rhs is not. If we can express \eta'(t) (in the rhs) as a power series in terms of \Gamma(t) then integrating rhs is also easily done. We can worry about the fact that power series exists later(see, we are in real world) lets follow this trail anyway.

Lets write
\displaystyle \eta'(t) = c_0 + c_1 \Gamma(t) + c_2 {\Gamma(t)}^2 + \ldots

then
\displaystyle \eta(t) = c_0 t + c_1 \int \Gamma(t) dt + c_2 \int \Gamma(t)^2 dt + \ldots

Our job is to find c_0 , c_1,\ldots such that both sides match. At this point we have to break off and do a hand calculation which reveals the following (or use as CAS tool to get these identities as you prefer)

\displaystyle c_0 = 1
\displaystyle c_1= -\frac{\Gamma ''(0)}{2 \Gamma '(0)^2}
\displaystyle c_2= \frac{7 \Gamma ''(0)^2-4 \Gamma ^{(3)}(0) \Gamma '(0)}{8 \Gamma '(0)^4}
\displaystyle c_3= \frac{-77 \Gamma ''(0)^3-12 \Gamma ^{(4)}(0) \Gamma '(0)^2+74 \Gamma ^{(3)}(0) \Gamma '(0) \Gamma ''(0)}{48 \Gamma '(0)^6}
\displaystyle c_4= \frac{1155 \Gamma ''(0)^4-32 \Gamma ^{(5)}(0) \Gamma '(0)^3+192 \Gamma ^{(3)}(0)^2 \Gamma '(0)^2+316 \Gamma ^{(4)}(0) \Gamma '(0)^2 \Gamma ''(0)-1526 \Gamma ^{(3)}(0) \Gamma '(0) \Gamma ''(0)^2}{384 \Gamma '(0)^8}

You can see additional terms here. Calculated using

pow

We can retrieve the arc length from this by using the identity

\displaystyle \gamma(t) = \gamma'(t) \eta(t) .
or

\displaystyle \gamma(t) = \gamma'(t) (c_0 t + c_1 \int \Gamma(t) \,dt + c_2 \int \Gamma(t)^2 \,dt + \ldots )

But wait, there is another way we can retrieve \gamma without involving additional integrals:
\displaystyle \gamma(t) = \frac{ \gamma'(t)^2 \cdot (1 - \eta'(t) )}{\gamma''(t)}
or

\displaystyle \gamma(t) = -{\gamma'(t)^2/\gamma''(t)} \cdot (c_1 \Gamma(t) + c_2 \Gamma(t)^2 + \ldots )

Lets unit test this with a basic function:

\displaystyle y = x^2/{\sqrt 2}

parametrized by x = t/{\sqrt 2}, y = \frac{t^2}{2{\sqrt 2}}

The closed form of the arc length of this parabola is:
\displaystyle (t \sqrt{1+t^2}+\sinh^{-1} {t})/(2 \sqrt{2})
(courtesy wolfram alpha) and this has the series expansion

\displaystyle t/\sqrt{2} + t^3/(6 \sqrt{2}) - t^5/(40 \sqrt{2}) + t^7/(112 \sqrt{2}) + \ldots

Using our calculations we get,
\displaystyle c_0 = 1, c_1 = 0, c_2 = - 4, c_3 = 0, c_4 = 32 from which

\gamma(t) = \displaystyle {\sqrt \frac{1 + t^2}{2}} (t - 4 (t^3/12 + t^5/30 + t^7/252)) + 32 ( t^5/80 + t^7/84 + O(t^9) )

which has the power series expansion:

t/\sqrt{2}+t^3/(6 \sqrt{2})-t^5/(40 \sqrt{2})+(3057 t^7)/(5040 \sqrt{2}) + O(t^9) )
if one uses the expansion that does not involve integrals we get

t/\sqrt{2}+t^3/(6 \sqrt{2})-(301 t^5)/(72 \sqrt{2})-(275 t^7)/(48 \sqrt{2})+O(t^9)
which has error an order of magnitude larger than the previous expansion.

Parametrization

Now, having come this far let us also look to see if we can parameterize the curve using arc length.
Given arc length s, our task is to find t such that s = \gamma(t)
Using our series expansion, we may write

\displaystyle t = s/\gamma'(t) - \int_{[0,t]}\sum_i c_i \Gamma_i(t)\,dt

which suggests fixed point algorithm using t_{j+1} = F(t_j) where
\displaystyle F(t) = s/\gamma'(t) - \int_{[0,t]}\sum_i c_i \Gamma_i(t)\,dt
this would converge provided \left|F'(t)\right| < 1 which would be true if
\displaystyle |\frac{s \gamma''(t)}{\gamma'(t)^2} + \sum_i c_i\Gamma_i(t) | < 1
or
\displaystyle \left |s - \gamma(t)\right | < \left |\frac{\gamma'(t)^2}{\gamma''(t)} \right |
The scheme will hence get us to the fixed point provided \gamma''(t)^2 stays away from zero.

Leave a comment