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

