clear all
close all

% disegna curva (in qs caso circonferenza)
t=linspace(0,2*pi,20);
gx=cos(t);
gy=sin(t);
gz(1:length(t))=0;
plot3(gx,gy,gz)

figure
clear all

% disegna tubo ellittico
[t,phi]=meshgrid(linspace(0,2*pi,20),linspace(0,2*pi,20));
% il vettore 2d dim assume le dimensioni della matrice t
dim=size(t)
gx=cos(t);
gy=sin(t);
gz(1:dim(1),1:dim(2))=0;
nx=-cos(t);
ny=-sin(t);
nz=gz;
bx=gz;
by=gz;
bz(1:dim(1),1:dim(2))=1;
cn=0.2.*cos(phi);
cb=0.8.*sin(phi);
x=gx+cn.*nx+cb.*bx;
y=gy+cn.*ny+cb.*by;
z=gz+cn.*nz+cb.*bz;
surf(x,y,z)
axis equal

figure
clear all

% disegna tubo diverso
[t,phi]=meshgrid(linspace(0,2*pi,20),linspace(0,2*pi,20));
% il vettore 2d dim assume le dimensioni della matrice t
dim=size(t)
gx=cos(t);
gy=sin(t);
gz(1:dim(1),1:dim(2))=0;
nx=-cos(t);
ny=-sin(t);
nz=gz;
bx=gz;
by=gz;
bz(1:dim(1),1:dim(2))=1;
cn=0.6.*cos(phi).*(1+0.2*sin(2*phi).^2);
cb=0.6.*sin(phi).*(1+0.2*sin(2*phi).^2);
x=gx+cn.*nx+cb.*bx;
y=gy+cn.*ny+cb.*by;
z=gz+cn.*nz+cb.*bz;
surf(x,y,z)
axis equal

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);
p = [0,0,0; 0,0,1; 1,1,1; 1,0,0];
[fx,fy,fz] = bezier(p,t);
plot3(fx,fy,fz)

figure
clear all

% superficie che interpola due curve di Bezier
[t,s]=meshgrid(linspace(0,1,20),linspace(0,1,5));
p = [0,0,0; 1,1,2; 3,2,2; 5,1,1];
[fx,fy,fz] = bezier(p,t);
q = [0,0,0; 1,-1,2; 3,-2,2; 5,-1,1];
[gx,gy,gz] = bezier(q,t);
x=fx+s.*(gx-fx);
y=fy+s.*(gy-fy);
z=fz+s.*(gz-fz);
surf(x,y,z)
axis equal


