1 function varargout = nrbdeval (nurbs, dnurbs, varargin)
3 % NRBDEVAL: Evaluation of the derivative and second derivatives of NURBS curve, surface or volume.
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})
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
25 % pnt - evaluated points.
26 % jac - evaluated first derivatives (Jacobian).
27 % hess - evaluated second derivatives (Hessian).
29 % Copyright (C) 2000 Mark Spink
30 % Copyright (C) 2010 Carlo de Falco
31 % Copyright (C) 2010, 2011 Rafael Vazquez
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.
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.
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/>.
49 dnurbs2 = varargin{1};
52 error ('nrbrdeval: wrong number of input parameters')
56 error('NURBS representation is not structure!');
59 if (~strcmp(nurbs.form,'B-NURBS'))
60 error('Not a recognised NURBS representation');
63 [cp,cw] = nrbeval(nurbs, tt);
65 if (iscell(nurbs.knots))
66 if (size(nurbs.knots,2) == 3)
67 % NURBS structure represents a volume
68 temp = cw(ones(3,1),:,:,:);
71 [cup,cuw] = nrbeval (dnurbs{1}, tt);
72 tempu = cuw(ones(3,1),:,:,:);
73 jac{1} = (cup-tempu.*pnt)./temp;
75 [cvp,cvw] = nrbeval (dnurbs{2}, tt);
76 tempv = cvw(ones(3,1),:,:,:);
77 jac{2} = (cvp-tempv.*pnt)./temp;
79 [cwp,cww] = nrbeval (dnurbs{3}, tt);
80 tempw = cww(ones(3,1),:,:,:);
81 jac{3} = (cwp-tempw.*pnt)./temp;
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
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
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
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
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
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
119 warning ('nrbdeval: dnurbs2 missing. The second derivative is not computed');
124 elseif (size(nurbs.knots,2) == 2)
125 % NURBS structure represents a surface
126 temp = cw(ones(3,1),:,:);
129 [cup,cuw] = nrbeval (dnurbs{1}, tt);
130 tempu = cuw(ones(3,1),:,:);
131 jac{1} = (cup-tempu.*pnt)./temp;
133 [cvp,cvw] = nrbeval (dnurbs{2}, tt);
134 tempv = cvw(ones(3,1),:,:);
135 jac{2} = (cvp-tempv.*pnt)./temp;
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;
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;
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};
153 warning ('nrbdeval: dnurbs2 missing. The second derivative is not computed');
162 temp = cw(ones(3,1),:);
166 [cup,cuw] = nrbeval (dnurbs,tt);
167 temp1 = cuw(ones(3,1),:);
168 jac = (cup-temp1.*pnt)./temp;
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;
190 %! title('First derivatives along a test curve.');
192 %! tt = linspace(0.0,1.0,9);
194 %! dcrv = nrbderiv(crv);
196 %! [p1, dp] = nrbdeval(crv,dcrv,tt);
201 %! plot(p1(1,:),p1(2,:),'ro');
202 %! h = quiver(p1(1,:),p1(2,:),p2(1,:),p2(2,:),0);
203 %! set(h,'Color','black');
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.');
214 %! tt = linspace(0.0,1.0,npts);
215 %! dsrf = nrbderiv(srf);
217 %! [p1, dp] = nrbdeval(srf, dsrf, {tt, tt});
219 %! up2 = vecnorm(dp{1});
220 %! vp2 = vecnorm(dp{2});
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');
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);
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)
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);
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)