version 1.3.0.0 (6.81 KB) by
John D'Errico

Distance based interpolation along a general curve in space

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

John D'Errico (2021). interparc (https://www.mathworks.com/matlabcentral/fileexchange/34874-interparc), MATLAB Central File Exchange. Retrieved .

Created with
R2011a

Compatible with any release

**Inspired:**
RivMAP - River Morphodynamics from Analysis of Planforms, Transform velocity field into wall coordinates, Active Figure Zoom for Selecting Points, jdugge/xy2sn, Probabilistic Earthquake Location in 2D

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Hami Bahadir BilirEge ErdemStefano PizziraniLuke HellwigMelonKingrubindanZITONG WANGsorry, please ignore my comment

ZITONG WANGdoes not work for ellipse: try to multiply sin(theta) in the example by 2

ZITONG WANGas per my comment, this can be easily seen by mulitple sin(theta) in the example by 2

SGBen CogswellActually works!!

NeelSushant PowarAngelina FrankThis function is so great, thank you!

Brandon LoboSaved me the time writing a function to smoothen an airfoil leading edge.

Brian FitzGibbonWorked perfectly for my problem! Thanks for the great function!

icylabThanks for this great toolbox. Is it possible to interpolate variables which contain NaNs? Essentially I have a variable xyz (x and y are coordinates) and z (depth) which is unevenly spaced, with very large gaps as NaNs. I want to interpolate the xyz variable onto a common interval spacing using your function but how can I deal with NaNs? Thank you again!

Dominik MattioliIs there an output for an equation of the produced curve, i.e. coefficients?

Matt JVery useful!

Stephen CobeldickA truly outstanding tool: very well conceived, very nicely implemented, very useful.

I would give this six star rating, if it were possible.

sri harsha arigelaIslam ElgamalElisabeth BöhmwalderWorks great with spline if you don't have the same values of data after each other. Otherwise "The data sites should be distinct."

DavidIs there any way to add endslope constraints (like in spline)?

Shuoqi ChenVibhav GaurTom DagesError raised if the spline is defined by only 3 points: spline returns the 1x1 struct rather than a cell of size 1, so spl{i} fails, although spl.coefs is what we want in this case.

Cell contents reference from a non-cell array object.

Error in distance2curve (line 226)

spl.coefs = [zeros(1,4-nc),spl{i}.coefs];

Ying TangShenbao ChenEdward BarnardThere is an error with the csape interpolation. Line 247 should read `n = n + 1` not `nt = nt + 1`.

Sathyanarayan RaoMehmet OZTURKRanjan MelpalJackie Zhaothanks for your code. Solved my problems.

Dominik RoszkowskiShan DououcKendra OudykHi John,

Thanks for posting this function!

I can't seem to get it to work properly for my data - the points are not evenly spaced along the interpolated line.

t=20;

px=[59150;58949;58822;58866;59115;59557];

py=[1182063;1178757;1175719;1172948;1170419;1168179];

pt= interparc(t,px,py,'pchip');

plot(px,py,'r*',pt(:,1),pt(:,2),'b-o')

Any idea what I'm doing wrong?

Thanks!

t4 t3Maral AmirAndrew DeFattaI am using this code to interpolate between experimental data points along a closed curve(aka not a single value function) so that I can match the size of data points in my theoretical closed curve.

Petr NechaevTakuma KitanishiAli MouradDear John

Thank you so much for your function. It is really useful and helpful for my work.

I have a question please.

how to get the the coefficients for each curve that pass through 2 points. i use cubic spline interpolation.

Thanks in advance

Tianyou ZhengI have to say thanks. Your code is helpful for me.

CarineNguyen Duy Minh PHANDear John

Thank you so much for your function. It is really useful for my work. Could you please give me the more detailed description about your interpolation method? In fact, Iam new in Matlab and also in the interpolation, so I don't understand completely the principe of your method. How can I cite your work if I use your function?

Nguyen Duy Minh PHANWillibald BremsSebastian WolffMarcus LööwBeatiful piece of code. Exactly what I needed. Thank you so much for this!

