% Disegno figura 1

clear all
close all
clc

puntiBezier = [0 1 2 1 0 1 2 1 0; 0 0 1 2 2.5 3 4 5 5]; % punti in uno spazio 3D
%puntiBezier = [0 1 1 0; 0 0 1 1]; % punti in uno spazio 3D
np = size(puntiBezier,2);       % numero di punti del poligono di controllo
dim = size(puntiBezier,1);      % dimensione
sudd = 100;                     % suddivisioni di [0,1]
t = linspace(0,1,sudd);         % t
curvaBezier = zeros(dim,sudd);  % succ. dei punti della curva di Bezier
% valutazione della curva di Bezier per ogni istante del parametro
for i=1:sudd,
    % valutazione del punto
    bt = deCasteljau(puntiBezier,t(i));
    % salvataggio del punto
    curvaBezier(:,i) = bt;
end

% disegno della curva come spezzata e del poligono di controllo
xBezier = curvaBezier(1,:);
yBezier = curvaBezier(2,:);

h = plot(xBezier,yBezier,'k');
set(h,'LineWidth',2);
hold on
h = plot(puntiBezier(1,:),puntiBezier(2,:),'ro--');
set(h,'LineWidth',2);
hold off


p = 3;
divisioni = np - p+1;
U = [repmat(0,1,p),linspace(0,1,divisioni),repmat(1,1,p)]
kmax = length(U)-p-1;
uvalue = linspace(0,1,101);
Nip = zeros(length(uvalue),kmax);
for k=1:kmax,
    for it=1:length(uvalue),
        u = uvalue(it);
        Nip(it,k) = BSplineEval(k,p,U,u);
    end
end

Cx = zeros(length(uvalue),1);
Cy = zeros(length(uvalue),1);
for k=1:kmax,
    Cx = Cx+Nip(:,k)*puntiBezier(1,k);
    Cy = Cy+Nip(:,k)*puntiBezier(2,k);
end
hold on;
h = plot(Cx,Cy,'b');
set(h,'LineWidth',2);
hold off;

% Valutazione dei punti di interruzione
uvaluei = linspace(0,1,np-p);
Nipi = zeros(length(uvaluei),kmax);
for k=1:kmax,
    for it=1:length(uvaluei),
        u = uvaluei(it);
        Nipi(it,k) = BSplineEval(k,p,U,u);
    end
end

Cxi = zeros(length(uvaluei),1);
Cyi = zeros(length(uvaluei),1);
for k=1:kmax,
    Cxi = Cxi+Nipi(:,k)*puntiBezier(1,k);
    Cyi = Cyi+Nipi(:,k)*puntiBezier(2,k);
end
hold on;
h = plot(Cxi,Cyi,'bo');
set(h,'LineWidth',2);
hold off;
axis([-0.1 2.2 -0.1 5.1])

