# Interpolate across arc length for constant distances between points

**URL:** <https://discourse.julialang.org/t/interpolate-across-arc-length-for-constant-distances-between-points/104520>\
**Category:** New to Julia\
**Tags:** question, interpolations\
**Created:** [October 3, 2023, 1:56am UTC](https://discourse.julialang.org/t/interpolate-across-arc-length-for-constant-distances-between-points/104520 "2023-10-03T01:56:01Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![modifiedbear](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/modifiedbear/32/50494_2.png) [@modifiedbear](https://discourse.julialang.org/u/modifiedbear)\
**Post date:** [October 3, 2023, 1:56am UTC](https://discourse.julialang.org/t/interpolate-across-arc-length-for-constant-distances-between-points/104520/1 "2023-10-03T01:56:01Z")

</div>

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.

```matlab
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.

 ![ellipse_ds](https://global.discourse-cdn.com/julialang/original/3X/4/1/4156b4a5ca256172421df1e44826b67382bad1ca.png)  
The differences between arc length in both discretizations.

---

<div class="post-metadata">

**Author:** ![skleinbo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skleinbo/32/36080_2.png) [@skleinbo](https://discourse.julialang.org/u/skleinbo)\
**Post date:** [October 5, 2023, 2:47pm UTC](https://discourse.julialang.org/t/interpolate-across-arc-length-for-constant-distances-between-points/104520/2 "2023-10-05T14:47:07Z")

</div>

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](https://github.com/dextorious/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.

```julia
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

```