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
