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.

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

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

and the test script

gives us the following plot :

You can download the code here.

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
if isempty(gradients) || numel(gradients) ~= 2
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
pStart = gradients(1);
else
pStart = (points(i,2) - points(i-2,2))/2;
end
if i == rows
pEnd = gradients(2);
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
if isempty(gradients) || numel(gradients) ~= 4
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 :

You can download the code here.

## No comments:

## Post a Comment