MATLAB Answers

0

How to create a spline for data without sorting?

Asked by David Winthrop on 15 Apr 2019
Latest activity Commented on by David Winthrop on 16 Apr 2019
Hello,
I have some data that represents trajectories through the 2D plane. My requirements:
  1. I need a pchip type spline to interpolate between the points. The default MATLAB function first sorts the points and this destroys the trajectory.
  2. I need to be able to divide up the splined path between successive points along the path into a given number of equal spaced segments.
Here are my points and the number of equal spaced points I need between each par (including the endpoints):
x = [218.0264 217.5283 217.4984 217.5386 217.1662 217.0156 218.0268 218.3743 218.3370 218.1877 218.0228 218.0729];
y = [276.9834 278.3856 278.6099 279.1224 280.0742 280.3255 279.0974 278.1296 277.3735 278.1461 275.6485 276.5042];
C = [13 13 15 19 29 31 9 11 3 35 31];
I found the excellent code by John D'Errico (INTERPARC) will solve the first requirement. The problem is that while it does fit a smooth spline between the points, and preserve the order, I can't figure out how to generate the equal spaced points between each pair. Rather, his program generates equal spaced points along the entire curve. So I need the curve his code generates if it passes through the dat points, but I need the equal spacing to include the x and y data as part of the sampling by including equal spaced pairs in between each.
Here is a plot of the spline fit returned by INTERPARC with 800 equal spaced points. I need to use this curve, fit to the entire data set, but I need to be able to segment between each pair of points along the way from point 1 to the final point according to vector C. So I need points along this curve but 13 points between (x(1),y(1)) and (x(2),y(2)) that include the endpoints, and so for each pair.
Any idea how I can accomplish this?

  0 Comments

Sign in to comment.

1 Answer

Answer by John D'Errico
on 16 Apr 2019
Edited by John D'Errico
on 16 Apr 2019
 Accepted Answer

I can certainly do it with interparc. But it gets tricky, because interparc is set up to find points that are designated in terms of given arclengths along the curve itself. Simplest is to just understand how to build an interpolant of this sort, then to use it directly.
Interparc uses a linear chordal arclength to parameterize the spline that it builds. I think that idea came from Eugene Lee, as I recall. Pick some random points to interpolate.
x = rand(1,5);
y = rand(1,5);
plot(x,y,'o-')
Now, suppose I decide to insert exactly 3 intermediate points between each (x,y) pair? First, compute a cumulative linear chordal arclength. Think of it as the distance along the line, as computed by a connect the dots interpolant.
t = [0,cumsum(sqrt(diff(x).^2 + diff(y).^2))];
>> t
t =
0 0.20299 0.53787 1.1588 1.5837
>> arclength(x,y,'linear')
ans =
1.5837
In fact you can see that my arclength function returns the same final number for the total arclength, when the "linear" method is used. So each of those arclength numbers are the cumulative lengths of the line segments connecting the points.
Now, I'll build a set of splines, representing x(t) and y(t).
x_t = pchip(t,x);
y_t = pchip(t,y);
Now, suppose I want to see 3 intermediate points between each pair. (I've chosen a small number here because I want to be easily able to count those points.)
nint = 3;
tint = cumsum([0,repelem(diff(t)/(nint+1),nint+1)]);
xint = ppval(x_t,tint);
yint = ppval(y_t,tint);
plot(x,y,'-ok')
hold on
plot(xint,yint,'--r*')
So 3 additional intermediate points between each of the original points, using a pchip interpolant on a completely general set of points. To get 13 points, just change nint to 13 instead of 3.
They don't have the property of being equally spaced in arclength along the curve, as interparc would generate. I know how to get interparc to do the work too, but it is late at night, and I'm feeling lazy if I can avoid actual thought.
Finally, see that I used repelem to build tint there. Repelem is a nice tool, introduced moderately recently in MATLAB. (In R2015a, according to the docs.) I could write it without repelem, but then it is late at night, as I said, and that one line seems elegant as it is.

  1 Comment

Thanks, John.
The thing I was missing was the (in hindsight) obvious parameterization of the splines by the distances. Here is how it looks in final form, for anyone else needing to do this in the future.
x = [218.0264 217.5283 217.4984 217.5386 217.1662 217.0156 218.0268 218.3743 218.3370 218.1877 218.0228 218.0729];
y = [276.9834 278.3856 278.6099 279.1224 280.0742 280.3255 279.0974 278.1296 277.3735 278.1461 275.6485 276.5042];
plot(x,y,'o-')
t = [0,cumsum(sqrt(diff(x).^2 + diff(y).^2))];
x_t = pchip(t,x);
y_t = pchip(t,y);
tp = linspace(t(1),t(end),1000);
hold on
plot(ppval(x_t,tp),ppval(y_t,tp),'r')
C = [13,13,15,19,29,31,9,11,3,35,31];
PTX = [];
PTY = [];
for k = 1:11
tn = linspace(t(k),t(k+1));
P = interparc(C(k),ppval(x_t,tn),ppval(y_t,tn));
PTX = [PTX;P(:,1)];
PTY = [PTY;P(:,2)];
end
PTX(cumsum(C)) = [];
PTY(cumsum(C)) = [];
plot(PTX,PTY,'*k')

Sign in to comment.