Is there a way to divide a discretized curve (or a set of points) into another where the segments are of equal length?

For example, in Matlab I can do the following to get an ellipse with equispaced points or discretization.

```
N = 50; mu0 = 0.5;
theta = (2/N) * (-N/2 : N/2 ) * pi;
% ellipse
X = cosh(mu0) * cos(theta); Y = sinh(mu0) * sin(theta);
dxdt = gradient(X, X(2) - X(1)); % get gradient of X
dydt = gradient(Y, X(2) - X(1)); % get gradient of Y
ds = sqrt(dxdt.^2 + dydt.^2) ; % get arc length (differential)
s = cumtrapz(ds); % calculate cumulative arc length (ds)
S = transpose(s);
V = [ transpose(X), transpose(Y) ];
N = length(s); % get number of total points
Xq = linspace(s(1),s(end),N);
Vq = interp1(S,V,Xq, 'linear'); % interpolate along equal length
% return equispaced points
xs = Vq(:,1);
ys = Vq(:,2);
```

I have tried using the `Interpolations.jl`

package to no avail. If there is another method to partition a collection of coordinates into segments of equal length, I’d be equally grateful.

The differences between arc length in both discretizations.

1 Like

I am not aware of any package that does this automatically. But we can easily reproduce the Matlab code. I’d recommend `DataInterpolations.jl`

for the interpolation step.

There is NumericalIntegration.jl which defines cumulative integration routines (and ironically depends on Interpolations.jl…), but for the problem at hand it’s easier to write `cumtrapz`

oneself.

```
using DataInterpolations
using LinearAlgebra
gradient(y::Vector, t::AbstractVector) = diff(y)./diff(t)
function cumtrapz(y, t)
cs = similar(y)
cs[begin] = 0
for i in 2:lastindex(y)
cs[i] = cs[i-1] + 0.5*(t[i]-t[i-1])*(y[i]+y[i-1])
end
cs
end
# sampling points
N = 2*50+1
t = range(-pi, pi, length=N)
# function to sample from
f(t) = [sin(t), cos(3t)]
y = f.(t)
dy = gradient(y, t)
norm_dy = norm.(dy)
S = cumtrapz(norm_dy, t)
Titp = LinearInterpolation(t, S)
S_subsamples = range(S[1], S[end], length=20)
titp = Titp.(S_subsamples)
yitp = f.(titp)
## -- Check against true arclengths between interpolated and original times
import ForwardDiff as FD
using QuadGK
norm_f′(t) = norm( FD.derivative(f, t) )
# integrate ||f'|| on each interval
S_itp = getindex.(quadgk.(norm_f′, titp[1:end-1], titp[2:end]), 1)
S_org = getindex.(quadgk.(norm_f′, t[1:end-1], t[2:end]), 1)
# maximum relative deviation
ΔS = step(S_subsamples)
maximum(x->abs(x-ΔS), S_itp)/ΔS
```