Spline Class

This class allows an intuitive definition of spline curves or surfaces and the calculation of important differential geometric characteristics.

Contents

Define a Spline

The definition of a spline with this class requires at least 3 parameters:

s = SPLINE(T,n,P,[type])
  1. Ti = [a b m]', Ti = [a b m r]', creates a uniform knot sequence with m intervals between a and b, and uniformly distanced outer knots . If r is set (as number) the knots a and b are taken r-times and it might be the case that the domain is not [a,b] (optimal choice r = order).
  2. Ti = [t1 ... tm], creates a (non)uniform knot sequence with the knots specified in Ti and a parameter interval according to the order.
  3. Ti = 'b', is equal to Ti = [0 1 1 order] and generates Bernstein polynomials as basis functions represented as B-splines with multiple knots at 0 and 1.

Either the order or the coefficients have to be passed to the class. For a missing order, it is calculated from the coefficients and for missing coefficients, a spline without coefficients is generated. Only for a periodic spline the order is required.

T = [2 7 5]'; % or T = [-1 0 1 2 3 4 5 6 7 8 9 10]
c = rand(2,8);
s = Spline(T,[],c); % Spline curve of order 4
figure
subplot(2,2,1);
s.Plot();
s = Spline({[0 1 2 3 4]},3,[1 0 -1 0],'o'); % Periodic Spline curve of order 3
subplot(2,2,2);
s.Plot(linspace(-4,8));

T = {[1 2 1]' [3 11 8]'};
C = rand(3,2,11);
s = Spline(T,[],C); %3-dimensional Spline Surface of order (2,4)
subplot(2,2,3);
s.Plot();
C = zeros(3,4,4);
for i=0:3
    for j=0:3
        C(:,i+1,j+1) = [i/3;sin(i*pi/6)*sin(2*j*pi/3);sin(i*pi/6)*cos(2*pi*j/3)];
    end
