[Math] arc length parameterization of planar curve in Matlab

differential-geometryMATLABplane-curves

Let $\gamma (t)$ be a planar curve parameterized by time $t$. For fun, let it be a limacon. In Matlab $\gamma (t)$ looks like this.

% t parameterizes the limacon x=(x(t),y(t)).
h=.1; t0 = 0; T = 2*pi+h;
t = t0:h:T;
a = 1; b = .5;
% This is the cartesian version of the polar curve r = b + a*cos(theta)
x = a/2 + b*cos(t) + (a/2)*cos(2*t);
y = b*sin(t) + (a/2)*sin(2*t);

Assume $\gamma '(t)$ is the time-derivative of $\gamma $ and the arc length $s$ satisfies $s=\int _{t0} ^T \| \gamma '(t) \| \mathrm{d}t $.
The vector s below gives the distance traveled along $\gamma $ as a function of time.

x_t = diff(x);
y_t = diff(y);
s = cumtrapz( sqrt(x_t.^2 + y_t.^2 ) );

I want to reparametrize $\gamma $ in terms of arc length $s$. Basically, I want a vector which gives the time $t$ as a function of $s$. I also need this to be a numerical solution, since in general my curves are quite arbitrary (experimental data), not neat little limacons like I am playing with here. So it will probably be necessary to interpolate between values of the vector s.

Best Answer

I've worked on this a little more, so I guess I'll take a stab at my own question. Here is our limacon again.

the limacon $\gamma(t)$

Notice what we begin with: a time parameter $t$ which runs in equal steps over the interval $[t_0,T]$, the cumulative distance $s$ traveled along the arc length at time $t$, and the Cartesian coordinates of the limacon $\gamma $ at time $t$. We set this much up in the question. What we want instead is an arc length parameter $s$ which runs in equal steps over the interval $[0,L]$, where $L$ is total arc length, together with the time $t(s)$ at $s$ and the Cartesian coordinates $\gamma (t(s))$.

Real quick, let's redo the cumulative distance $s$ traveled along the arc length with derivative (File Exchange link) so that we don't lose samples. 

x_t = derivative(x);
y_t = derivative(y);
s = cumtrapz( sqrt(x_t.^2 + y_t.^2 ) );

We will be using the Matlab function interp1 in order to interpolate the timepoints $t(s_i)$ and coordinates $\gamma (t(s_i))$ which correspond to equally spaced steps $s_i$ along the arc length. Accordingly, we declare the arguments to interp1 with the variable names given in its documentation. The vector X is the original $s_i$, i.e. the distance traveled along the arc length at time $t_i$. The matrix V has column vectors whose rows are $t_i$ (our original time parameter in equal step size), and the coordinates $\gamma (t_i)$.

X = s.';
V = [ t.', x.', y.' ];

We need to interpolate to find the time points $t(s_i)$ corresponding to equally spaced arc length steps along the curve. These equally spaced steps run from $s_0=0$ to $s_N=L$, where $L$ is the total arc length.

L = s(end);
s0 = s(1);
N = length(s);
Xq = linspace(s0,L,N); % equally spaced indices

We can verify that s(N) == sum(diff(Xq)), as should be the case. We use the multidimensional piece-wise linear interpolation capability of Matlab's interpolation function interp1 to find interpolated values for time $t(s_i)$ and coordinates $\gamma (t(s_i))$ which correspond to the equally spaced steps of $s_i$ given in Xq.

Vq = interp1(X,V,Xq);

A call to diff(Xq) shows that we are moving in steps of constant arc length along $\gamma $, while a call to diff(Vq) shows that the time steps are no longer equal (first column). Thus, the result of our interpolation is:

  • arc length $s_i$ in constant steps (Xq)
  • corresponding timepoints $t_i$, no longer in constant steps (first column of Vq)
  • coordinates $\gamma(t(s))$ as function of arc length (second and third columns of Vq)

So the arc length parametrization of the limacon $\gamma (t(s))$ is expressed in Cartesian coordinates as xs and ys. 

xs = Vq(:,2);
ys = Vq(:,3);

These vectors give the position of a particle moving along $\gamma $ as a function of $s$, where $s$ increments in steps of 0.1058. Here are the two parameterizations of the curve $\gamma $, namely $\gamma (t)$ in blue and $\gamma (t(s))$ in red.

two parametrizations

I'll have to do some more work if I want to allow periodic solutions, since the current Matlab script does not in general. To see this, just look at the distance between the last and first sample of $\gamma (t(s))$ in red above. If we were to advance one more step of 0.1058 units along the limacon we would not land on a point we already visited, even though we know limacons are $2\pi$-periodic. I think this basically amounts to choosing a different step size for $s$ than 0.1058.

Related Question