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