interparc

Distance based interpolation along a general curve in space

Description

A common request is to interpolate a set of points at fixed distances along some curve in space (2 or more dimensions.) The user typically has a set of points along a curve, some of which are closely spaced, others not so close, and they wish to create a new set which is uniformly spaced along the same curve.

When the interpolation is assumed to be piecewise linear, this is easy. However, if the curve is to be a spline, perhaps interpolated as a function of chordal arclength between the points, this gets a bit more difficult. A nice trick is to formulate the problem in terms of differential equations that describe the path along the curve. Then the interpolation can be done using an ODE solver.

As an example of use, I'll pick a random set of points around a circle in the plane, then generate a new set of points that are equally spaced in terms of arc length along the curve, so around the perimeter of the circle.

theta = sort(rand(15,1))*2*pi;
theta(end+1) = theta(1);
px = cos(theta);
py = sin(theta);
 
100 equally spaced points, using a spline interpolant.

pt = interparc(100,px,py,'spline');

% Plot the result
plot(px,py,'r*',pt(:,1),pt(:,2),'b-o')
axis([-1.1 1.1 -1.1 1.1])
axis equal
grid on
xlabel X
ylabel Y
title 'Points in blue are uniform in arclength around the circle'

You can now also return a function handle to evaluate the curve itself at any point. As well, CSAPE is an option for periodic (closed) curves, as long as it is available in your matlab installation.

[~,~,foft] = interparc([],px,py,'spline');
foft(0:0.25:1)
ans =
      0.98319 0.18257
     -0.19064 0.98151
     -0.98493 -0.17486
      0.18634 -0.98406
      0.98319 0.18257

Acknowledgements

This file inspired Xy2sn.

• Required Products: MATLAB
• MATLAB release: MATLAB 7.12 (R2011a)
• Other requirements: Earlier releases will work as well.

[Evgeny Pr]

 

[Jon]

John,

Thanks for the response! I've been playing with the linear interpolation and you're right--it's much faster. It also seems to be accurate enough for my model. Thanks for your input!

[John D'Errico]

Jon - you can gain considerable speed by use of the 'linear' method in the code, if you can afford the loss in accuracy. It will be MUCH faster (as much as 1000 times faster), because I do the interpolation for that method in a different way from the other methods.

Those smoother methods in this tool are forced to call an odsolver for the interpolation, so the time taken is much higher for those methods. The only possibility is that I might have offered a different ode solver than ode45, but it is the fastest solver for these problems in my tests.

[Jon]

Excellent code.

I have only one complaint, which is the code takes awhile to run. For some applications this doesn't matter much, but I'm trying to resample a stream centerline after each time step in a 100's of time steps numerical model, and run times have become too large for me to use this (at least with each iteration). Do you have any suggestions for speeding it up? I tried going line-by-line through the code to find the bottleneck, but I was unsuccessful.

[John D'Errico]

Iain - You may find pchip a better choice for that class of curve. However, you have a valid point that there is a problem for this curve when a spline is used. I've just sent in a new version that should do better. It should appear in the morning.

[Iain]

Hi John. It is only a problem when pushing the spline fit quite hard, i.e. when there is a big gap between data. Then t and ti are not necessarily close. Test data and program at: http://www.srcf.ucam.org/~icw26/interparc_temp/

Try with and without change suggested. Thanks again for the code, hope that helps.

[John D'Errico]

Iain - correct, there was a problem in recognizing contractions. I accidentally used strcmpi, not strncmpi. Fixed, and uploaded to appear when they process it.

As far as derivatives, there is no need to move the derivative evaluation into the loop. ppval is fully vectorized, so no explicit loop is needed for that part. In fact, putting it inside a loop will just slow things down. If you have found a case where this fails please let me know, but the code does run correctly in my tests, as I would expect.

[Iain]

Excellent function, just what I needed. Thanks.

Slight bug in evaluating dudt - the (nargout > 1) block should be moved into the i=1:nt loop, and evaluated at ti, i.e.

for i=1:nt
%stuff - then
if nargout > 1
for j = 1:ndim
dudt(i,j) = ppval(spld{j},ti);
end
end
end
Also does not accept contractions of 'spline' etc. as equivalent arclength function does.

[John D'Errico]

Prabhu - I don't know what you are doing to generate nans. Please send me an example that causes it to happen, as none of my tests generated NaNs.

With one exception. There is actually one case I can find where interparc returns NaN. That is when I passed in a vector with NaNs in it already. In that case, how would one expect it not to return NaN? What logical action should it do when one asks it to interpolate a curve that is undefined at some set of points?

[Prabhu]

John:
Thanks for this excellent code.
I'm facing one problem when I use this code. It works fine in most cases. But sometimes its churning NaNs as one of the points. Can you help me out with this issue?
Thanks