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