NaturalOpenCubicSplineCurves.m

%******************************************************
%***       NATURAL OPEN CUBIC SPLINE CURVES         ***
%***              Nguyen Quoc Duan                  ***
%******************************************************
% This is a practise exercise on CAD & Computer Graphics (CAO)

%  University of Liege - EMMC 
%  ( European Master in Mechanics of Construction )
%  Professor : Pierre Beckers   ( Aerospace Laboratory - ULg, Belgium )
%  Student   : Nguyen Quoc Duan ( EMMC11 - Ho Chi Minh City, Vietnam )

%  PROBLEM : *************************************************************
% Given a set of n points defined in the 2D Cartesian space.            **
% Build the open cubic spline curve using uniform spacing of knots      **
% (uniform parameterization)and natural spline definition.              **
% Build a second curve by using a new knot vector with spacing of knots **
% proportional to the distances between the points.                     **
% Draw the 2 curves and the points.                                     **
% ************************************************************************
clear all; close all; clc;


% ----------------------
% INPUT DATA 
%-----------------------
% Four data lines can be changed
n=30;            % number of points 
x=[1 5 3 1 3 5 7 7 7 9 9 9 11 11 13 13 13 15 15 15 17 17 19 21 23 21 19 19 21 23];    
y=[4 6 7 4 1 3 7 1 4 7 1 4 7   1  7  1  4 7  1  4  7  1  5  7  5  7  5  3  1  3];
uk=10 ;                 % number of uniform knots (-1) for a local segment


u=[0:1/uk:1]   ;        % uniform spacing of knots vector
lu=length(u);           % number (-1) of local knots    
% verify input data 
if length(x)==length(y)     
    n=length(x);
else
    disp('Unvalid values ! The number of x and y coordinates must be the same !')
    break
end

Y=[x' y'];      % matrix of point coordinates
m=n-1;              % number of segments
Ch=zeros(1,m);       % matrix of the length of the chords
for i=1:m
    Ch(i)=sqrt((Y(i+1,1)-Y(i,1))^2+(Y(i+1,2)-Y(i,2))^2);
end

if max(Ch)==0
    disp('Unvalid input data! There must be at least two different points!')
    break
else
    Chmax=max(Ch);      % the longest value of chords
end

% ----------------------    
%   [M]*[D]=[B]
% ----------------------
% Initial values 
M=zeros(n,n);   % three-diagonal matrix for natural cubic spline
B=zeros(n,2);   % right-hand side matrix
D=zeros(n,2);   % matrix of derivatives

% Compute values of M,B and D
for i=1:n       % loop for definition of M matrix
   if i==1 | i==n
       M(i,i)=2;
       M(1,2)=1;
       M(n,n-1)=1;
   else
       M(i,i)=4;
       M(i,i-1)=1;
       M(i,i+1)=1;
   end
end

for i=1:n       % loop for computing B matrix
    if i==1 
        B(i,:)=3*(Y(i+1,:)-Y(i,:));
    elseif i==n
        B(i,:)=3*(Y(i,:)-Y(i-1,:));
    else
        B(i,:)=3*(Y(i+1,:)-Y(i-1,:));
    end
end

D=M\B;          % solve equation M*D=B for D 

% ----------------------------------------------------
%           COMPUTE NODAL VALUES AT KNOTS 
%           with uniform spacing of knots 
% ----------------------------------------------------
U=zeros(lu,4);          % initial knot matrix
for j=1:lu              % loop for computing U matrix
    U(j,:)=[u(j)^3 u(j)^2 u(j) 1];
end

% Definition of Hermite matrix for cubic spline : P=U*N*G
 N=[2 -2 1 1;-3 3 -2 -1;0 0 1 0;1 0 0 0];   % numerical matrix 
 
 G=zeros(4*m,2);    % Initial Hermite geometry matrix
 P=zeros(m*lu,2);   % Initial nodal matrix of m segments
 for i=1:m
     for j=1:lu
         % Hermite geometry matrix of the i-th segment
         G(4*i-3:4*i,:)=[Y(i,:);Y(i+1,:);D(i,:);D(i+1,:)]; 
         % nodal matrix of the i-th segment
         P((i-1)*lu+j,:)=U(j,:)*N*G(4*i-3:4*i,:);
     end
 end

 % Delete repeat row to obtain matrix of nodal points 
 P(lu:lu:(m-1)*lu,:)=[];    
 
% ------------------------------------------------------------
%           COMPUTE NODAL VALUES AT KNOTS 
%      Spacing of knots proportional to the chords
% -------------------------------------------------------------
nuk=ones(1,m);      % matrix of numbers of uniform local spacing of knots
                    % of  the segments
integer=[1:1:uk];    
for i=1:m           % loop for computing number of spacing
    nuk(i)=uk*Ch(i)/Chmax;
    for j=1:uk              % loop for finding integer part of nuk(i)
        if integer(j+1)< nuk(i)
           j=j+1;
        else
            nuk(i)=integer(j);
            break
        end
    end
end

Pnuk=[];         % matrix of nodal values with non uniform knot 
for i=1:m
    NUK=[];                 % initial non uniform knot matrix
    up=[0:1/nuk(i):1];
    for j=1:length(up)      % loop for computing Un matrix
    NUK(j,:)=[up(j)^3 up(j)^2 up(j) 1];
    end
    % Hermite geometry matrix of the i-th segment
    iG=[Y(i,:);Y(i+1,:);D(i,:);D(i+1,:)]; 
    % nodal matrix of the i-th segment
    iPnuk=NUK*N*iG;
    if i==m
        iPnuk=iPnuk;
    else
        iPnuk(length(up),:)=[]; % delete the end knots of the (m-1) first
    end                         % segments.
    Pnuk=[Pnuk;iPnuk];
end

% -------------------------------------------------------
%                        PLOT GRAPHIC 
% -------------------------------------------------------
figure(1);
h=gcf;
set(h,'name','NATURAL OPEN CUBIC SPLINE CURVES');
axis equal;
plot(P(:,1),P(:,2),Pnuk(:,1),Pnuk(:,2),'r');
legend('uniform knots','proportional knots');
xlabel('x'),ylabel('y');
title('Uniform knots & knots proportional to the distances');
for i=1:(lu-1):length(P(:,1))
    text(P(i,1),P(i,2),num2str((i-1)/(lu-1)+1));
end

% --------------------------------------------------
%                   OUTPUT DATA
% --------------------------------------------------

disp('          NATURAL OPEN CUBIC SPLINE CURVES ');

disp('number of given points : '); n=n
disp('number of given segments : '); m=m
disp(' The coordinates of the given points :');
number=1:1:length(Y(:,1));
disp('point    x      y');
[number' Y]
disp('chord length :');
number=1:1:length(Ch);
disp('segment     length');
[number' Ch']
disp('Uniform knot vector: each segment is divided by :');
uk
disp('Non-uniform knot vector with spacing proportional to the distances.');
disp('Number of divided segments :');
number=1:1:length(nuk);
disp('segment   divided-into-pieces     length');
[number' nuk' Ch']
disp('The coordinates of all nodes with uniform knots:');
disp('  x      y');
P
disp('The coordinates of all nodes with non uniform knots:');
disp('  x      y');
Pnuk

disp('  THANK YOU VERY MUCH FOR YOUR INTERESTING LECTURES, PROFESSOR !');