Cubic Spline Interpolation in Matlab

This was a homework in my geometric modeling class. Later I modified it to receive input with mouse click and a button to clear and redraw, and also a drop down to choose parametrization method.

Take a look at the spiral I made!!

Bezier representation:
$$
p(t)=\sum_{i=0}^{n} c_i B^n_i(t)
$$

function [ ] = draw( )
figure
hold on; box on;
grid;
uicontrol('Style', 'pushbutton', 'String', 'Try me again!',...
        'Position', [20 20 50 20],...
        'Callback', @draw_new);  
uicontrol('Style', 'popup',...
           'String', 'uniform|centripetal|chordal',...
           'Position', [20 340 100 50],...
           'Callback', @setmethod); 
draw_new ('','');
end

function [] = setmethod (hObj,event)
axis([0 10 0 10]);
global PARAMETER;
PARAMETER=get(hObj,'Value');
draw_spline();
end

function[]=draw_spline()
cla;
global XY;
global PARAMETER;
if PARAMETER ==1
    interpolate(XY, 0, 'r');
elseif PARAMETER == 2
    interpolate(XY, 0.5, 'g');
elseif PARAMETER == 3
    interpolate(XY, 1, 'b');
else
    interpolate(XY, 0, 'r');   
end
end

function [] = draw_new (hObj,event) %#ok
global XY;
cla;
axis([0 10 0 10]);
% Initially, the list of points is empty.
XY = [];
n = 0;
% Loop, picking up the points.
disp('Left mouse button picks points.')
disp('Right mouse button picks last point.')
but = 1;
while but == 1
    [xi,yi,but] = ginput(1);
    plot(xi,yi,'ro')
    n = n+1;
    XY(:,n) = [xi;yi];
    if (length(XY(1,:)) > 1)
        draw_spline;
    end
end
draw_spline();
end

function [] = interpolate(x, mu, color)
n=length(x);
t=parameterize(x,mu);
h(1:n-1) = t(2:n)-t(1:n-1);
mi=get_cubic_slope(x,t);
for i=1:n-1
    a=x(1,i);
    b=x(1,i+1);

    p(1,:)=x(:,i);
    p(2,:)=x(:,i)+(h(i)).*mi(:,i)/3;
    p(3,:)=x(:,i+1)-(h(i)).*mi(:,i+1)/3;
    p(4,:)=x(:,i+1);

    [q_bez] = deCasteljau(0,1,p,linspace(0,1,100));

    plot(q_bez(:,1),q_bez(:,2),[color,'-'],'LineWidth',3);
    plot(p(:,1),p(:,2),[color,'--o'], 'MarkerSize',3);
end

plot(x(1,:),x(2,:),'ro', 'MarkerSize',5,'MarkerFaceColor','r');
xlabel('x'), ylabel('y'),title(['Cubic Spline Interpolation']);
end

function[m] = get_cubic_slope(x,t)
%cublic spline interpolation
n=length(x);
h(2:n) = t(2:n)-t(1:n-1);
A=zeros(n-1,n-1);
b=zeros(2,n-1);
b(:,1)= 3*((h(1)/h(2)*(x(:,2)-x(:,1))));

for i=1:n-1
    A(i,i)=2*(h(i)+h(i+1));
end

for i=2:n-1
    b(:,i)=3*((h(i+1)/h(i)*(x(:,i)-x(:,i-1)))+(h(i)/h(i+1)*(x(:,i+1)-x(:,i))));
    A(i,i-1)=h(i+1);
    A(i-1,i)=h(i);
end

m=A\b';
m=m';
m(:,n)=0;
end

function[m] = get_hermite_slope(x,t)
n=length(x);
%piecewise hermite interpolation
m(:,1)=(x(:,2)-x(:,1))./(t(2)-t(1));
m(:,n)=(x(:,n)-x(:,n-1))./(t(n)-t(n-1));
m(1,2:n-1)=(x(1,3:n)-x(1,1:n-2))./(t(3:n)-t(1:n-2));
m(2,2:n-1)=(x(2,3:n)-x(2,1:n-2))./(t(3:n)-t(1:n-2));
end

function [t] = parameterize(x, mu)
t=zeros(1,length(x));
for i=2:length(x)
    t(i)=sqrt((x(1,i) - x(1,i-1))^2 + (x(2,i) - x(2,i-1))^2)^mu + t(i-1);
end
end

function [p_t] = deCasteljau(a,b,p,t)
n = size(p,1);
m = length(t);
p_t = zeros(m,2);
X(:,1) = p(:,1);
Y(:,1) = p(:,2);

for j = 1:m
    for i = 2:n
        X(i:n,i) = (t(j)-a)/(b-a)*X(i-1:n-1,i-1) + (b-t(j)-a)/(b-a)*X(i:n,i-1);
        Y(i:n,i) = (t(j)-a)/(b-a)*Y(i-1:n-1,i-1) + (b-t(j)-a)/(b-a)*Y(i:n,i-1);
    end

    p_t(j,1) = X(n,n);
    p_t(j,2) = Y(n,n);
end

end