octave_packages/nurbs-1.3.6/nrbeval.m

   1 function [p,w] = nrbeval(nurbs,tt)
   2 %
   3 % NRBEVAL: Evaluate a NURBS at parametric points.
   4 %
   5 % Calling Sequences:
   6 %
   7 %   [p,w] = nrbeval(crv,ut)
   8 %   [p,w] = nrbeval(srf,{ut,vt})
   9 %   [p,w] = nrbeval(vol,{ut,vt,wt})
  10 %   [p,w] = nrbeval(srf,pts)
  11 %
  12 % INPUT:
  13 %
  14 %   crv         : NURBS curve, see nrbmak.
  15 %
  16 %   srf         : NURBS surface, see nrbmak.
  17 %
  18 %   vol         : NURBS volume, see nrbmak.
  19 %
  20 %   ut          : Parametric evaluation points along U direction.
  21 %
  22 %   vt          : Parametric evaluation points along V direction.
  23 %
  24 %   wt          : Parametric evaluation points along W direction.
  25 %
  26 %   pts     : Array of scattered points in parametric domain
  27 %
  28 % OUTPUT:
  29 %
  30 %   p           : Evaluated points on the NURBS curve, surface or volume as
  31 %               Cartesian coordinates (x,y,z). If w is included on the lhs argument
  32 %               list the points are returned as homogeneous coordinates (wx,wy,wz).
  33 %
  34 %   w           : Weights of the homogeneous coordinates of the evaluated
  35 %               points. Note inclusion of this argument changes the type
  36 %               of coordinates returned in p (see above).
  37 %
  38 % Description:
  39 %
  40 %   Evaluation of NURBS curves, surfaces or volume at parametric points along
  41 %   the U, V and W directions. Either homogeneous coordinates are returned
  42 %   if the weights are requested in the lhs arguments, or as Cartesian coordinates.
  43 %   This function utilises the 'C' interface bspeval.
  44 %
  45 % Examples:
  46 %
  47 %   Evaluate the NURBS circle at twenty points from 0.0 to 1.0
  48 %
  49 %   nrb = nrbcirc;
  50 %   ut = linspace(0.0,1.0,20);
  51 %   p = nrbeval(nrb,ut);
  52 %
  53 % See also:
  54 %
  55 %     bspeval
  56 %
  57 % Copyright (C) 2000 Mark Spink
  58 % Copyright (C) 2010 Carlo de Falco
  59 % Copyright (C) 2010, 2011 Rafael Vazquez
  60 %
  61 %    This program is free software: you can redistribute it and/or modify
  62 %    it under the terms of the GNU General Public License as published by
  63 %    the Free Software Foundation, either version 2 of the License, or
  64 %    (at your option) any later version.
  65
  66 %    This program is distributed in the hope that it will be useful,
  67 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  68 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  69 %    GNU General Public License for more details.
  70 %
  71 %    You should have received a copy of the GNU General Public License
  72 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  73
  74 if (nargin < 2)
  75   error('Not enough input arguments');
  76 end
  77
  78 foption = 1;    % output format 3D cartesian coordinates
  79 if (nargout == 2)
  80   foption = 0;  % output format 4D homogenous coordinates
  81 end
  82
  83 if (~isstruct(nurbs))
  84   error('NURBS representation is not structure!');
  85 end
  86
  87 if (~strcmp(nurbs.form,'B-NURBS'))
  88   error('Not a recognised NURBS representation');
  89 end
  90
  91 if (iscell(nurbs.knots))
  92   if (size(nurbs.knots,2) == 3)
  93     %% NURBS structure represents a volume
  94
  95     num1 = nurbs.number(1);
  96     num2 = nurbs.number(2);
  97     num3 = nurbs.number(3);
  98     degree = nurbs.order-1;
  99
 100     if (iscell(tt))
 101       nt1 = numel (tt{1});
 102       nt2 = numel (tt{2});
 103       nt3 = numel (tt{3});
 104
 105       %% evaluate along the w direction
 106       val = reshape (nurbs.coefs, 4*num1*num2, num3);
 107       val = bspeval (degree(3), val, nurbs.knots{3}, tt{3});
 108       val = reshape (val, [4 num1 num2 nt3]);
 109
 110       %% Evaluate along the v direction
 111       val = permute (val, [1 2 4 3]);
 112       val = reshape (val, 4*num1*nt3, num2);
 113       val = bspeval (degree(2), val, nurbs.knots{2}, tt{2});
 114       val = reshape (val, [4 num1 nt3 nt2]);
 115       val = permute (val, [1 2 4 3]);
 116
 117       %% Evaluate along the u direction
 118       val = permute (val, [1 3 4 2]);
 119       val = reshape (val, 4*nt2*nt3, num1);
 120       val = bspeval (degree(1), val, nurbs.knots{1}, tt{1});
 121       val = reshape (val, [4 nt2 nt3 nt1]);
 122       val = permute (val, [1 4 2 3]);
 123       pnts = val;
 124
 125       p = pnts(1:3,:,:,:);
 126       w = pnts(4,:,:,:);
 127       if (foption)
 128         p = p./repmat(w,[3 1 1 1]);
 129       end
 130
 131     else
 132
 133       %% Evaluate at scattered points
 134       %% tt(1,:) represents the u direction
 135       %% tt(2,:) represents the v direction
 136       %% tt(3,:) represents the w direction
 137
 138       %% evaluate along the w direction
 139       nt = size(tt,2);
 140       val = reshape(nurbs.coefs,4*num1*num2,num3);
 141       val = bspeval(degree(3),val,nurbs.knots{3},tt(3,:));
 142       val = reshape(val,[4 num1 num2 nt]);
 143
 144       %% evaluate along the v direction
 145       val2 = zeros(4*num1,nt);
 146       for v = 1:nt
 147         coefs = reshape(val(:,:,:,v),4*num1,num2);
 148         val2(:,v) = bspeval(degree(2),coefs,nurbs.knots{2},tt(2,v));
 149       end
 150       val2 = reshape(val2,[4 num1 nt]);
 151
 152       %% evaluate along the u direction
 153       pnts = zeros(4,nt);
 154       for v = 1:nt
 155         coefs = reshape (val2(:,:,v), [4 num1]);
 156         pnts(:,v) = bspeval(degree(1),coefs,nurbs.knots{1},tt(1,v));
 157       end
 158
 159       w = pnts(4,:);
 160       p = pnts(1:3,:);
 161       if (foption)
 162         p = p./repmat(w,[3, 1]);
 163       end
 164     end
 165
 166   elseif (size(nurbs.knots,2) == 2)
 167     %% NURBS structure represents a surface
 168
 169     num1 = nurbs.number(1);
 170     num2 = nurbs.number(2);
 171     degree = nurbs.order-1;
 172
 173     if (iscell(tt))
 174       %% Evaluate over a [u,v] grid
 175       %% tt{1} represents the u direction
 176       %% tt{2} represents the v direction
 177
 178       nt1 = length(tt{1});
 179       nt2 = length(tt{2});
 180
 181       %% Evaluate along the v direction
 182       val = reshape(nurbs.coefs,4*num1,num2);
 183       val = bspeval(degree(2),val,nurbs.knots{2},tt{2});
 184       val = reshape(val,[4 num1 nt2]);
 185
 186       %% Evaluate along the u direction
 187       val = permute(val,[1 3 2]);
 188       val = reshape(val,4*nt2,num1);
 189       val = bspeval(degree(1),val,nurbs.knots{1},tt{1});
 190       val = reshape(val,[4 nt2 nt1]);
 191       val = permute(val,[1 3 2]);
 192
 193       w = val(4,:,:);
 194       p = val(1:3,:,:);
 195       if (foption)
 196         p = p./repmat(w,[3 1 1]);
 197       end
 198
 199     else
 200
 201       %% Evaluate at scattered points
 202       %% tt(1,:) represents the u direction
 203       %% tt(2,:) represents the v direction
 204
 205       nt = size(tt,2);
 206
 207       val = reshape(nurbs.coefs,4*num1,num2);
 208       val = bspeval(degree(2),val,nurbs.knots{2},tt(2,:));
 209       val = reshape(val,[4 num1 nt]);
 210
 211
 212       %% evaluate along the u direction
 213       pnts = zeros(4,nt);
 214       for v = 1:nt
 215         coefs = reshape (val(:,:,v), [4 num1]);
 216         pnts(:,v) = bspeval(degree(1),coefs,nurbs.knots{1},tt(1,v));
 217       end
 218
 219       w = pnts(4,:);
 220       p = pnts(1:3,:);
 221       if (foption)
 222         p = p./repmat(w,[3, 1]);
 223       end
 224
 225     end
 226
 227   end
 228 else
 229
 230   %% NURBS structure represents a curve
 231   %%  tt represent a vector of parametric points in the u direction
 232
 233   val = bspeval(nurbs.order-1,nurbs.coefs,nurbs.knots,tt);
 234
 235   w = val(4,:);
 236   p = val(1:3,:);
 237   if foption
 238     p = p./repmat(w,3,1);
 239   end
 240
 241 end
 242
 243 end
 244
 245 %!demo
 246 %! srf = nrbtestsrf;
 247 %! p = nrbeval(srf,{linspace(0.0,1.0,20) linspace(0.0,1.0,20)});
 248 %! h = surf(squeeze(p(1,:,:)),squeeze(p(2,:,:)),squeeze(p(3,:,:)));
 249 %! title('Test surface.');
 250 %! hold off
 251
 252 %!test
 253 %! knots{1} = [0 0 0 1 1 1];
 254 %! knots{2} = [0 0 0 .5 1 1 1];
 255 %! knots{3} = [0 0 0 0 1 1 1 1];
 256 %! cx = [0 0.5 1]; nx = length(cx);
 257 %! cy = [0 0.25 0.75 1]; ny = length(cy);
 258 %! cz = [0 1/3 2/3 1]; nz = length(cz);
 259 %! coefs(1,:,:,:) = repmat(reshape(cx,nx,1,1),[1 ny nz]);
 260 %! coefs(2,:,:,:) = repmat(reshape(cy,1,ny,1),[nx 1 nz]);
 261 %! coefs(3,:,:,:) = repmat(reshape(cz,1,1,nz),[nx ny 1]);
 262 %! coefs(4,:,:,:) = 1;
 263 %! nurbs = nrbmak(coefs, knots);
 264 %! x = rand(5,1); y = rand(5,1); z = rand(5,1);
 265 %! tt = [x y z]';
 266 %! points = nrbeval(nurbs,tt);
 267 %!
 268 %! assert(points,tt,1e-10)
 269 %!
 270 %!test
 271 %! knots{1} = [0 0 0 1 1 1];
 272 %! knots{2} = [0 0 0 0 1 1 1 1];
 273 %! knots{3} = [0 0 1 1];
 274 %! cx = [0 0 1]; nx = length(cx);
 275 %! cy = [0 0 0 1]; ny = length(cy);
 276 %! cz = [0 1]; nz = length(cz);
 277 %! coefs(1,:,:,:) = repmat(reshape(cx,nx,1,1),[1 ny nz]);
 278 %! coefs(2,:,:,:) = repmat(reshape(cy,1,ny,1),[nx 1 nz]);
 279 %! coefs(3,:,:,:) = repmat(reshape(cz,1,1,nz),[nx ny 1]);
 280 %! coefs(4,:,:,:) = 1;
 281 %! nurbs = nrbmak(coefs, knots);
 282 %! x = rand(5,1); y = rand(5,1); z = rand(5,1);
 283 %! tt = [x y z]';
 284 %! points = nrbeval(nurbs,tt);
 285 %! assert(points,[x.^2 y.^3 z]',1e-10);
 286 %!
 287 %!test
 288 %! knots{1} = [0 0 0 1 1 1];
 289 %! knots{2} = [0 0 0 0 1 1 1 1];
 290 %! knots{3} = [0 0 1 1];
 291 %! cx = [0 0 1]; nx = length(cx);
 292 %! cy = [0 0 0 1]; ny = length(cy);
 293 %! cz = [0 1]; nz = length(cz);
 294 %! coefs(1,:,:,:) = repmat(reshape(cx,nx,1,1),[1 ny nz]);
 295 %! coefs(2,:,:,:) = repmat(reshape(cy,1,ny,1),[nx 1 nz]);
 296 %! coefs(3,:,:,:) = repmat(reshape(cz,1,1,nz),[nx ny 1]);
 297 %! coefs(4,:,:,:) = 1;
 298 %! coefs = coefs([2 1 3 4],:,:,:);
 299 %! nurbs = nrbmak(coefs, knots);
 300 %! x = rand(5,1); y = rand(5,1); z = rand(5,1);
 301 %! tt = [x y z]';
 302 %! points = nrbeval(nurbs,tt);
 303 %! [y.^3 x.^2 z]';
 304 %! assert(points,[y.^3 x.^2 z]',1e-10);