## Sunday, March 25, 2012

### Cubic spline interpolation

It's been a long since I actually coded any interpolation method.
Matlab is notorious for making you lazy since it's so easy to get things done and you tend to stop looking under the hood. A friend asked me for help on a cubic interpolation problem and since that was too easy I expanded it so I can use it on my projects.

The math behind cubic spline is really simple. You piecewise fit cubic polynomials using 4 data values (two points and two tangents) in order to create a smooth spline that passes from all given points. The wikipedia sources are really good so I won't dive into the math. Instead I'll provide some matlab code for doing the dirty deed. Matlab (as always) has a command for this (spline) but we wont be using it because I like getting my hands dirty.
``````function [ yy, xx] = cubicSpline(N,points,gradients)
% CUBICSPLINE - returns N interpolation points using cubic spline
%   interpolation method, this method cannot be used for space curves
%
% Input
% N         : number of interpolation points
% points    : given points [x y]
% gradients : gradient at first and last point
%
% Output
% xx        : uniform spaced function interpolation at N points
% yy        : uniform spaced N points
%
%-------------------------------------------------------------------------
%% Validate input arguments
if isempty(N) || N < 1
error('N must be >= 1');
end

if isempty(points) || size(points,1) < 1
error('point must have >= 1 rows');
end

if isempty(points) || size(points,2) ~= 2
error('point must have 2 collumns');
end

error('gradients must have 2 elements');
end
%% coefficient calculation part

% get number of points
[rows ~] = size(points);

% compute inverse matrix to be used
matrix = inv([1 0 0 0 ; 0 1 0 0 ; 1 1 1 1; 0 1 2 3]);

% initialize coefficients structure
coefficients = zeros(rows-1,4);

% given n points we must calculate n-1 polynomials
for i = 2 : rows
pEnd = [];
pStart = [];

% calculate gradient using finite central differences
if (i-1) == 1
else
pStart = (points(i,2) - points(i-2,2))/2;
end

if i == rows
else
pEnd = (points(i+1,2) - points(i-1,2))/2;
end

% create vector [Pi P'i Pi+1 P'i+1]'
vector = [points(i-1,2);pStart;points(i,2);pEnd];

% calculate polynomial coefficients
coefficients(i-1,:) = (matrix * vector)';
end

%% interpolation part

% get max X and min X and interval
minX = points(1,1);
maxX = points(end,1);
intervalX = (maxX - minX) / (N - 1);
xx = minX : intervalX : maxX;

% interpolate at given locations
yy = zeros(1,N);
splineIndex = 1;
for i = 2 : N-1
x = xx(i);
% find the index of the used spline
for j = splineIndex : rows
if x >= points(j,1) && x < points(j+1,1)
splineIndex = j;
break;
end
end
splineCoeffs = coefficients(splineIndex,:);
% compute m
m = (xx(i) - points(splineIndex,1))/...
(points(splineIndex+1,1) - points(splineIndex,1));
% compute value with given spline and m
yy(i) = splineCoeffs(1) + splineCoeffs(2) * m + ...
splineCoeffs(3) * m^2 + splineCoeffs(4) * m^3;
end
yy(1) = points(1,2);
yy(end) = points(end,2);
end

``````

This code can be used to interpolate y=f(x) functions. For example :

``````%% Demonstration of cubic splines
N = 100;
x = [0:1:10];
y = sin(x);
xOriginal = [0:0.1:10];
yOriginal = sin(xOriginal);
gradient = [0 0];
[yy xx] = cubicSpline(N,[x' y'],gradient);
figure;
plot(xOriginal,yOriginal,xx,yy,'--',x,y,'o','LineWidth',2);
legend('Original function','Interpolation spline','Given points');
%%

``````

gives us the following graph :

any errors at the beginning and the end are due to the fact that I entered zero gradient at those points but provided the correct gradients the result should be much more precise.

To to make it somewhat useful in my projects I should use this function as a basis for calculating space curves. This excellent source explains that space curves are functions of u such as y = f(u) and x = f(u).

``````function [ yy,xx ] = cubicSpline2d(N, points, gradients )
% CUBICSPLINE - returns N interpolation points using cubic spline
%   interpolation method, this method cann be used for space curves
%
% Input
% N         : number of interpolation points
% points    : given points [x y]
% gradients : gradient at first and last point
%
% Output
% xx        : uniform spaced function interpolation at N points
% yy        : uniform spaced N points
%
%-------------------------------------------------------------------------
%% Validate input arguments
if isempty(N) || N < 1
error('N must be >= 1');
end

if isempty(points) || size(points,1) < 1
error('point must have >= 1 rows');
end

if isempty(points) || size(points,2) ~= 2
error('point must have 2 collumns');
end

error('gradients must have 4 elements');
end
%%
% get number of points
[rows ~] = size(points);

% get total length of points
u = [1 : rows]';
x = [points(:,1)];
y = [points(:,2)];

[xx,~] = cubicSpline(N, [u x],gradients(:,1));
[yy,~] = cubicSpline(N, [u y],gradients(:,2));
end

``````

and the test script

``````%% Demonstration of cubic splines 2d
u = 0:0.5:2*pi;
N = numel(u)*10;
y = sin(u);
x = sin(u) + cos(u);
gradient = [0 0; 0 0];
[yy xx] = cubicSpline2d(N,[x' y'],gradient);
plot(xx,yy,'--r',x,y,'ob','LineWidth',2);
legend('Interpolation spline','Given points');
title('Cubic space curve interpolation')
xlabel('x-axis');
ylabel('y-axis');
%%

``````

gives us the following plot :