Daniel MorenoYour function saved my day! Thank you so much =]

Alexandr BuyvalJohn,

Thank you for this function. It exactly what I need.

Could you answer me how to get a curvature values along curve?

John D'Errico:)

Toby WenzThis is my favorite external Matlab function. I use it all the time in my research to resample datasets in order to make sure I can combine, plot and analyse them together.

John D'ErricoIf for some x-values, you may or may not get multiple y values, then the problem becomes complicated.

You may want to break the problem down into separate sub-curves, where the branch of the curve is chosen so the function is now single valued on that domain. Then you can use any tool.

Jianwei TuJohn,

Thank you for your reply. I actually have a relationship where some x have multiple y values. I was looking for a function similar to interp1 that does:

yinterp = thefunction(x,y,xgiven);

but I guess it does not make sense and could not be done according to what you just said.

Thanks anyway!

Jianwei

John D'ErricoJianwei,

If you have a problem with a single valued function of x, then you can just use a call to spline, pchip, or even just interp1. Just pass in the x values you want to interpolate at, and those tools can predict a result nicely.

If your relationship has multiple y values for any given x, then the above tools won't work. But then it won't make sense to be able to interpolate as a function of x then either.

So I'm not really sure why you would need interparc here.

Having said that, I suppose I can come up with a scenario. Consider a path in three dimensions, of the form

(x,y(x),z(x))

So y and z are both functions of x. Now, you might wish to interpolate points along that curve for given values of x. But that is still easy to do. Again, just use spline models for each of y(x) and z(x).

So I'm not sure what you might be looking for.

Jianwei TuJohn, is it possible to use interparc to generate a series of data points with user defined x values? Thanks, Jianwei

Mario MarínJoe PassmanJackAnother terrific and very useful routine from John D'Errico!

NadimI need a similar function that runs on Fortran. Does anyone know where I can find such a function. Thanks!

Kevin J. DelaneyFantastic! Just what I needed.

Christopher ThissenVery helpful - thanks!

LeonardoDear John,

thank you for your code. I would like to ask you..

It is in your plan to implement an upgrading also for the interpolating method "cscvn()" for closed curves ?

Thank you in advance.

navid alemiThanks alot.

HarelMaximilianWorks perfect! Thanks!

JustinThis is amazing. Basically matlab magic.

Running on MATLAB 7.8(R2009a)

Stevevery nice, works great for my 2D-points!

ThomasVery good code BTW.

John D'ErricoIt is never obvious, because the error often results from some language construct of which the older MATLAB release has no understanding. So it throws some strange error.

Blaisehey,

I have upgraded to 2014A and that fixes it. sorry to bother you, should have checked this myself

blaise

John D'ErricoEvery time someone tells me that some code of mine fails to run on their machine with a strange looking error, especially when they have called it in a normal way, my first guess is they are trying to run it with too old of a MATLAB release.

INTERPARC was written on MATLAB 2011a. If your release is significantly older than that, I imagine this is the problem. If not, then I can only guess you have some conflicting code on your search path, perhaps you wrote something called spline.m?

Blaisehello,

I have just downloaded the code and tried to run it with the code taken from above:

theta = sort(rand(15,1))*2*pi;

theta(end+1) = theta(1);

px = cos(theta);

py = sin(theta);

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

however i get the error message:

??? Reference to non-existent field 'coefs'.

Error in ==> interparc at 317

nc = numel(spl{i}.coefs);

Error in ==> test at 13

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

Just before executing the line spl is a 1*2 cell array, with the second cell being empty and the first one is a structure (with one field, dim=1).

I was surprised since the code looks working well for everyone else.

Anybody knows what goes wrong?

Best,

Blaise

John D'ErricoAlias, sorry. Never published, not by me at least. The idea is a basic one though.

1. Compute the cumulative linear arclength (t) along the curve. This is simply the distance along the curve using a connect-the-dots (piecewise linear) function.

2. Fit an interpolating spline through the points in each dimension, as a function of that cumulative linear arclength. This is the classic approach to fitting a curve through a completely general set of points in any number of dimensions, especially if that curve may be a closed one. Thus I generate the curves {x(t), y(t), z(t),...}. Many tools in MATLAB use this same scheme to fit interpolating splines through arbitrarily shaped curves.

