function P = BezierCurve(Q, np)
  % -----------------------------------------------------------------------
  % BezierCurve - blend the n+1 points in Q according to 
  %                 P(t) = sum(i=0:n){Q(:,i+1)*C(n,i)*t^i*(1-t)^(n-i)}
  %               where
  %                 C(n,i) = n!/(i!*(n-i)!)
  %               and t goes from 0 to 1 in np-1 intervals
  %
  % Args: Q - 4 x n+1 matrix of Bezier control points in columns
  %       np - number of values of t in 0..1 at which to evaluate P(t)
  %
  % Rtns: P - 4 x np matrix of points on the curve
  % -----------------------------------------------------------------------

  n = size(Q, 2)-1; % the order of the polynomials (blending functions)
  N = BezierMatrix(n); % matrix combinatorial coefficients
  % row i of T is [t_i^n t_i^(n-1) ... t_i^0] where t_i = (i-1)/(np-1)
  T = repmat(linspace(0,1,np)',1,n+1) .^ repmat([n:-1:0],np,1);
  % each column of T * N is a blending function evaluated at np points
  P = T * N * Q';
end
