clear all
close all

function [bx,by,bz] = bezier(p,t)
  [p1,p2] = size(p);
  [t1,t2] = size(t);
  n = p1-1;
  if p2 ~= 3
    error('p must be an n x 3 matrix');
  end
  bx = zeros(t1,t2);
  by = zeros(t1,t2);
  bz = zeros(t1,t2);
  o = ones(t1,t2);
  for k=0:n
    bx = bx + bincoeff(n,k)*t.^k.*(o-t).^(n-k)*p(k+1,1);
    by = by + bincoeff(n,k)*t.^k.*(o-t).^(n-k)*p(k+1,2);
    bz = bz + bincoeff(n,k)*t.^k.*(o-t).^(n-k)*p(k+1,3);
  end
end

figure
clear all

t=linspace(0,1,30);
p1 = [0,0,0; 0,0,1; 1,1,1; 1,0,0];
[g1x,g1y,g1z] = bezier(p1,t);

p2 = [0,1,0; 0,1,2; 1,0,1; 2,1,0];
[g2x,g2y,g2z] = bezier(p2,t);

hold on
plot3(g1x,g1y,g1z)
plot3(g2x,g2y,g2z)

% disegna superficie rigata di raccordo fra due curve
clear all
[t,s]=meshgrid(linspace(0,1,30),linspace(0,1,20));

% matrice di soli 1
o=ones(size(t));

p1 = [0,0,0; 0,0,1; 1,1,1; 2,0,0];
[g1x,g1y,g1z] = bezier(p1,t);

p2 = [0,1,0; 0,1,2; 1,0,1; 2,1,0];
[g2x,g2y,g2z] = bezier(p2,t);

x=(o-s).*g1x+s.*g2x;
y=(o-s).*g1y+s.*g2y;
z=(o-s).*g1z+s.*g2z;

figure
surf(x,y,z)

% aggiunge fattore correttivo (sempre uguale) sull'asse z
corr=1.5*(o-s).*s;

x1=(o-s).*g1x+s.*g2x;
y1=(o-s).*g1y+s.*g2y;
z1=(o-s).*g1z+s.*g2z+corr;

figure
surf(x1,y1,z1)

% aggiunge fattore correttivo sull'asse z, in funzione della distanza
% fra i punti da connettere
dist=sqrt((g1x-g2x).^2+(g1y-g2y).^2+(g1z-g2z).^2);

x2=(o-s).*g1x+s.*g2x;
y2=(o-s).*g1y+s.*g2y;
z2=(o-s).*g1z+s.*g2z+dist.*corr;

figure
surf(x2,y2,z2)
