# Spline Class

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

## Define a Spline

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

s = SPLINE(T,n,P,[type])

• T = {T1,...} is a cell array containing the knot sequence for each dimension (up to two in the current version). For an easier application, the cell array brackets are only necessary for the bivariate case. Each Ti allows three formats
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.
• n is the order as 1x1 or 1x2 array.
• P are the coefficients. They use the same convention as the result of the evaluation: The first dimension corresponds to the space dimension, e.g. 1-,2-,3-dimensional coefficients, and the last dimensions correspond to the indices of operation elements, e.g. B-splines, evaluation points. It is possible to omit the first dimension for scalar spline surfaces (use a mxn matrix instead of a 1xmxn matrix as argument).
• the optional argument type is a character string that indicates if the spline is full periodic in all dimensions ('o','oo'), periodic only in a specific dimension ('-o','o-'), or non-periodic ('-','--'). The default cases are '-' and '--', respectively.

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 . 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:

• two row vectors representing the x- and y-components of the points,
• a column vector for the x- and a row vector for the y-components. In this case, the spline is evaluated at the grid points, as they are generated by the ndgrid command.

The optional argument 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.

• : Evaluation of the (partial) derivative of the spline at the given points. (Only if there are coefficients defined, otherwise see next section.)
• : For spline curves, the tangent vectors, normal vectors, binormal vectors, curvatures and torsions are calculated (binormal and torsion only for space curves).
• : As for spline curves, the first indicator calculates the normal vectors of the surface. The result of the second indicator is a struct that contains the 3x3 matrices E representing the embedded Weingarten map, the mean curvature (H), the Gaussian curvature (K) and the first (G) and second fundamental (B) form at the evaluation points.
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 that turns the spline into a general smoothing spline for given data. It possesses 4 arguments which correspond to

1. a cell array of cell arrays. Each of them has the form , with the same convention as for the evaluation. Additionally, it is possible to interpolate normal vectors and curvatures. Therefore, it requires entries and .
2. a cell array consisting of 3 or 4 entries . 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 is possible, which allows the separation of the points by their coordinates. (The same is possible for the interpolation points.)
3. a matrix with entries representing weights for the derivatives with for .
4. a scalar .

The algorithm replaces the coefficients of the Spline object while keeping the order and the knots such that it solves the following minimization problem    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

• approximation data with the same syntax as described in the previous section. The approximation is performed without any smoothing terms which may result in bad results for an inappropriate choice of points.
• a function handle or a cell array of function handles: the spline is created by performing an approximation of the function on the parameter interval. In case of a cell array, each component of the spline represents the corresponding function.
• a Spline object: this case has the same result as in the previous case but it allows the insertion or removal of knots. As can be seen in the following example, the knot vectors do not have to be presorted. Please check the last chapter for the syntax of the additionally available functions for knot insertion/removal.
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 in order to visualize it. Its syntax is similar to the evaluation and the output argument is the handle to the graphic.

• handle = Plot() plots the Spline in its domain at uniformly generated (grid) points.
• handle = s.Plot(C1,C2,...,Cn) uses the general plotting parameters (LineWidth, Color, ...) specified in C1,...,Cn similarly to the Matlab functions plot, plot3, surf.
• handle = s.Plot(k,C1,C2,...,Cn) plots the k-th derivative of a spline curve s, with k a 1x2 vector for surfaces. k can also be chosen as the characters known from the evaluation ('N','T','B') to additionally plot the corresponding vectors.
• handle = s.Plot(x,k,C1,C2,...,Cn) plots the k-th derivative of the curve s at the points specified in x. The call Plot(x) has the same result as Plot(x,0).
• handle = s.Plot(x1,x2,k,C1,C2,...,Cn) plots the k-th derivative of a spline surface s at points generated by the vectors x1 and x2. The evaluation points follow the convention as for the evaluation of a Spline object. Additionally, the same rules for k appplies as for a univariate spline. For k = 'E' the plots are coloured by the Gaussian and the Mean curvature.
• the parameter 'mesh', followed by a number m, plots a surface with m parameter lines per dimension. Alternatively, m can be a 1x2 vector to add more parameterlines in one Dimension, or m can be a 1x2 cell array that contains the specified grid lines. Additional plot parameters can be used but must be compatible with the command trisurf.

## 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):

• uniSubd(k) performs k uniform subdivision steps of the knot vector for a uniform spline curve.
• insertKnot(t0,[dim]) inserts the knot t0 in the [dim]-th variable of a spline. The second parameter is only required for bivariate splines.
• deleteKnot(index,[dim]) deletes the index-th knot in the [dim]-th variable of a spline. The second parameter is only required for bivariate splines.
• gramian(k,[dim]) returns the Gramian matrix for the k-th derivative of the basis functions. For a multivariate Spline, the second parameter can be set to get the Gramian matrix for the univariate basis functions in this specific dimension.
• [P,x,[y]] = sample(m,[k]) returns m uniform sample values P of the k-th derivative of the spline object at the points (x,[y]). m can be a scalar or an array for a bivariate spline.
• the well known transpose ' operator for matrices performs an analogue operation for a bivariate spline by interchanging the dimensions of the parameter domain.
• set_P(newP) sets the coefficients of a Spline object to the values specified in newP.
• interpolate(x,[y],F) interpolates the data points (x,[y],F), with x and y given as for the evaluation (gridded or non-gridded).
• schoenberg(f) calculates the Schoenberg approximation of the function handle f and overrides the current coefficients.

## 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:

• Type
• Knots
• Domain
• Order
• Coefficients
• Dimension

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:

• type: a string that indicates the type of a spline, starting with a 'c' for a curve or 's' for a surface, followed by up to two indicators '-' or 'o' as described above.
• ktype: a string for each dimension that indicates a uniform ('u') or nonuniform ('n') knot sequence.
• T: a cell array with the complete knot sequence for each dimension. Especially for periodic splines, the periodic interval is extended by additional knots at both sides to get a 'normal' knot sequence.
• n: the order for each dimension.
• P: the coefficients of the spline with the same extension as for the knots for periodic splines.
• h: the lengths of the intervals of the knot sequences.
• dim: the dimension of the spline space.
• grev: the Greville abscissae corresponding to the spline space.
• dom: the domain of the spline space.
• monoms: this property contains the different polynomial segments of the basis functions for each knot interval.
• local_gramian: this property contains the local Gramians, which are the Gramians reduced to each knot interval and only referring to the active basis functions.