Arc length parametrization: catenary [part 1]

mathematics / gamedev

Wed, Dec 29, 21

General case

Arbitrary curve in ${\mathbb{R}}^n$

Let a parametric curve $C$ be described by a continuously differentiable and injective function $\gamma :I\to {\mathbb{R}}^n$ for $I=[a,b]\subset \mathbb{R}$. Let $t\in I$ be the curve parameter. Then, we can define the arc length of the curve up to the value $t$ of the parameter as:

The full length of the curve is obtained by integrating over the whole interval $I$:

If we want to find the value of $t$ such that the arc length $s$ has a given value, then all we need to do in theory, is to compute the inverse function $t(s)$. Then the function $\gamma \circ t$ is parametrized in terms of the arc length. If we need the value that $\gamma$ takes exactly at the midpoint of the curve $C$, for example, then we calculate $(\gamma \circ t)(\frac{L}{2})$. We can additionaly compose this function with $\iota :u\mapsto Lu$, so that we just need to feed it the length fraction as a number in $[0,1]$:

And the value of $\gamma$ at the midpoint of $C$ is simply $\bar{\gamma }(\frac{1}{2})$.

Planar curves defined by function graphs

To set these ideas on a special case, let’s consider the curve obtained by graphing some arbitrary $C^1$-function. Let a curve $C$ be specified as the graph of a function $\gamma :I\to \mathbb{R}$. Then all the points $\overrightarrow{r}\in {\mathbb{R}}^2$ belonging to this graph are of the form:

Then $\text{1}$ takes the form:

If we’re looking for a closed-form expression, then $\text{1}$ must be expressible in terms of elementary functions, and the expression must be invertible. This is not possible in general as $\text{1}$ is often an elliptic integral, as you might suspect from seing that bad boy radical in the integral $\text{5}$. And if it is not the case, then there is always the possibility that no closed-form expression exists for the inverse of $s$. In fact, there are only a handful of curves for which a closed-form arc length expression exists, and we’ll have a look at one of them in the following section.

On a side note, we have some adjustments to make to $\text{4}$ if we intend to plot $\bar{\gamma }$ using the normalized arc length $\bar{s}\in [0,1]$:

Arc length parametrization of the catenary

A catenary is the shape assumed by a hanging chain or rope, supported only at its both ends, under the effect of its own weight. If you ever saw an electric cable hanging between two posts, or the Gateway Arch in St. Louis, you know what a catenary looks like. It resembles a parabola, but a more rounded version of it. I won’t prove to you why hanging cables assume this shape, but I encourage you to have a look at a proof involving variational calculus if you’re interested in physics, it’s worth it. We will assume that the catenary is given by the following general expression:

That’s right, a good ol’ hyperbolic cosine! Here, $p$, $q$ and $a$ are curve parameters. It is more practical, in terms of game development, to think of the anchor points $P_1$ and $P_2$ and the total length $L$ as the input parameters to this problem, as it makes more sense to control them in an editor, rather than the curve parameters themselves. So $p$, $q$ and $a$ will be expressed later on as functions of $P_1$, $P_2$ and $L$.

Solving for the arc length

Let $\overrightarrow{P_1}=(x_1,y_1)$ and $\overrightarrow{P_2}=(x_2,y_2)$ be the position of the anchor points in the plane. We suppose that $x_1<x_2$ and $y_1\le y_2$. If it’s not the case, you can always swap $P_1$ and $P_2$ and reflect the curve about the midline. We also note $v=y_2-y_1$ and $h=x_2-x_1$. We suppose that the total length $L$ of the catenary verifies the inequality $L>\Vert \overrightarrow{P_2}-\overrightarrow{P_1}\Vert$ as the curve cannot be shorter than the straight line connecting the anchor points.

By the chain rule, we have that:

Thus, according to $\text{5}$, the arc length of the catenary is given by:

And because $s(x_1)=0$ we have a value for the integration constant:

So the final expression for $s$ is:

And substituting $x_2$ gives an expression for the total length $L$:

Because we have the insolent good fortune to work with such a nicely behaved function, we can invert $\text{11}$ to express $x$ as a function of $s$:

with $\mathrm{arcsinh}(u)\equiv \ln (u+\sqrt{1+u^2})$. And with this we can easily obtain an arc length parametrized $\bar{\gamma }$ using $\text{3}$:

Determining parameters

If we want to find the value of $a$ such that the catenary between $P_1$ and $P_2$ is of length $L$, then we need to solve the following equation:

Now, the issue with this one, is that it’s a transcendental equation in $a$, so it can only be solved numerically. A few iterations of Newton-Raphson will do the trick:

Our update scheme $\text{20}$ requires us to input an initial guess $a_0$ for the value of $a$, and the convergence is quite sensitive so we would better not be too far off. In practice, I found that $a_0$ can be selected by iteratively incrementing $a$ by a geometrically growing step size until the sign of $f(a)$ changes. This works because$f$ only has a single zero, as it’s monotonically decreasing. More on that later. Then, $p$ and $q$ are given by:

Implementation

My original implementation is scattered across multiple C++ files in a project called Kibble, located on my personal GitHub account. You can have a look at this folder, more specifically the files named catenary.h/cpp and numeric.h/cpp, if you’re looking for a C++ implementation. There is also an example here on how to use the implementation. Kibble is a utility library I’ve been working on for some time, that features many software components I end up always reusing in my game development projects.

But I find it more practical to plot stuff in Python, so let’s have some good fun with numpy, scipy and matplotlib:

import numpy as np
from scipy import optimize
from matplotlib import pyplot as plt

First, we implement our catenary function from $\text{7}$:

def catenary(x, a, p, q):
    return a * np.cosh((x-p)/a) + q

The arc length parametrization is handled by another short function that makes use of $\text{13}$:

def arclen_remap(s_bar, L, a, p, c):
    return a * np.arcsinh((s_bar * L - c) / a) + p

We move on to determining the values of $a$, $p$, $q$ and $c$ from the input parameters. We will use scipy.optimize.newton() to find $a$ by solving $\text{15}$, so we must implement functions for $\text{18}$ and $\text{19}$:

def f(a, h, v, L):
    return 2 * a * np.sinh(0.5 * h/a) - np.sqrt(np.power(L, 2)-np.power(v, 2))

def fprime(a, h):
    return 2 * np.sinh(0.5 * h/a) - (h/a) * np.cosh(0.5 * h/a)

We also need an initial guess for the Newton-Raphson algorithm to work:

def nr_first_guess(ff, start_x, start_step, alpha):
    xx = start_x
    step = start_step
    yy = ff(xx)
    yy_prev = yy

    while yy * yy_prev > 0:
        yy_prev = yy
        xx += step
        step *= alpha
        yy = ff(xx)

    # Backtrack a bit
    return xx - 0.5 * step / alpha

This function iteratively increments the value of the argument xx fed to an input function ff. The steps themselves grow geometrically with a common factor alpha which is bigger than $1$. We break from the loop when the sign of the function has changed. At this point, we know we’re past the zero of the input function, so we interpolate between the current and previous value of xx to output a more accurate value for the zero. This is a quick and dirty solution, you can get fancier with a binary search if you need to.

Then, using the previous functions, as well as formulas $\text{21}$, $\text{27}$ and $\text{10}$, we are able to compute the values of all the parameters that the catenary() function expects as inputs:

def get_params(p1, p2, L):
    hv = p2 - p1
    m = p1 + p2
    def f_bind(a): return f(a, *hv, L)
    def fprime_bind(a): return fprime(a, hv[0])

    a0 = nr_first_guess(f_bind, 0.1, 0.01, 1.8)
    a = optimize.newton(f_bind, a0, fprime_bind)

    p = 0.5 * (m[0] - a * np.log((L+hv[1])/(L-hv[1])))
    q = 0.5 * (m[1] - L / np.tanh(0.5 * hv[0]/a))
    c = -a * np.sinh((p1[0]-p)/a)

    return a, p, q, c

The values fed to nr_first_guess() are somewhat arbitrary, but we need to keep in mind that we should not start too close to $0$ as $f(a)$ diverges there. And basically, we’re done! We can compute catenary parameters like so:

p1 = np.array([0.5, 0.6])
p2 = np.array([4.1, 2.5])
min_len = np.linalg.norm(p2-p1)
L = 2 * min_len
a, p, q, c = get_params(p1, p2, L)

Note that we calculate the minimal length as the Euclidean distance between $P_1$ and $P_2$, so we can make sure the input distance L will be greater than that. To get points on the catenary using the initial parametrization we would do:

xp = np.linspace(p1[0], p2[0], 100)
yp = catenary(xp, a, p, q)

And to use the arc length parametrization we would do:

s_bar = np.linspace(0, 1, 100)
xp_arclen = arclen_remap(s_bar, L, a, p, c)
yp_arclen = catenary(xp_arclen, a, p, q)

So let’s plot a few catenary curves with varying lengths, and compare the two parametrizations:

You can see how the “default” parametrization leads to uneven segment lengths, the effect is more pronounced for longer curves. The segments obtained using the arc length parametrization on the other hand, are of equal lengths. If you were to skin a curve like this in a game, the arc length parametrization would completly get rid of texture stretching artifacts. Do note however that the segments tend to deviate more from the curve, so a higher number of sampling points must be chosen in order to retain a rounded shape. Have a look at the full source code (link below) if you want more detail on how to plot this.

In the next article, we are going to implement an arc length parametrization for cubic Hermite splines. The problem will be solved entirely numerically as no closed-form expression exists for these curves. Hope you enjoyed, see you on the next one!