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.

List all files

Suppose you want to list all files which match a given regular expression  in a directory along its subdirectories. Maybe you want to gather your training samples, images, text files or whatever.

Here I provide a small code snippet to quickly do all the dirty work and give you back a cell array with the full names.

function [files] = ListFiles(directory, regex)
% LISTFILES - returns all files located in dir and its subdirectories 
%             whose names matches the regex.
%  
% For example ListFiles('.', '.*.avi') will return all avi files in the 
% current directory and all of it's subdirectories. 
%%
    % check if directory is valid otherwise return
    if isempty(directory) | isnan(directory)
        files = {};
        return;
    end
    
    % get all files and subdirectories in directory
    allDirectoryEntries = dir(directory);
    
    % get files in current directory
    n = length(allDirectoryEntries);
    
    % init empty files list and directories list
    files = cell(n,1);
    directories = cell(n,1);
    k = 0;
    d = 0;
    
    for i = 1 : n     
        file = allDirectoryEntries(i);
        name = fullfile(directory, file.name);
        % if file is directory add to directories otherwise add to files
        if file.name(1) ~= '.' && file.isdir && ~strcmp(name,directory)
            d = d + 1;
            directories{d} = fullfile(directory, file.name);
        elseif length(regexp(name, regex)) > 0 && ~file.isdir
            k = k + 1;
            files{k} = name;
        end
    end
    
    % keep only the non empty cells of files
    files = {files{1:k}};
    
    % keep only the non empty cells of directories
    directories = {directories{1:d}};
    
    % run ListFiles subdirectories recursively in all subdirectories
    for i = 1 : d
        files = [files ListFiles(directories{i}, regex)];
    end
end

And a simple example to gather all images file in the inputDirectory and its subdirectories :

%------------------------------------------------------------------------
    % Find all image files in input directory
    files = ListFiles(inputDirectory,'.*.(jpg|jpeg|tif|tiff|png|gif|bmp|PNG)');
%------------------------------------------------------------------------

Gaussian 2d

I needed to compute a 2-dimensional Gaussian distribution which is very common when using Gabor filters.
First google result Custom 2D Gauss  provided a quick solution but upon first look the implementation didn't take advantage of any of matlab's features (i.e matrix manipulation) or included functions so it is a bit slow.

I ended rigging it up a bit making it very fast.

One problem I didn't address is adding a scale for the pixels but that can be fixed inside applications by scaling the standard deviation.

The code in Matlab:

function mat = Gaussian2d(...
    gsize, center, sigmax, sigmay, theta, offset, factor)
% GAUSSIAN2d - Produces a 2-dimensional gaussian matrix
%
%   Input
%   gsize   :  Size of the output 'gauss', should be a 1x2 vector
%   center  :  The center position of the gaussian [row col]
%   sigmax  :  Std. dev. in the X direction
%   sigmay  :  Std. dev. in the Y direction
%   theta   :  Rotation in degrees
%   offset  :  Minimum value in output
%   factor  :  Related to maximum value of output, should be 
%                    different from zero  
%
%   Output
%   mat     : probability vector for each angle bin 
%-------------------------------------------------------------------------
%%  Validate input arguments
    if ~exist('gsize','var')
        error('Error : gsize argument not set');
    end
    
    if ~exist('center','var')
        center = gsize / 2;
    end
    
    if ~exist('sigmax','var')
        sigmax = 1;
    end
    
    if ~exist('sigmay','var')
        sigmay = 1;
    end
    
    if ~exist('theta','var')
        theta = 0;
    end
     
    if ~exist('offset','var')
        offset = 0;
    end
    
    if ~exist('factor','var')
        factor = 1;
    end
    theta  = (theta/180)*pi;  
%-------------------------------------------------------------------------
%%  Compute
    [X,Y] = meshgrid(1:gsize(1),1:gsize(2));
    Xc  = X - center(1);
    Yc  = Y - center(2);
    xm  = (Xc)*cos(theta) - (Yc)*sin(theta);
    ym  = (Xc)*sin(theta) + (Yc)*cos(theta);
    u   = (xm / sigmax).^2 + (ym / sigmay).^2;
    mat = offset + factor * exp(-u/2);
end

The code in C# in the OpenCv framework :

        /// Create a rotated 2d Gaussian distribution 
        /// </summary>
        /// <param name="_width">Width of the result image in pixels</param>
        /// <param name="_height">Height of the result image in pixels</param>
        /// <param name="_centerX">X-coordinate of the center of the distribution in pixels</param>
        /// <param name="_centerY">Y-coordinate of the center of the distribution in pixels</param>
        /// <param name="_sigmaX">Standard deviation of X-axis in pixels</param>
        /// <param name="_sigmaY">Standard deviation of Y-axis in pixels</param>
        /// <param name="_theta">Rotation of the distribution in degrees</param>
        /// <param name="_offset">Resulting value offset</param>
        /// <param name="_factor">Resulting value multiplier</param>
        /// <returns></returns>
        public static Image<Gray, float> Gaussian2D(
            int   _width    = 10, 
            int   _height   = 10,
            float _centerX  = 5, 
            float _centerY  = 5,
            float _sigmaX   = 1.0f, 
            float _sigmaY   = 1.0f,
            float _theta    = 0.0f, 
            float _offset   = 0.0f, 
            float _factor   = 1.0f)
        {
            float[, ,] data = new float[_height, _width, 1];
            float xc, yc, xm, ym , tmpX, tmpY, u,
                cosTheta = (float)Math.Cos(_theta / 180.0f * Math.PI), 
                sinTheta = (float)Math.Sin(_theta / 180.0f * Math.PI);
            _sigmaX = 1.0f / _sigmaX;
            _sigmaY = 1.0f / _sigmaY;
            for (int i = 0; i < _width; i++)
            {
                xc = i - _centerX;
                for (int j = 0; j < _height; j++)
                {
                    yc = j - _centerY;
                    xm = xc * cosTheta - yc * sinTheta;
                    ym = xc * sinTheta + yc * cosTheta;
                    tmpX = xm * _sigmaX;
                    tmpY = ym * _sigmaY;
                    u = tmpX * tmpX + tmpY * tmpY;
                    data[j, i, 0] = _offset + _factor * (float)Math.Exp(-u / 2);
                }
            }
            return new Image<Gray,float>(data);
        }