1function 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, dvol2, {tu tv tw}) 11% [pnt, jac, hess, third_der] = nrbdeval (crv, dcrv, dcrv2, dcrv3, tt) 12% [pnt, jac, hess, third_der] = nrbdeval (srf, dsrf, dsrf2, dsrf3, {tu tv}) 13% [pnt, jac, hess, third_der, fourth_der] = nrbdeval (crv, dcrv, dcrv2, dcrv3, dcrv4, tt) 14% [pnt, jac, hess, third_der, fourth_der] = nrbdeval (srf, dsrf, dsrf2, dsrf3, dsrf4, {tu tv}) 15% 16% INPUTS: 17% 18% crv, srf, vol - original NURBS curve, surface or volume. 19% dcrv, dsrf, dvol - NURBS derivative representation of crv, srf 20% or vol (see nrbderiv2) 21% dcrv2, dsrf2, dvol2 - NURBS second derivative representation of crv, 22% srf or vol (see nrbderiv2) 23% dcrv3, dsrf3, dvol3 - NURBS third derivative representation of crv, 24% srf or vol (see nrbderiv) 25% dcrv4, dsrf4, dvol4 - NURBS fourth derivative representation of crv, 26% srf or vol (see nrbderiv) 27% 28% tt - parametric evaluation points 29% If the nurbs is a surface or a volume then tt is a cell 30% {tu, tv} or {tu, tv, tw} are the parametric coordinates 31% 32% OUTPUT: 33% 34% pnt - evaluated points. 35% jac - evaluated first derivatives (Jacobian). 36% hess - evaluated second derivatives (Hessian). 37% third_der - evaluated third derivatives 38% fourth_der - evaluated fourth derivatives 39% 40% Copyright (C) 2000 Mark Spink 41% Copyright (C) 2010 Carlo de Falco 42% Copyright (C) 2010, 2011 Rafael Vazquez 43% Copyright (C) 2018 Luca Coradello 44% 45% This program is free software: you can redistribute it and/or modify 46% it under the terms of the GNU General Public License as published by 47% the Free Software Foundation, either version 3 of the License, or 48% (at your option) any later version. 49 50% This program is distributed in the hope that it will be useful, 51% but WITHOUT ANY WARRANTY; without even the implied warranty of 52% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 53% GNU General Public License for more details. 54% 55% You should have received a copy of the GNU General Public License 56% along with this program. If not, see <http://www.gnu.org/licenses/>. 57 58if (nargin == 3) 59 tt = varargin{1}; 60elseif (nargin == 4) 61 dnurbs2 = varargin{1}; 62 tt = varargin{2}; 63elseif (nargin == 5) 64 dnurbs2 = varargin{1}; 65 dnurbs3 = varargin{2}; 66 tt = varargin{3}; 67elseif (nargin == 6) 68 dnurbs2 = varargin{1}; 69 dnurbs3 = varargin{2}; 70 dnurbs4 = varargin{3}; 71 tt = varargin{4}; 72else 73 error ('nrbrdeval: wrong number of input parameters') 74end 75 76if (~isstruct(nurbs)) 77 error('NURBS representation is not structure!'); 78end 79 80if (~strcmp(nurbs.form,'B-NURBS')) 81 error('Not a recognised NURBS representation'); 82end 83 84[cp,cw] = nrbeval(nurbs, tt); 85 86if (iscell(nurbs.knots)) 87 if (size(nurbs.knots,2) == 3) 88 % NURBS structure represents a volume 89 temp = cw(ones(3,1),:,:,:); 90 pnt = cp./temp; 91 92 [cup,cuw] = nrbeval (dnurbs{1}, tt); 93 tempu = cuw(ones(3,1),:,:,:); 94 jac{1} = (cup-tempu.*pnt)./temp; 95 96 [cvp,cvw] = nrbeval (dnurbs{2}, tt); 97 tempv = cvw(ones(3,1),:,:,:); 98 jac{2} = (cvp-tempv.*pnt)./temp; 99 100 [cwp,cww] = nrbeval (dnurbs{3}, tt); 101 tempw = cww(ones(3,1),:,:,:); 102 jac{3} = (cwp-tempw.*pnt)./temp; 103 104% second derivatives 105 if (nargout >= 3) 106 if (exist ('dnurbs2')) 107 [cuup, cuuw] = nrbeval (dnurbs2{1,1}, tt); 108 tempuu = cuuw(ones(3,1),:,:,:); 109 hess{1,1} = (cuup - (2*cup.*tempu + cp.*tempuu)./temp + 2*cp.*tempu.^2./temp.^2)./temp; 110 clear cuup cuuw tempuu 111 112 [cvvp, cvvw] = nrbeval (dnurbs2{2,2}, tt); 113 tempvv = cvvw(ones(3,1),:,:,:); 114 hess{2,2} = (cvvp - (2*cvp.*tempv + cp.*tempvv)./temp + 2*cp.*tempv.^2./temp.^2)./temp; 115 clear cvvp cvvw tempvv 116 117 [cwwp, cwww] = nrbeval (dnurbs2{3,3}, tt); 118 tempww = cwww(ones(3,1),:,:,:); 119 hess{3,3} = (cwwp - (2*cwp.*tempw + cp.*tempww)./temp + 2*cp.*tempw.^2./temp.^2)./temp; 120 clear cwwp cwww tempww 121 122 [cuvp, cuvw] = nrbeval (dnurbs2{1,2}, tt); 123 tempuv = cuvw(ones(3,1),:,:,:); 124 hess{1,2} = (cuvp - (cup.*tempv + cvp.*tempu + cp.*tempuv)./temp + 2*cp.*tempu.*tempv./temp.^2)./temp; 125 hess{2,1} = hess{1,2}; 126 clear cuvp cuvw tempuv 127 128 [cuwp, cuww] = nrbeval (dnurbs2{1,3}, tt); 129 tempuw = cuww(ones(3,1),:,:,:); 130 hess{1,3} = (cuwp - (cup.*tempw + cwp.*tempu + cp.*tempuw)./temp + 2*cp.*tempu.*tempw./temp.^2)./temp; 131 hess{3,1} = hess{1,3}; 132 clear cuwp cuww tempuw 133 134 [cvwp, cvww] = nrbeval (dnurbs2{2,3}, tt); 135 tempvw = cvww(ones(3,1),:,:,:); 136 hess{2,3} = (cvwp - (cvp.*tempw + cwp.*tempv + cp.*tempvw)./temp + 2*cp.*tempv.*tempw./temp.^2)./temp; 137 hess{3,2} = hess{2,3}; 138 clear cvwp cvww tempvw 139 else 140 warning ('nrbdeval: dnurbs2 missing. The second derivative is not computed'); 141 hess = []; 142 end 143 if (nargout > 3) 144 warning ('nrbdeval: 3rd and 4th order derivatives not implemented for volumes'); 145 third_der = []; 146 fourth_der = []; 147 end 148 end 149 150 elseif (size(nurbs.knots,2) == 2) 151 % NURBS structure represents a surface 152 temp = cw(ones(3,1),:,:); 153 pnt = cp./temp; 154 155 [cup,cuw] = nrbeval (dnurbs{1}, tt); 156 tempu = cuw(ones(3,1),:,:); 157 jac{1} = (cup-tempu.*pnt)./temp; 158 159 [cvp,cvw] = nrbeval (dnurbs{2}, tt); 160 tempv = cvw(ones(3,1),:,:); 161 jac{2} = (cvp-tempv.*pnt)./temp; 162 163% second derivatives 164 if (nargout >= 3) 165 if (exist ('dnurbs2')) 166 [cuup, cuuw] = nrbeval (dnurbs2{1,1}, tt); 167 tempuu = cuuw(ones(3,1),:,:); 168 hess{1,1} = (cuup - (2*cup.*tempu + cp.*tempuu)./temp + 2*cp.*tempu.^2./temp.^2)./temp; 169 170 [cvvp, cvvw] = nrbeval (dnurbs2{2,2}, tt); 171 tempvv = cvvw(ones(3,1),:,:); 172 hess{2,2} = (cvvp - (2*cvp.*tempv + cp.*tempvv)./temp + 2*cp.*tempv.^2./temp.^2)./temp; 173 174 [cuvp, cuvw] = nrbeval (dnurbs2{1,2}, tt); 175 tempuv = cuvw(ones(3,1),:,:); 176 hess{1,2} = (cuvp - (cup.*tempv + cvp.*tempu + cp.*tempuv)./temp + 2*cp.*tempu.*tempv./temp.^2)./temp; 177 hess{2,1} = hess{1,2}; 178 else 179 warning ('nrbdeval: dnurbs2 missing. The second derivative is not computed'); 180 hess = []; 181 end 182 end 183 if (nargout >= 4) 184 if (exist ('dnurbs3')) 185 [cuuup, cuuuw] = nrbeval (dnurbs3{1,1,1}, tt); 186 tempuuu = cuuuw(ones(3,1),:,:); 187 third_der{1,1,1} = (cuuup - 3*jac{1}.*tempuu - cp.*tempuuu - 3*hess{1,1}.*tempu)./temp; 188 189 [cvvvp, cvvvw] = nrbeval (dnurbs3{2,2,2}, tt); 190 tempvvv = cvvvw(ones(3,1),:,:); 191 third_der{2,2,2} = (cvvvp - 3*jac{2}.*tempvv - cp.*tempvvv - 3*hess{2,2}.*tempv)./temp; 192 193 [cuuvp, cuuvw] = nrbeval (dnurbs3{1,1,2}, tt); 194 tempuuv = cuuvw(ones(3,1),:,:); 195 third_der{1,1,2} = (cuuvp - 2*jac{1}.*tempuv - jac{2}.*tempuu - cp.*tempuuv - 2*hess{1,2}.*tempu - hess{1,1}.*tempv)./temp; 196 third_der{1,2,1} = third_der{1,1,2}; 197 third_der{2,1,1} = third_der{1,1,2}; 198 199 [cuvvp, cuvvw] = nrbeval (dnurbs3{1,2,2}, tt); 200 tempuvv = cuvvw(ones(3,1),:,:); 201 third_der{1,2,2} = (cuvvp - 2*jac{2}.*tempuv - jac{1}.*tempvv - cp.*tempuvv - 2*hess{1,2}.*tempv - hess{2,2}.*tempu)./temp; 202 third_der{2,2,1} = third_der{1,2,2}; 203 third_der{2,1,2} = third_der{1,2,2}; 204 205 else 206 warning ('nrbdeval: dnurbs3 missing. The third derivative is not computed'); 207 third_der = []; 208 end 209 end 210 211 if (nargout == 5) 212 if (exist ('dnurbs4')) 213 [cuuuup, cuuuuw] = nrbeval (dnurbs4{1,1,1,1}, tt); 214 tempuuuu = cuuuuw(ones(3,1),:,:); 215 fourth_der{1,1,1,1} = (cuuuup - cp.*tempuuuu - 4*jac{1}.*tempuuu -6*hess{1,1}.*tempuu - 4*third_der{1,1,1}.*tempu)./temp; 216 217 [cvvvvp, cvvvvw] = nrbeval (dnurbs4{2,2,2,2}, tt); 218 tempvvvv = cvvvvw(ones(3,1),:,:); 219 fourth_der{2,2,2,2} = (cvvvvp - cp.*tempvvvv - 4*jac{2}.*tempvvv -6*hess{2,2}.*tempvv - 4*third_der{2,2,2}.*tempv)./temp; 220 221 [cuuuvp, cuuuvw] = nrbeval (dnurbs4{1,1,1,2}, tt); 222 tempuuuv = cuuuvw(ones(3,1),:,:); 223 fourth_der{1,1,1,2} = (cuuuvp - cp.*tempuuuv - 3*jac{1}.*tempuuv - jac{2}.*tempuuu -3*hess{1,2}.*tempuu -3*hess{1,1}.*tempuv ... 224 - third_der{1,1,1}.*tempv - 3*third_der{1,1,2}.*tempu)./temp; 225 226 fourth_der{1,1,2,1} = fourth_der{1,1,1,2}; 227 fourth_der{1,2,1,1} = fourth_der{1,1,1,2}; 228 fourth_der{2,1,1,1} = fourth_der{1,1,1,2}; 229 230 [cuvvvp, cuvvvw] = nrbeval (dnurbs4{1,2,2,2}, tt); 231 tempuvvv = cuvvvw(ones(3,1),:,:); 232 fourth_der{1,2,2,2} = (cuvvvp - cp.*tempuvvv - 3*jac{2}.*tempuvv - jac{1}.*tempvvv -3*hess{1,2}.*tempvv -3*hess{2,2}.*tempuv ... 233 - third_der{2,2,2}.*tempu - 3*third_der{1,2,2}.*tempv)./temp; 234 235 fourth_der{2,2,1,2} = fourth_der{1,2,2,2}; 236 fourth_der{2,1,2,2} = fourth_der{1,2,2,2}; 237 fourth_der{2,2,2,1} = fourth_der{1,2,2,2}; 238 239 [cuuvvp, cuuvvw] = nrbeval (dnurbs4{1,1,2,2}, tt); 240 tempuuvv = cuuvvw(ones(3,1),:,:); 241 fourth_der{1,1,2,2} = (cuuvvp - cp.*tempuuvv - 2*jac{1}.*tempuvv - 2*jac{2}.*tempuuv -hess{1,1}.*tempvv -hess{2,2}.*tempuu ... 242 -4*hess{1,2}.*tempuv - 2*third_der{1,1,2}.*tempv - 2*third_der{1,2,2}.*tempu)./temp; 243 244 fourth_der{1,2,2,1} = fourth_der{1,1,2,2}; 245 fourth_der{2,2,1,1} = fourth_der{1,1,2,2}; 246 fourth_der{1,2,1,2} = fourth_der{1,1,2,2}; 247 fourth_der{2,1,2,1} = fourth_der{1,1,2,2}; 248 fourth_der{2,1,1,2} = fourth_der{1,1,2,2}; 249 250 clear cuup cuuw tempuu cvvp cvvw tempvv cuvp cuvw tempuv 251 clear cuuup cuuuw tempuuu cvvvp cvvvw tempvvv cuuvp cuuvw tempuuv cuvvp cuvvw tempuvv 252 clear cuuuup cuuuuw tempuuuu cvvvvp cvvvvw tempvvvv cuuuvp cuuuvw tempuuuv cuuvvp cuuvvw tempuuvv cvvvup cvvvuw tempvvvu 253 254 else 255 warning ('nrbdeval: dnurbs4 missing. The fourth derivative is not computed'); 256 fourth_der = []; 257 end 258 end 259 260 261 end 262else 263 264 % NURBS is a curve 265 temp = cw(ones(3,1),:); 266 pnt = cp./temp; 267 268 % first derivative 269 [cup,cuw] = nrbeval (dnurbs,tt); 270 temp1 = cuw(ones(3,1),:); 271 jac = (cup-temp1.*pnt)./temp; 272 273 % second derivative 274 if (nargout >= 3 && exist ('dnurbs2')) 275 [cuup,cuuw] = nrbeval (dnurbs2, tt); 276 temp2 = cuuw(ones(3,1),:); 277 hess = (cuup - (2*cup.*temp1 + cp.*temp2)./temp + 2*cp.*temp1.^2./temp.^2)./temp; 278 if (nargout >= 4 && exist ('dnurbs3')) 279 280 [cuuup, cuuuw] = nrbeval (dnurbs3, tt); 281 temp3 = cuuuw(ones(3,1),:); 282 third_der = (cuuup - 3*jac.*temp2 - cp.*temp3 - 3*hess.*temp1)./temp; 283 if (nargout >= 5 && exist ('dnurbs4')) 284 285 [cuuuup, cuuuuw] = nrbeval (dnurbs4, tt); 286 temp4 = cuuuuw(ones(3,1),:); 287 fourth_der = (cuuuup - cp.*temp4 - 4*jac.*temp3 -6*hess.*temp2 - 4*third_der.*temp1)./temp; 288 if (iscell (tt)) 289 fourth_der = {fourth_der}; 290 end 291 292 end 293 294 if (iscell (tt)) 295 third_der = {third_der}; 296 end 297 end 298 299 if (iscell (tt)) 300 hess = {hess}; 301 end 302 end 303 304 if (iscell (tt)) 305 jac = {jac}; 306 end 307 308end 309 310varargout{1} = pnt; 311varargout{2} = jac; 312if (nargout >= 3) 313 varargout{3} = hess; 314end 315if (nargout >= 4) 316 varargout{4} = third_der; 317end 318if (nargout == 5) 319 varargout{5} = fourth_der; 320end 321 322end 323 324%!demo 325%! crv = nrbtestcrv; 326%! nrbplot(crv,48); 327%! title('First derivatives along a test curve.'); 328%! 329%! tt = linspace(0.0,1.0,9); 330%! 331%! dcrv = nrbderiv(crv); 332%! 333%! [p1, dp] = nrbdeval(crv,dcrv,tt); 334%! 335%! p2 = vecnormalize(dp); 336%! 337%! hold on; 338%! plot(p1(1,:),p1(2,:),'ro'); 339%! h = quiver(p1(1,:),p1(2,:),p2(1,:),p2(2,:),0); 340%! set(h,'Color','black'); 341%! hold off; 342 343%!demo 344%! srf = nrbtestsrf; 345%! p = nrbeval(srf,{linspace(0.0,1.0,20) linspace(0.0,1.0,20)}); 346%! h = surf(squeeze(p(1,:,:)),squeeze(p(2,:,:)),squeeze(p(3,:,:))); 347%! set(h,'FaceColor','blue','EdgeColor','blue'); 348%! title('First derivatives over a test surface.'); 349%! 350%! npts = 5; 351%! tt = linspace(0.0,1.0,npts); 352%! dsrf = nrbderiv(srf); 353%! 354%! [p1, dp] = nrbdeval(srf, dsrf, {tt, tt}); 355%! 356%! up2 = vecnormalize(dp{1}); 357%! vp2 = vecnormalize(dp{2}); 358%! 359%! hold on; 360%! plot3(p1(1,:),p1(2,:),p1(3,:),'ro'); 361%! h1 = quiver3(p1(1,:),p1(2,:),p1(3,:),up2(1,:),up2(2,:),up2(3,:)); 362%! h2 = quiver3(p1(1,:),p1(2,:),p1(3,:),vp2(1,:),vp2(2,:),vp2(3,:)); 363%! set(h1,'Color','black'); 364%! set(h2,'Color','black'); 365%! 366%! hold off; 367 368%!test 369%! knots{1} = [0 0 0 1 1 1]; 370%! knots{2} = [0 0 0 .5 1 1 1]; 371%! knots{3} = [0 0 0 0 1 1 1 1]; 372%! cx = [0 0.5 1]; nx = length(cx); 373%! cy = [0 0.25 0.75 1]; ny = length(cy); 374%! cz = [0 1/3 2/3 1]; nz = length(cz); 375%! coefs(1,:,:,:) = repmat(reshape(cx,nx,1,1),[1 ny nz]); 376%! coefs(2,:,:,:) = repmat(reshape(cy,1,ny,1),[nx 1 nz]); 377%! coefs(3,:,:,:) = repmat(reshape(cz,1,1,nz),[nx ny 1]); 378%! coefs(4,:,:,:) = 1; 379%! nurbs = nrbmak(coefs, knots); 380%! x = rand(5,1); y = rand(5,1); z = rand(5,1); 381%! tt = [x y z]'; 382%! ders = nrbderiv(nurbs); 383%! [points,jac] = nrbdeval(nurbs,ders,tt); 384%! assert(points,tt,1e-10) 385%! assert(jac{1}(1,:,:),ones(size(jac{1}(1,:,:))),1e-12) 386%! assert(jac{2}(2,:,:),ones(size(jac{2}(2,:,:))),1e-12) 387%! assert(jac{3}(3,:,:),ones(size(jac{3}(3,:,:))),1e-12) 388%! 389%!test 390%! knots{1} = [0 0 0 1 1 1]; 391%! knots{2} = [0 0 0 0 1 1 1 1]; 392%! knots{3} = [0 0 0 1 1 1]; 393%! cx = [0 0 1]; nx = length(cx); 394%! cy = [0 0 0 1]; ny = length(cy); 395%! cz = [0 0.5 1]; nz = length(cz); 396%! coefs(1,:,:,:) = repmat(reshape(cx,nx,1,1),[1 ny nz]); 397%! coefs(2,:,:,:) = repmat(reshape(cy,1,ny,1),[nx 1 nz]); 398%! coefs(3,:,:,:) = repmat(reshape(cz,1,1,nz),[nx ny 1]); 399%! coefs(4,:,:,:) = 1; 400%! coefs = coefs([2 1 3 4],:,:,:); 401%! nurbs = nrbmak(coefs, knots); 402%! x = rand(5,1); y = rand(5,1); z = rand(5,1); 403%! tt = [x y z]'; 404%! dnurbs = nrbderiv(nurbs); 405%! [points, jac] = nrbdeval(nurbs,dnurbs,tt); 406%! assert(points,[y.^3 x.^2 z]',1e-10); 407%! assert(jac{2}(1,:,:),3*y'.^2,1e-12) 408%! assert(jac{1}(2,:,:),2*x',1e-12) 409%! assert(jac{3}(3,:,:),ones(size(z')),1e-12) 410