3. Differentiate those spline functions, to get {x'(t),y'(t),z'(t), ...}. Since the functions are simply cubic polynomial segments, differentiation is trivial.

4. Integrate the cumulative arclength differential element,

sqrt(x'^2 + y'^2 + z'2 + ...)

Numerical integration is employed in this step, in this case an ODE solver is used, because ODE45 has the ability to integrate a function, finding where that function crosses a list of points of interest. That is important for interparc, since it must locate points along the curve at specific distances from the start point.

Merzouk YOUNSIIt's a good work Mr. John, this is exactly what I need. But can you give us some papers/references about the used interpolation methods ?

PeteSuperb.

*Extremely* useful, especially in conjunction with distance2curve.m and/or slmengine.m

I owe you many beers(!)

MayarThank You Mr.John, you found my problem very well.

With regards

John D'ErricoMayar -

I'm sorry, but you need to provide an example that fails. In fact, if I try to emulate what you SAID you did that fails, I have no problems at all.

x = [200 221 221 240];

y = [100 102 107 110];

interparc(10,x,y)

ans =

200 100

205.69 97.159

211.63 94.893

217.88 94.532

220.9 99.706

220.98 106.06

222.27 112.2

228.01 114.09

234.15 112.49

240 110

It is true that if there are two consecutive points that are identical for all the variables, then it will fail. That is, this sequence of points will cause an error:

x = [200 221 221 240];

y = [100 102 102 110];

If this is indeed your issue, and NOT what you actually stated, just drop the replicates of that point! The error message indicated exactly that. Learn how to use unique to solve that, or you can find my consolidator tool on the file exchange to help with floating point near replicates.

MayarDear:

this function is not working while we have two value for the same station. for instance i have two value of Z for the same station of x. ex: x=221.0 , z= 102 next value in the data set is x=221.0 , z= 107,

my data set contain such data, while i run this function, it give me error massage of " the data sites should be distinct", for t=1, and while i change t to vector suppose t=0:5:556 ( 556 end of data set), it give me again error of " all element of t must be 0<= t <=1 " . even I set the t as t=0:0.5:556.

wish you give me some instruction soon

with regards

John D'ErricoThere is no decoupling of the axes. Basic 2nd semester calc, or something like that. I won't claim it is obvious at all, since so many nice tricks in mathematics are only obvious after you see the trick. The arc length of a parametric curve is an integral of the sqrt of the sum of the squares of the first derivatives.

http://en.wikipedia.org/wiki/Arc_length

Since interparc wants to locate a specific point along that curve, it simply integrates that function using an ODE solver. A nice thing about the ODE solvers in MATLAB is they can identify the location where the integrated function crosses a point of interest.

The only decoupling of the axes is the use of splines on each dimension independently as a function of the linear chordal arc length. This allows me to define a parametric curve in an arbitrary number of dimensions, (x(t),y(t),z(t),...), and then form differentials at any point along that curve.

So we have the arc length as:

L = int(sqrt(x'^2 +y'^2 + ...))

where the differentials are with respect to t, the linear chordal arc length.

MathewThis program is really impressive. I wanted to understand the math a bit more since it seems to be unique relative to many other interpolation methods. From what I can understand for the spline case, is it using a cubic spline interpolation over each dimension independently, then using differential equations to solve for the length of the curve for each axis? For a curve in 3 dimensions, xyz, I am confused if it is decoupling each axis or it is solving them together.

If anyone understands how this actually works it would be great to hear.

DineshThanks so much John. It really helped me in my hard time. This program is very nice, it is fast also. Even if I input smaller number of points it adjust them to the required number of equally spaced points.

Thank you again.

Dinesh

Evgeny PrJonJohn,

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'ErricoJon - 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.

JonExcellent 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'ErricoIain - 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.

IainHi 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'ErricoIain - 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.

IainExcellent 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'ErricoPrabhu - 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?

PrabhuJohn:

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