%
% MATLAB for DUMMIES 11-12
% Exercice 3
%
% Solution detaillee
%  Vincent Legat, Damien Bouvy, Denis-Gabriel Caprace, Gilles De Neyer, 
%  Kevin Degraux, Mohammad Iranmanesh, Christophe Lorant, Remy Moulin, 
%  Matthew Philippe, Karim Slaoui, Olivier Thiry, Laurent Van Eesbeeck,
%  Thomas Vanderstraeten, Antoine Weynants
%
% -------------------------------------------------------------------------
%
function matlab3()
[x y z] = canoe(0.05);

close all; figure;
ax = axes('Visible','off',...
    'PlotBoxAspectRatio',[450 350 450],...
    'DataAspectRatio',[1 1 1],...
    'CameraViewAngle',6);
view(ax,[-143 13]);  hold(ax,'all');
surf(x,z,y,'Parent',ax,'BackFaceLighting','unlit',...
    'FaceColor',[1 0 0],'LineStyle','none'); 
light('Parent',ax,'Position',[1 0 0]);
light('Parent',ax,'Position',[10 4 3]);

end

function [x,y,z] = canoe(dt)
%CANOE Provides data for drawing a canoe with "surf"
%   [X,Y,Z] = CANOE(DT) returns points on the surface of a canoe with 
%   spatial spacing DT. The canoe is computed using NURBS with second order 
%   B-splines.
%   

if (nargin == 0) dt = 0.05; end

%
% -1- Geometrical modelling
% -1.1- Parametrizing the cross-section of the canoe
%

S =  [0 0 0 1 1 2 2 3 3 4 4 5];
Xc = [0.7  0.7  0.5  0.2  0.0  -0.2 -0.5 -0.7 -0.7]; 
Yc = [0.0 -0.3 -0.7 -1.0 -1.0  -1.0 -0.7 -0.3  0.0];
Zc = ones(size(Xc));
Wc = [1 1 1 1 1 1 1 1 1];

%
% -1.2- Extruding the cross-section along the axis
%

T =  [0 0 0 1 1 2 2 3];
X = Xc' * [0.0 0.7 0.7 0.7 0.0];
Y = Yc' * [0.0 0.7 0.7 0.7 0.0];
Z = Zc' * [-1 -0.8 0.5 1.8 2];
W = Wc' * [1 0.5 0.5 0.5 1];

%
% -1.3- Bending both ends of the boat
%

Y = Y + ones(size(Yc'))*[0.3 0.0 0.0 0.0 0.6];

%
% -2- Generating all points 
%     Note the pre-allocation of Bt and Bs 
%       (nicely suggested by some tutors :-)
%       (not present in the reference program vase.m !)
%

p = 2;
nt = length(T) - 1;
t = [T(p+1):dt:T(nt-p+1)];
Bt = zeros(nt-p,numel(t));
for i=0:nt-p-1
  Bt(i+1,:) = b(t,T,i,p);
end

ns = length(S) - 1;
s = [S(p+1):dt:S(ns-p+1)];
Bs = zeros(ns-p,numel(s));
for i=0:ns-p-1
  Bs(i+1,:) = b(s,S,i,p);
end

w = Bs' * W * Bt;
x = Bs' * (W .* X) * Bt ./ w;
y = Bs' * (W .* Y) * Bt ./ w;
z = Bs' * (W .* Z) * Bt ./ w;

end

% 
% ===== B-Splines mini-library =====================================
%

function u = b(t,T,j,p)
%
% A simple (but not efficient) way of computing B-splines functions
%
i = j+1;
if p==0
    u = (t>= T(i) & t < T(i+p+1)); return 
end

u = zeros(size(t));
if T(i) ~= T(i+p)
    u = u + ((t-T(i)) ./ (T(i+p) -T(i))) .* b(t,T,j,p-1);
end
if T(i+1) ~= T(i+p+1)
    u = u + ((T(i+p+1)-t) ./ (T(i+p+1) -T(i+1))) .* b(t,T,j+1,p-1);
end
end