We implement a simple spline scheme in Euler Math Toolbox (EMT) just to show how this can be done.

EMT already has a spline function which computes the natural cubic spline through data points. The natural spline is a cubic polynomial between the data points with smooth second derivative in the date points. In the first and last data point the second derivative is assumed to be 0. Thus the spline can be continued with a straight line outside of the interpolation range.

Let us demonstrate this by computing a spline through 4 points in the plane.

>P=[-1,-1;1,-1;1,1;-1,1;-1,-1]'; ...
>s1=spline(0:4,P[1]); s2=spline(0:4,P[2]); ...
>x=-0.1:0.01:4.1;  ...
>plot2d(splineval(x,0:4,P[1],s1),splineval(x,0:4,P[2],s2),r=1.5); ...
>plot2d(P[1],P[2],>points,>add):

The points cp from 0 to 4 are the control points. The spline is evaluated from -0.1 to 4.1 then. The function spline() is only for one dimensional splines. So we must use it for each component once. From the plot, we see that the spline is not periodic. There are other formulas for periodic spline which could be easily implemented in EMT in the same way.

There are more schemes to generate splines. One is to use the first derivative in the control points. This is also called a Catmull-Rom spline.

To implement this, we first make a function which evaluates a spline on points t=0 to t=m given the values and the derivatives in the control points. To get a formula for this, we solve the problem to find a cubic polynomial in [0,1] with given values y0, y1 and derivatives v0, v1 in 0 and 1 respectively.

>function s(t) &= y0+v0*t+a*t^2+b*t^3

                                       3      2
                        y0 + t v0 + b t  + a t

>sol &= solve([s(1)=y1,diffat(s(t),t=1)=v1],[a,b])

         [[a = 3 y1 - 3 y0 - v1 - 2 v0, b = - 2 y1 + 2 y0 + v1 + v0]]

With this information, we design the evaluating function.

>function sfdeval (x,y,v) ...
$  a=fold(y,[-3,3])+fold(v,[-2,-1]);
$  b=fold(y,[2,-2])+fold(v,[1,1]);
$  res=zeros(rows(y),cols(x));
$  k=floor(x); t=x-k;
$  for j=1 to cols(x);
$    m=k[j]+1;
$    if m>cols(a) then res[,j]=y[,-1];
$    else
$      res[,j] = y[,m]+v[,m]*t[j]+a[,m]*t[j]^2+b[,m]*t[j]^3;
$    endif;
$  end;
$  return res;
$endfunction
>plot2d("sfdeval(x,[0,1,0.5],[1,0,1])",0,2): // quick test for one dimension

It should only be applied in [0,m] for m+1 control points 0 to m. The test is only in one dimension. But the function already works in several dimensions. To compute the coefficients a and b in a quick way for all points, we use fold(). See the documentation for more information. Moreover, we already vectorize over a row vector x and return columns of points, one for each element in the vector x.

To compute the derivatives, we use fold() again. The derivatives in the points should satisfy

\(v_k = \dfrac{p_{k+1}-p_{k-1}}{2}\)

where the vector p is interpreted in a periodic way.

>function makeder (P) ...
$  return fold(P[,-1]|P|P[,1],[-0.5,0,0.5])
$endfunction
>P=[-1,-1;1,-1;1,1;-1,1]';
>V=makeder(P)
             1             1            -1            -1 
            -1             1             1            -1 
>P=[-1,-1;1,-1;1,1;-1,1]';x=0:0.01:4; Y=sfdeval(x,P|P[,1],V|V[,1]);
>plot2d(Y[1],Y[2]); plot2d(P[1],P[2],>points,>add):

Now the spline is periodic. Here is another example.

>function sfdcirc (x,P) ...
$  V=makeder(P);
$  return sfdeval(x,P|P[,1],V|V[,1]);
$endfunction
>P=[0,0;1,1;0,0;1,-1;0,0;-1,-1;0,0;-1,1]';  ...
>x=0:0.01:cols(P); Y=sfdcirc(x,P); plot2d(Y[1],Y[2],r=1.5):

Here the derivatives in the corners are 0 by definition of the spline. This makes the spline return in a sharp way from the corners.