How to create a spline for data without sorting?

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?

Answer 1

[John D'Errico]

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.

[David Winthrop]

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')