octave_packages/nurbs-1.3.6/nrbdeval.m

   1 function varargout = nrbdeval (nurbs, dnurbs, varargin)
   2
   3 % NRBDEVAL: Evaluation of the derivative and second derivatives of NURBS curve, surface or volume.
   4 %
   5 %     [pnt, jac] = nrbdeval (crv, dcrv, tt)
   6 %     [pnt, jac] = nrbdeval (srf, dsrf, {tu tv})
   7 %     [pnt, jac] = nrbdeval (vol, dvol, {tu tv tw})
   8 %     [pnt, jac, hess] = nrbdeval (crv, dcrv, dcrv2, tt)
   9 %     [pnt, jac, hess] = nrbdeval (srf, dsrf, dsrf2, {tu tv})
  10 %     [pnt, jac, hess] = nrbdeval (vol, dvol, {tu tv tw})
  11 %
  12 % INPUTS:
  13 %
  14 %   crv,   srf,   vol   - original NURBS curve, surface or volume.
  15 %   dcrv,  dsrf,  dvol  - NURBS derivative representation of crv, srf
  16 %                          or vol (see nrbderiv2)
  17 %   dcrv2, dsrf2, dvol2 - NURBS second derivative representation of crv,
  18 %                          srf or vol (see nrbderiv2)
  19 %   tt     - parametric evaluation points
  20 %            If the nurbs is a surface or a volume then tt is a cell
  21 %            {tu, tv} or {tu, tv, tw} are the parametric coordinates
  22 %
  23 % OUTPUT:
  24 %
  25 %   pnt  - evaluated points.
  26 %   jac  - evaluated first derivatives (Jacobian).
  27 %   hess - evaluated second derivatives (Hessian).
  28 %
  29 % Copyright (C) 2000 Mark Spink
  30 % Copyright (C) 2010 Carlo de Falco
  31 % Copyright (C) 2010, 2011 Rafael Vazquez
  32 %
  33 %    This program is free software: you can redistribute it and/or modify
  34 %    it under the terms of the GNU General Public License as published by
  35 %    the Free Software Foundation, either version 2 of the License, or
  36 %    (at your option) any later version.
  37
  38 %    This program is distributed in the hope that it will be useful,
  39 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  40 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  41 %    GNU General Public License for more details.
  42 %
  43 %    You should have received a copy of the GNU General Public License
  44 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  45
  46 if (nargin == 3)
  47   tt = varargin{1};
  48 elseif (nargin == 4)
  49   dnurbs2 = varargin{1};
  50   tt = varargin{2};
  51 else
  52   error ('nrbrdeval: wrong number of input parameters')
  53 end
  54
  55 if (~isstruct(nurbs))
  56   error('NURBS representation is not structure!');
  57 end
  58
  59 if (~strcmp(nurbs.form,'B-NURBS'))
  60   error('Not a recognised NURBS representation');
  61 end
  62
  63 [cp,cw] = nrbeval(nurbs, tt);
  64
  65 if (iscell(nurbs.knots))
  66   if (size(nurbs.knots,2) == 3)
  67   % NURBS structure represents a volume
  68     temp = cw(ones(3,1),:,:,:);
  69     pnt = cp./temp;
  70
  71     [cup,cuw] = nrbeval (dnurbs{1}, tt);
  72     tempu = cuw(ones(3,1),:,:,:);
  73     jac{1} = (cup-tempu.*pnt)./temp;
  74
  75     [cvp,cvw] = nrbeval (dnurbs{2}, tt);
  76     tempv = cvw(ones(3,1),:,:,:);
  77     jac{2} = (cvp-tempv.*pnt)./temp;
  78
  79     [cwp,cww] = nrbeval (dnurbs{3}, tt);
  80     tempw = cww(ones(3,1),:,:,:);
  81     jac{3} = (cwp-tempw.*pnt)./temp;
  82
  83 % second derivatives
  84     if (nargout == 3)
  85       if (exist ('dnurbs2'))
  86         [cuup, cuuw] = nrbeval (dnurbs2{1,1}, tt);
  87         tempuu = cuuw(ones(3,1),:,:,:);
  88         hess{1,1} = (cuup - (2*cup.*tempu + cp.*tempuu)./temp + 2*cp.*tempu.^2./temp.^2)./temp;
  89         clear cuup cuuw tempuu
  90
  91         [cvvp, cvvw] = nrbeval (dnurbs2{2,2}, tt);
  92         tempvv = cvvw(ones(3,1),:,:,:);
  93         hess{2,2} = (cvvp - (2*cvp.*tempv + cp.*tempvv)./temp + 2*cp.*tempv.^2./temp.^2)./temp;
  94         clear cvvp cvvw tempvv
  95
  96         [cwwp, cwww] = nrbeval (dnurbs2{3,3}, tt);
  97         tempww = cwww(ones(3,1),:,:,:);
  98         hess{3,3} = (cwwp - (2*cwp.*tempw + cp.*tempww)./temp + 2*cp.*tempw.^2./temp.^2)./temp;
  99         clear cwwp cwww tempww
 100
 101         [cuvp, cuvw] = nrbeval (dnurbs2{1,2}, tt);
 102         tempuv = cuvw(ones(3,1),:,:,:);
 103         hess{1,2} = (cuvp - (cup.*tempv + cvp.*tempu + cp.*tempuv)./temp + 2*cp.*tempu.*tempv./temp.^2)./temp;
 104         hess{2,1} = hess{1,2};
 105         clear cuvp cuvw tempuv
 106
 107         [cuwp, cuww] = nrbeval (dnurbs2{1,3}, tt);
 108         tempuw = cuww(ones(3,1),:,:,:);
 109         hess{1,3} = (cuwp - (cup.*tempw + cwp.*tempu + cp.*tempuw)./temp + 2*cp.*tempu.*tempw./temp.^2)./temp;
 110         hess{3,1} = hess{1,3};
 111         clear cuwp cuww tempuw
 112
 113         [cvwp, cvww] = nrbeval (dnurbs2{2,3}, tt);
 114         tempvw = cvww(ones(3,1),:,:,:);
 115         hess{2,3} = (cvwp - (cvp.*tempw + cwp.*tempv + cp.*tempvw)./temp + 2*cp.*tempv.*tempw./temp.^2)./temp;
 116         hess{3,2} = hess{2,3};
 117         clear cvwp cvww tempvw
 118       else
 119         warning ('nrbdeval: dnurbs2 missing. The second derivative is not computed');
 120         hess = [];
 121       end
 122     end
 123
 124   elseif (size(nurbs.knots,2) == 2)
 125   % NURBS structure represents a surface
 126     temp = cw(ones(3,1),:,:);
 127     pnt = cp./temp;
 128
 129     [cup,cuw] = nrbeval (dnurbs{1}, tt);
 130     tempu = cuw(ones(3,1),:,:);
 131     jac{1} = (cup-tempu.*pnt)./temp;
 132
 133     [cvp,cvw] = nrbeval (dnurbs{2}, tt);
 134     tempv = cvw(ones(3,1),:,:);
 135     jac{2} = (cvp-tempv.*pnt)./temp;
 136
 137 % second derivatives
 138     if (nargout == 3)
 139       if (exist ('dnurbs2'))
 140         [cuup, cuuw] = nrbeval (dnurbs2{1,1}, tt);
 141         tempuu = cuuw(ones(3,1),:,:);
 142         hess{1,1} = (cuup - (2*cup.*tempu + cp.*tempuu)./temp + 2*cp.*tempu.^2./temp.^2)./temp;
 143
 144         [cvvp, cvvw] = nrbeval (dnurbs2{2,2}, tt);
 145         tempvv = cvvw(ones(3,1),:,:);
 146         hess{2,2} = (cvvp - (2*cvp.*tempv + cp.*tempvv)./temp + 2*cp.*tempv.^2./temp.^2)./temp;
 147
 148         [cuvp, cuvw] = nrbeval (dnurbs2{1,2}, tt);
 149         tempuv = cuvw(ones(3,1),:,:);
 150         hess{1,2} = (cuvp - (cup.*tempv + cvp.*tempu + cp.*tempuv)./temp + 2*cp.*tempu.*tempv./temp.^2)./temp;
 151         hess{2,1} = hess{1,2};
 152       else
 153         warning ('nrbdeval: dnurbs2 missing. The second derivative is not computed');
 154         hess = [];
 155       end
 156     end
 157
 158   end
 159 else
 160
 161   % NURBS is a curve
 162   temp = cw(ones(3,1),:);
 163   pnt = cp./temp;
 164
 165   % first derivative
 166   [cup,cuw] = nrbeval (dnurbs,tt);
 167   temp1 = cuw(ones(3,1),:);
 168   jac = (cup-temp1.*pnt)./temp;
 169
 170   % second derivative
 171   if (nargout == 3 && exist ('dnurbs2'))
 172     [cuup,cuuw] = nrbeval (dnurbs2, tt);
 173     temp2 = cuuw(ones(3,1),:);
 174     hess = (cuup - (2*cup.*temp1 + cp.*temp2)./temp + 2*cp.*temp1.^2./temp.^2)./temp;
 175   end
 176
 177 end
 178
 179 varargout{1} = pnt;
 180 varargout{2} = jac;
 181 if (nargout == 3)
 182   varargout{3} = hess;
 183 end
 184
 185 end
 186
 187 %!demo
 188 %! crv = nrbtestcrv;
 189 %! nrbplot(crv,48);
 190 %! title('First derivatives along a test curve.');
 191 %!
 192 %! tt = linspace(0.0,1.0,9);
 193 %!
 194 %! dcrv = nrbderiv(crv);
 195 %!
 196 %! [p1, dp] = nrbdeval(crv,dcrv,tt);
 197 %!
 198 %! p2 = vecnorm(dp);
 199 %!
 200 %! hold on;
 201 %! plot(p1(1,:),p1(2,:),'ro');
 202 %! h = quiver(p1(1,:),p1(2,:),p2(1,:),p2(2,:),0);
 203 %! set(h,'Color','black');
 204 %! hold off;
 205
 206 %!demo
 207 %! srf = nrbtestsrf;
 208 %! p = nrbeval(srf,{linspace(0.0,1.0,20) linspace(0.0,1.0,20)});
 209 %! h = surf(squeeze(p(1,:,:)),squeeze(p(2,:,:)),squeeze(p(3,:,:)));
 210 %! set(h,'FaceColor','blue','EdgeColor','blue');
 211 %! title('First derivatives over a test surface.');
 212 %!
 213 %! npts = 5;
 214 %! tt = linspace(0.0,1.0,npts);
 215 %! dsrf = nrbderiv(srf);
 216 %!
 217 %! [p1, dp] = nrbdeval(srf, dsrf, {tt, tt});
 218 %!
 219 %! up2 = vecnorm(dp{1});
 220 %! vp2 = vecnorm(dp{2});
 221 %!
 222 %! hold on;
 223 %! plot3(p1(1,:),p1(2,:),p1(3,:),'ro');
 224 %! h1 = quiver3(p1(1,:),p1(2,:),p1(3,:),up2(1,:),up2(2,:),up2(3,:));
 225 %! h2 = quiver3(p1(1,:),p1(2,:),p1(3,:),vp2(1,:),vp2(2,:),vp2(3,:));
 226 %! set(h1,'Color','black');
 227 %! set(h2,'Color','black');
 228 %!
 229 %! hold off;
 230
 231 %!test
 232 %! knots{1} = [0 0 0 1 1 1];
 233 %! knots{2} = [0 0 0 .5 1 1 1];
 234 %! knots{3} = [0 0 0 0 1 1 1 1];
 235 %! cx = [0 0.5 1]; nx = length(cx);
 236 %! cy = [0 0.25 0.75 1]; ny = length(cy);
 237 %! cz = [0 1/3 2/3 1]; nz = length(cz);
 238 %! coefs(1,:,:,:) = repmat(reshape(cx,nx,1,1),[1 ny nz]);
 239 %! coefs(2,:,:,:) = repmat(reshape(cy,1,ny,1),[nx 1 nz]);
 240 %! coefs(3,:,:,:) = repmat(reshape(cz,1,1,nz),[nx ny 1]);
 241 %! coefs(4,:,:,:) = 1;
 242 %! nurbs = nrbmak(coefs, knots);
 243 %! x = rand(5,1); y = rand(5,1); z = rand(5,1);
 244 %! tt = [x y z]';
 245 %! ders = nrbderiv(nurbs);
 246 %! [points,jac] = nrbdeval(nurbs,ders,tt);
 247 %! assert(points,tt,1e-10)
 248 %! assert(jac{1}(1,:,:),ones(size(jac{1}(1,:,:))),1e-12)
 249 %! assert(jac{2}(2,:,:),ones(size(jac{2}(2,:,:))),1e-12)
 250 %! assert(jac{3}(3,:,:),ones(size(jac{3}(3,:,:))),1e-12)
 251 %!
 252 %!test
 253 %! knots{1} = [0 0 0 1 1 1];
 254 %! knots{2} = [0 0 0 0 1 1 1 1];
 255 %! knots{3} = [0 0 0 1 1 1];
 256 %! cx = [0 0 1]; nx = length(cx);
 257 %! cy = [0 0 0 1]; ny = length(cy);
 258 %! cz = [0 0.5 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 %! coefs = coefs([2 1 3 4],:,:,:);
 264 %! nurbs = nrbmak(coefs, knots);
 265 %! x = rand(5,1); y = rand(5,1); z = rand(5,1);
 266 %! tt = [x y z]';
 267 %! dnurbs = nrbderiv(nurbs);
 268 %! [points, jac] = nrbdeval(nurbs,dnurbs,tt);
 269 %! assert(points,[y.^3 x.^2 z]',1e-10);
 270 %! assert(jac{2}(1,:,:),3*y'.^2,1e-12)
 271 %! assert(jac{1}(2,:,:),2*x',1e-12)
 272 %! assert(jac{3}(3,:,:),ones(size(z')),1e-12)