end
s = Spline({[0 1 1 2 3 3 4] [0 4 4]'},[3 4],C,'-o'); %3-dimensional Spline Surface of order (3,4) that is periodic in the second dimension
subplot(2,2,4);
s.Plot();

Evaluation

For a simplified evaluation, the bracket notation of a spline object has been overloaded. Instead of calling the evaluation method, it is sufficient to write $$ {\tt Spline\_object(Pts1,(Pts2),[Ind])} $$. The first argument (for a curve), or the first and second (surface case) arguments, represent the evaluation points. In the second case, there are two possibilities:

The optional argument $${\tt ind} $$ either indicates the order of the derivative of the spline or a differential geometric characteristic. If this argument is left out, the spline itself is evaluated.

close all
x1 = linspace(1.1,3,100);
x2 = linspace(0,4,100);
yf = s(x1',x2,[0 0]); %Evaluation at the grid generated by ndgrid(x1,x2)
E = s(x1',x2,'G');

% plot of the surface colored by its Gaussian curvature and the normal
% vectors
surf(squeeze(yf(1,:,:)),squeeze(yf(2,:,:)),squeeze(yf(3,:,:)),E.K,'EdgeColor','none','LineStyle','none','FaceLighting','phong');
axis equal
x1 = linspace(1,3,20);
x2 = linspace(0,4,20);
yf = s(x1',x2,[0 0]);
N = s(x1',x2,'N');
hold on
quiver3(squeeze(yf(1,:,:)),squeeze(yf(2,:,:)),squeeze(yf(3,:,:)),squeeze(N(1,:,:)),squeeze(N(2,:,:)),squeeze(N(3,:,:)),'red');
colorbar;

Calculation of Collocation Matrices

Some applications require a fast calculation of collocation matrices which are the values of the B-spline basis functions or their derivatives at given (grid) points. For this purpose, this spline class permits the calculation of those matrices. In case of an already defined spline WITH coefficients, the curly bracket notation has been overloaded and an equivalent call as for the evaluation provides the matrix whose first dimension corresponds to the basis functions.

close all;
T = [2 7 6]';
x = linspace(2,7,10000);
c = rand(1,9);
s = Spline(T,[],c);
% evaluate all basis functions and their first derivative on the paramter interval
y = s{x,0};
dy = s{x,1};
plot(x,y(5,:),'blue');
hold on
plot(x,dy(5,:),'red');

Alternatively, these matrices can be generated by performing an evaluation of a spline object WITHOUT coefficients.

close all;
T = [0 0 1 1 4 4 6 9 10 10 10];
x = linspace(1,9,10000);
s = Spline(T,4,[]);
% evaluate all basis functions and their first derivative on the paramter interval
y = s(x,0);
dy = s(x,1);
plot(x,y(4,:),'blue');
hold on
plot(x,dy(4,:),'red');

By using the 'b'-notation, Bernstein polynomials as special case of B-splines can be created.

close all
x = linspace(0,1,10000);
b = Spline('b',4,[]); %creates all cubic Bernstein polynomials
y = b(x);
for i = 1:4
   plot(x,y(i,:));
   hold on
end

General Smoothing Spline

Each Spline object provides a method $$ {\tt Smooth} $$ that turns the spline into a general smoothing spline for given data. It possesses 4 arguments $$ ({\tt INTER}, {\tt APPROX}, {\tt SMOOTH}, {\tt alpha}) $$ which correspond to

  1. a cell array $$ {\tt INTER} $$ of cell arrays. Each of them has the form $$ (Pts_I1,(Pts_I2),Values_I,Order_I) $$, with the same convention as for the evaluation. Additionally, it is possible to interpolate normal vectors and curvatures. Therefore, it requires entries $$ (Pts_N1,(Pts_N2),Normalvectors_N,'N') $$ and $$ (Pts_E1,(Pts_E2),Normalvectors_E,kappa/E,'k'/'E') $$ .
  2. a cell array $$ {\tt APPROX} $$ consisting of 3 or 4 entries $$ (Pts_A1,(Pts_A2),Values_A,Weights_A) $$. The matrix for the weights has to be of the same size as the values, obviously only scalars instead of vectors, or it can be only one scalar that is applied to all points. Alternatively, an input of the form $$ (Pts_A1,(Pts_A2),(Values_A)_x,(Values_A)_y,(Values_A)_z,Weights_A) $$ is possible, which allows the separation of the points by their coordinates. (The same is possible for the interpolation points.)
  3. a matrix $$ {\tt SMOOTH} $$ with entries $$ (i,j) $$ representing weights for the derivatives $$ \partial^{i-1}_\alpha $$ with $$ \alpha = (i-1,0) $$ for $$ j = 1 $$.
  4. a scalar $$ {\tt alpha} $$.

The algorithm replaces the coefficients of the Spline object while keeping the order and the knots such that it solves the following minimization problem

$$ \min_f \left(\sum_{j=1}^{\#Pts_A} Weights_A^j |f(Pts_A^j)-Values_A^j|^2 + \alpha \int_{ParDom} \sum_\alpha (\partial_\alpha f)^2 \cdot SMOOTH(|\alpha|,\alpha) \right) $$

$$ s.t.\qquad \partial_{Order_I} f(Pts_I) = Values_I$$

$$ s.t.\qquad \partial^\perp f(Pts_N) = Normalvectors_N $$

$$ s.t.\qquad kappa/E f(Pts_E) = kappa/E $$

It should be mentioned that the interpolation of normalvectors and curvatures can only be applied for multi-dimensional spline curves or three-dimensional spline surfaces.

T = [2 8 20]';
s = Spline(T,4,[]);
INT = {{[2.5:5.5],[-1 2 -2 1],0}};
APPROX = {linspace(2,7,10),5*rand(1,10), 1};
SMOOTH = [0 0 1]';
alpha = 1/1000;
s = s.Smooth(INT,APPROX,SMOOTH,1/1000);
close all
s.Plot();
hold on
plot(APPROX{1},APPROX{2},'ored');
plot(INT{1}{1},INT{1}{2},'oblue');
s = s.Smooth(INT,APPROX,SMOOTH,1000);
s.Plot();

Create a Spline from Data/Function Handle/Spline

In addition to the methods described in the first section, this class provides some special ways to create a spline. All of them use the already presented syntax for the constructor with the difference that the coefficients are replaced by

close all
s1 = Spline({[0 5 10]'},3,@(x) exp(x)); % univariate approximation of exponential function
s1.Plot();
s2 = Spline({[s1.T{1} 2.75]},3,s1); % the knot t0 = 2.75 is inserted, no presorting is required
hold on
s2.Plot();

Plot a Spline

Each Spline object s provides a method $$ {\tt Plot} $$ in order to visualize it. Its syntax is similar to the evaluation and the output argument is the handle to the graphic.

Audi

This class supports the evaluation with audi-variables (https://de.mathworks.com/matlabcentral/fileexchange/61849-autodiff_r2016b) to apply further functions of that class to a spline (like the Laplace beltrami operator). The syntax is the same as for the usual evaluation with the exception that x (and y) have to be of type audi. The output always contains derivatives up to order 2, as it is the usual range for further applications. For that reason, the indicator for the derivative is not required for spline curves. However, it has been replaced for spline surfaces by a simple 0 or 1 value to indicate that the points are located on a grid (1) or not (0), so the programm can extract that grid and apply the faster evaulation routine.

Useful Methods/Commands

In addition to the functionalities discribed above, this class provides the following methods (arguments in []-notation are additional parameter):

Shown/Hidden Properties

Each class object possesses two types of properties. At first, the following properties are shown when the object is accessed over the workspace and are self-explanatory:

It is important to know that these properties can not be changed, as they are set in the constructor. The only exception are the coefficients, as they can be changed with the method set_P.

In addition to these properties, there are hidden ones that are used for internal calculations. By changing these, especially the last two of them, no proper work flow can be guarantueed as the only consistency check for matching order, knots and coefficients is made before the evaluation: