function [p,w] = nrbeval(nurbs,tt)
% 
% NRBEVAL: Evaluate a NURBS at parametric points.
% 
% Calling Sequences:
% 
%   [p,w] = nrbeval(crv,ut)
%   [p,w] = nrbeval(srf,{ut,vt})
%   [p,w] = nrbeval(vol,{ut,vt,wt})
%   [p,w] = nrbeval(srf,pts)
% 
% INPUT:
% 
%   crv		: NURBS curve, see nrbmak.
% 
%   srf		: NURBS surface, see nrbmak.
%
%   vol		: NURBS volume, see nrbmak.
% 
%   ut		: Parametric evaluation points along U direction.
%
%   vt		: Parametric evaluation points along V direction.
% 
%   wt		: Parametric evaluation points along W direction.
%
%   pts     : Array of scattered points in parametric domain
% 
% OUTPUT:
%
%   p		: Evaluated points on the NURBS curve, surface or volume as 
% 		Cartesian coordinates (x,y,z). If w is included on the lhs argument
% 		list the points are returned as homogeneous coordinates (wx,wy,wz).
% 
%   w		: Weights of the homogeneous coordinates of the evaluated
% 		points. Note inclusion of this argument changes the type 
% 		of coordinates returned in p (see above).
% 
% Description:
% 
%   Evaluation of NURBS curves, surfaces or volume at parametric points along  
%   the U, V and W directions. Either homogeneous coordinates are returned
%   if the weights are requested in the lhs arguments, or as Cartesian coordinates.
%   This function utilises the 'C' interface bspeval.
% 
% Examples:
% 
%   Evaluate the NURBS circle at twenty points from 0.0 to 1.0
% 
%   nrb = nrbcirc;
%   ut = linspace(0.0,1.0,20);
%   p = nrbeval(nrb,ut);
% 
% See also:
%  
%     bspeval
%
% Copyright (C) 2000 Mark Spink 
% Copyright (C) 2010 Carlo de Falco
% Copyright (C) 2010, 2011 Rafael Vazquez
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.

%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

if (nargin < 2)
  error('Not enough input arguments');
end

foption = 1;    % output format 3D cartesian coordinates
if (nargout == 2)
  foption = 0;  % output format 4D homogenous coordinates 
end
   
if (~isstruct(nurbs))
  error('NURBS representation is not structure!');
end

if (~strcmp(nurbs.form,'B-NURBS'))
  error('Not a recognised NURBS representation');
end

if (iscell(nurbs.knots))
  if (size(nurbs.knots,2) == 3)
    %% NURBS structure represents a volume

    num1 = nurbs.number(1);
    num2 = nurbs.number(2);
    num3 = nurbs.number(3);
    degree = nurbs.order-1;

    if (iscell(tt))
      nt1 = numel (tt{1});
      nt2 = numel (tt{2});
      nt3 = numel (tt{3});

      %% evaluate along the w direction
      val = reshape (nurbs.coefs, 4*num1*num2, num3);
      val = bspeval (degree(3), val, nurbs.knots{3}, tt{3});
      val = reshape (val, [4 num1 num2 nt3]);

      %% Evaluate along the v direction
      val = permute (val, [1 2 4 3]);
      val = reshape (val, 4*num1*nt3, num2);
      val = bspeval (degree(2), val, nurbs.knots{2}, tt{2});
      val = reshape (val, [4 num1 nt3 nt2]);
      val = permute (val, [1 2 4 3]);

      %% Evaluate along the u direction
      val = permute (val, [1 3 4 2]);
      val = reshape (val, 4*nt2*nt3, num1);
      val = bspeval (degree(1), val, nurbs.knots{1}, tt{1});
      val = reshape (val, [4 nt2 nt3 nt1]);
      val = permute (val, [1 4 2 3]);
      pnts = val;

      p = pnts(1:3,:,:,:);
      w = pnts(4,:,:,:);
      if (foption)
        p = p./repmat(w,[3 1 1 1]);
      end

    else

      %% Evaluate at scattered points
      %% tt(1,:) represents the u direction
      %% tt(2,:) represents the v direction
      %% tt(3,:) represents the w direction

      st = size(tt);
      nt = prod(st(2:end));

      tt = reshape (tt, [3, nt]);

      %% evaluate along the w direction
      val = reshape(nurbs.coefs,4*num1*num2,num3);
      val = bspeval(degree(3),val,nurbs.knots{3},tt(3,:));
      val = reshape(val,[4 num1 num2 nt]);

      %% evaluate along the v direction
      val2 = zeros(4*num1,nt);
      for v = 1:nt
        coefs = reshape(val(:,:,:,v),4*num1,num2);
        val2(:,v) = bspeval(degree(2),coefs,nurbs.knots{2},tt(2,v));
      end
      val2 = reshape(val2,[4 num1 nt]);

      %% evaluate along the u direction
      pnts = zeros(4,nt);
      for v = 1:nt
        coefs = reshape (val2(:,:,v), [4 num1]);
        pnts(:,v) = bspeval(degree(1),coefs,nurbs.knots{1},tt(1,v));
      end

      w = pnts(4,:);
      p = pnts(1:3,:);
      if (foption)
        p = p./repmat(w,[3, 1]);
      end

      if (numel(st) ~= 2)
        w = reshape (w, [st(2:end)]);
        p = reshape (p, [3, st(2:end)]);
      end
    end

  elseif (size(nurbs.knots,2) == 2)
    %% NURBS structure represents a surface
  
    num1 = nurbs.number(1);
    num2 = nurbs.number(2);
    degree = nurbs.order-1;

    if (iscell(tt))
      %% Evaluate over a [u,v] grid
      %% tt{1} represents the u direction
      %% tt{2} represents the v direction

      nt1 = length(tt{1});
      nt2 = length(tt{2});
    
      %% Evaluate along the v direction
      val = reshape(nurbs.coefs,4*num1,num2);
      val = bspeval(degree(2),val,nurbs.knots{2},tt{2});
      val = reshape(val,[4 num1 nt2]);
    
      %% Evaluate along the u direction
      val = permute(val,[1 3 2]);
      val = reshape(val,4*nt2,num1);
      val = bspeval(degree(1),val,nurbs.knots{1},tt{1});
      val = reshape(val,[4 nt2 nt1]);
      val = permute(val,[1 3 2]);

      w = val(4,:,:);
      p = val(1:3,:,:);
      if (foption)
	p = p./repmat(w,[3 1 1]);
      end

    else

      %% Evaluate at scattered points
      %% tt(1,:) represents the u direction
      %% tt(2,:) represents the v direction

      st = size(tt);
      nt = prod(st(2:end));

      tt = reshape (tt, [2, nt]);

      val = reshape(nurbs.coefs,4*num1,num2);
      val = bspeval(degree(2),val,nurbs.knots{2},tt(2,:));
      val = reshape(val,[4 num1 nt]);


      %% evaluate along the u direction
      pnts = zeros(4,nt);
      for v = 1:nt
	coefs = reshape (val(:,:,v), [4 num1]);
	pnts(:,v) = bspeval(degree(1),coefs,nurbs.knots{1},tt(1,v));
      end

      w = pnts(4,:);
      p = pnts(1:3,:);
      if (foption)
	p = p./repmat(w,[3, 1]);
      end

      if (numel(st) ~= 2)
        w = reshape (w, [st(2:end)]);
        p = reshape (p, [3, st(2:end)]);
      end
        
    end

  end
else

  %% NURBS structure represents a curve
  %%  tt represent a vector of parametric points in the u direction

  if (iscell (tt) && numel (tt) == 1)
    tt = cell2mat (tt);
  end
  
  st = size (tt);
  
  val = bspeval(nurbs.order-1,nurbs.coefs,nurbs.knots,tt(:)');

  w = val(4,:);
  p = val(1:3,:);
  if foption
    p = p./repmat(w,3,1);
  end

  if (st(1) ~= 1 || numel(st) ~= 2)
    w = reshape (w, st);
    p = reshape (p, [3, st]);
  end

end

end

%!demo
%! srf = nrbtestsrf;
%! p = nrbeval(srf,{linspace(0.0,1.0,20) linspace(0.0,1.0,20)});
%! h = surf(squeeze(p(1,:,:)),squeeze(p(2,:,:)),squeeze(p(3,:,:)));
%! title('Test surface.');
%! hold off

%!test
%! knots{1} = [0 0 0 1 1 1];
%! knots{2} = [0 0 0 .5 1 1 1];
%! knots{3} = [0 0 0 0 1 1 1 1];
%! cx = [0 0.5 1]; nx = length(cx);
%! cy = [0 0.25 0.75 1]; ny = length(cy);
%! cz = [0 1/3 2/3 1]; nz = length(cz);
%! coefs(1,:,:,:) = repmat(reshape(cx,nx,1,1),[1 ny nz]);
%! coefs(2,:,:,:) = repmat(reshape(cy,1,ny,1),[nx 1 nz]);
%! coefs(3,:,:,:) = repmat(reshape(cz,1,1,nz),[nx ny 1]);
%! coefs(4,:,:,:) = 1;
%! nurbs = nrbmak(coefs, knots);
%! x = rand(5,1); y = rand(5,1); z = rand(5,1);
%! tt = [x y z]';
%! points = nrbeval(nurbs,tt);
%!
%! assert(points,tt,1e-10)
%!
%!test
%! knots{1} = [0 0 0 1 1 1];
%! knots{2} = [0 0 0 0 1 1 1 1];
%! knots{3} = [0 0 1 1];
%! cx = [0 0 1]; nx = length(cx);
%! cy = [0 0 0 1]; ny = length(cy);
%! cz = [0 1]; nz = length(cz);
%! coefs(1,:,:,:) = repmat(reshape(cx,nx,1,1),[1 ny nz]);
%! coefs(2,:,:,:) = repmat(reshape(cy,1,ny,1),[nx 1 nz]);
%! coefs(3,:,:,:) = repmat(reshape(cz,1,1,nz),[nx ny 1]);
%! coefs(4,:,:,:) = 1;
%! nurbs = nrbmak(coefs, knots);
%! x = rand(5,1); y = rand(5,1); z = rand(5,1);
%! tt = [x y z]';
%! points = nrbeval(nurbs,tt);
%! assert(points,[x.^2 y.^3 z]',1e-10);
%!
%!test
%! knots{1} = [0 0 0 1 1 1];
%! knots{2} = [0 0 0 0 1 1 1 1];
%! knots{3} = [0 0 1 1];
%! cx = [0 0 1]; nx = length(cx);
%! cy = [0 0 0 1]; ny = length(cy);
%! cz = [0 1]; nz = length(cz);
%! coefs(1,:,:,:) = repmat(reshape(cx,nx,1,1),[1 ny nz]);
%! coefs(2,:,:,:) = repmat(reshape(cy,1,ny,1),[nx 1 nz]);
%! coefs(3,:,:,:) = repmat(reshape(cz,1,1,nz),[nx ny 1]);
%! coefs(4,:,:,:) = 1;
%! coefs = coefs([2 1 3 4],:,:,:);
%! nurbs = nrbmak(coefs, knots);
%! x = rand(5,1); y = rand(5,1); z = rand(5,1);
%! tt = [x y z]';
%! points = nrbeval(nurbs,tt);
%! [y.^3 x.^2 z]';
%! assert(points,[y.^3 x.^2 z]',1e-10);