|
8215
|
1 function varargout = nrbdeval (nurbs, dnurbs, varargin) |
|
5552
|
2 |
|
8215
|
3 % NRBDEVAL: Evaluation of the derivative and second derivatives of NURBS curve, surface or volume. |
|
5552
|
4 % |
|
8215
|
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}) |
|
5552
|
11 % |
|
|
12 % INPUTS: |
|
|
13 % |
|
8215
|
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) |
|
5552
|
19 % tt - parametric evaluation points |
|
6910
|
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: |
|
5552
|
24 % |
|
|
25 % pnt - evaluated points. |
|
|
26 % jac - evaluated first derivatives (Jacobian). |
|
8215
|
27 % hess - evaluated second derivatives (Hessian). |
|
6910
|
28 % |
|
7114
|
29 % Copyright (C) 2000 Mark Spink |
|
|
30 % Copyright (C) 2010 Carlo de Falco |
|
8215
|
31 % Copyright (C) 2010, 2011 Rafael Vazquez |
|
6910
|
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 |
|
11634
|
35 % the Free Software Foundation, either version 3 of the License, or |
|
6910
|
36 % (at your option) any later version. |
|
5552
|
37 |
|
6910
|
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/>. |
|
5552
|
45 |
|
8215
|
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 |
|
8214
|
55 if (~isstruct(nurbs)) |
|
5552
|
56 error('NURBS representation is not structure!'); |
|
|
57 end |
|
|
58 |
|
8214
|
59 if (~strcmp(nurbs.form,'B-NURBS')) |
|
5552
|
60 error('Not a recognised NURBS representation'); |
|
|
61 end |
|
|
62 |
|
8215
|
63 [cp,cw] = nrbeval(nurbs, tt); |
|
5552
|
64 |
|
8214
|
65 if (iscell(nurbs.knots)) |
|
|
66 if (size(nurbs.knots,2) == 3) |
|
6910
|
67 % NURBS structure represents a volume |
|
|
68 temp = cw(ones(3,1),:,:,:); |
|
|
69 pnt = cp./temp; |
|
|
70 |
|
8214
|
71 [cup,cuw] = nrbeval (dnurbs{1}, tt); |
|
6910
|
72 tempu = cuw(ones(3,1),:,:,:); |
|
|
73 jac{1} = (cup-tempu.*pnt)./temp; |
|
5552
|
74 |
|
8214
|
75 [cvp,cvw] = nrbeval (dnurbs{2}, tt); |
|
6910
|
76 tempv = cvw(ones(3,1),:,:,:); |
|
|
77 jac{2} = (cvp-tempv.*pnt)./temp; |
|
|
78 |
|
8214
|
79 [cwp,cww] = nrbeval (dnurbs{3}, tt); |
|
6910
|
80 tempw = cww(ones(3,1),:,:,:); |
|
|
81 jac{3} = (cwp-tempw.*pnt)./temp; |
|
|
82 |
|
8286
|
83 % second derivatives |
|
8215
|
84 if (nargout == 3) |
|
8286
|
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 |
|
8215
|
122 end |
|
|
123 |
|
8214
|
124 elseif (size(nurbs.knots,2) == 2) |
|
6910
|
125 % NURBS structure represents a surface |
|
|
126 temp = cw(ones(3,1),:,:); |
|
|
127 pnt = cp./temp; |
|
5552
|
128 |
|
8214
|
129 [cup,cuw] = nrbeval (dnurbs{1}, tt); |
|
6910
|
130 tempu = cuw(ones(3,1),:,:); |
|
|
131 jac{1} = (cup-tempu.*pnt)./temp; |
|
|
132 |
|
8214
|
133 [cvp,cvw] = nrbeval (dnurbs{2}, tt); |
|
6910
|
134 tempv = cvw(ones(3,1),:,:); |
|
|
135 jac{2} = (cvp-tempv.*pnt)./temp; |
|
5552
|
136 |
|
8215
|
137 % second derivatives |
|
8286
|
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; |
|
8215
|
143 |
|
8286
|
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; |
|
8215
|
147 |
|
8286
|
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 |
|
8215
|
156 end |
|
|
157 |
|
6910
|
158 end |
|
5552
|
159 else |
|
|
160 |
|
|
161 % NURBS is a curve |
|
|
162 temp = cw(ones(3,1),:); |
|
|
163 pnt = cp./temp; |
|
|
164 |
|
|
165 % first derivative |
|
8214
|
166 [cup,cuw] = nrbeval (dnurbs,tt); |
|
5552
|
167 temp1 = cuw(ones(3,1),:); |
|
|
168 jac = (cup-temp1.*pnt)./temp; |
|
12567
|
169 if (iscell (tt)) |
|
|
170 jac = {jac}; |
|
|
171 end |
|
5552
|
172 |
|
8215
|
173 % second derivative |
|
|
174 if (nargout == 3 && exist ('dnurbs2')) |
|
|
175 [cuup,cuuw] = nrbeval (dnurbs2, tt); |
|
|
176 temp2 = cuuw(ones(3,1),:); |
|
|
177 hess = (cuup - (2*cup.*temp1 + cp.*temp2)./temp + 2*cp.*temp1.^2./temp.^2)./temp; |
|
12567
|
178 if (iscell (tt)) |
|
|
179 hess = {hess}; |
|
|
180 end |
|
8215
|
181 end |
|
12567
|
182 |
|
8215
|
183 end |
|
|
184 |
|
|
185 varargout{1} = pnt; |
|
|
186 varargout{2} = jac; |
|
|
187 if (nargout == 3) |
|
|
188 varargout{3} = hess; |
|
5552
|
189 end |
|
|
190 |
|
6910
|
191 end |
|
5552
|
192 |
|
6910
|
193 %!demo |
|
|
194 %! crv = nrbtestcrv; |
|
|
195 %! nrbplot(crv,48); |
|
|
196 %! title('First derivatives along a test curve.'); |
|
|
197 %! |
|
|
198 %! tt = linspace(0.0,1.0,9); |
|
|
199 %! |
|
|
200 %! dcrv = nrbderiv(crv); |
|
|
201 %! |
|
|
202 %! [p1, dp] = nrbdeval(crv,dcrv,tt); |
|
|
203 %! |
|
|
204 %! p2 = vecnorm(dp); |
|
|
205 %! |
|
|
206 %! hold on; |
|
|
207 %! plot(p1(1,:),p1(2,:),'ro'); |
|
|
208 %! h = quiver(p1(1,:),p1(2,:),p2(1,:),p2(2,:),0); |
|
|
209 %! set(h,'Color','black'); |
|
|
210 %! hold off; |
|
|
211 |
|
|
212 %!demo |
|
|
213 %! srf = nrbtestsrf; |
|
|
214 %! p = nrbeval(srf,{linspace(0.0,1.0,20) linspace(0.0,1.0,20)}); |
|
|
215 %! h = surf(squeeze(p(1,:,:)),squeeze(p(2,:,:)),squeeze(p(3,:,:))); |
|
|
216 %! set(h,'FaceColor','blue','EdgeColor','blue'); |
|
|
217 %! title('First derivatives over a test surface.'); |
|
|
218 %! |
|
|
219 %! npts = 5; |
|
|
220 %! tt = linspace(0.0,1.0,npts); |
|
|
221 %! dsrf = nrbderiv(srf); |
|
|
222 %! |
|
|
223 %! [p1, dp] = nrbdeval(srf, dsrf, {tt, tt}); |
|
|
224 %! |
|
|
225 %! up2 = vecnorm(dp{1}); |
|
|
226 %! vp2 = vecnorm(dp{2}); |
|
|
227 %! |
|
|
228 %! hold on; |
|
|
229 %! plot3(p1(1,:),p1(2,:),p1(3,:),'ro'); |
|
|
230 %! h1 = quiver3(p1(1,:),p1(2,:),p1(3,:),up2(1,:),up2(2,:),up2(3,:)); |
|
|
231 %! h2 = quiver3(p1(1,:),p1(2,:),p1(3,:),vp2(1,:),vp2(2,:),vp2(3,:)); |
|
|
232 %! set(h1,'Color','black'); |
|
|
233 %! set(h2,'Color','black'); |
|
|
234 %! |
|
|
235 %! hold off; |
|
|
236 |
|
|
237 %!test |
|
|
238 %! knots{1} = [0 0 0 1 1 1]; |
|
|
239 %! knots{2} = [0 0 0 .5 1 1 1]; |
|
|
240 %! knots{3} = [0 0 0 0 1 1 1 1]; |
|
|
241 %! cx = [0 0.5 1]; nx = length(cx); |
|
|
242 %! cy = [0 0.25 0.75 1]; ny = length(cy); |
|
|
243 %! cz = [0 1/3 2/3 1]; nz = length(cz); |
|
|
244 %! coefs(1,:,:,:) = repmat(reshape(cx,nx,1,1),[1 ny nz]); |
|
|
245 %! coefs(2,:,:,:) = repmat(reshape(cy,1,ny,1),[nx 1 nz]); |
|
|
246 %! coefs(3,:,:,:) = repmat(reshape(cz,1,1,nz),[nx ny 1]); |
|
|
247 %! coefs(4,:,:,:) = 1; |
|
|
248 %! nurbs = nrbmak(coefs, knots); |
|
|
249 %! x = rand(5,1); y = rand(5,1); z = rand(5,1); |
|
|
250 %! tt = [x y z]'; |
|
|
251 %! ders = nrbderiv(nurbs); |
|
|
252 %! [points,jac] = nrbdeval(nurbs,ders,tt); |
|
|
253 %! assert(points,tt,1e-10) |
|
|
254 %! assert(jac{1}(1,:,:),ones(size(jac{1}(1,:,:))),1e-12) |
|
|
255 %! assert(jac{2}(2,:,:),ones(size(jac{2}(2,:,:))),1e-12) |
|
|
256 %! assert(jac{3}(3,:,:),ones(size(jac{3}(3,:,:))),1e-12) |
|
|
257 %! |
|
|
258 %!test |
|
|
259 %! knots{1} = [0 0 0 1 1 1]; |
|
|
260 %! knots{2} = [0 0 0 0 1 1 1 1]; |
|
|
261 %! knots{3} = [0 0 0 1 1 1]; |
|
|
262 %! cx = [0 0 1]; nx = length(cx); |
|
|
263 %! cy = [0 0 0 1]; ny = length(cy); |
|
|
264 %! cz = [0 0.5 1]; nz = length(cz); |
|
|
265 %! coefs(1,:,:,:) = repmat(reshape(cx,nx,1,1),[1 ny nz]); |
|
|
266 %! coefs(2,:,:,:) = repmat(reshape(cy,1,ny,1),[nx 1 nz]); |
|
|
267 %! coefs(3,:,:,:) = repmat(reshape(cz,1,1,nz),[nx ny 1]); |
|
|
268 %! coefs(4,:,:,:) = 1; |
|
|
269 %! coefs = coefs([2 1 3 4],:,:,:); |
|
|
270 %! nurbs = nrbmak(coefs, knots); |
|
|
271 %! x = rand(5,1); y = rand(5,1); z = rand(5,1); |
|
|
272 %! tt = [x y z]'; |
|
|
273 %! dnurbs = nrbderiv(nurbs); |
|
|
274 %! [points, jac] = nrbdeval(nurbs,dnurbs,tt); |
|
|
275 %! assert(points,[y.^3 x.^2 z]',1e-10); |
|
|
276 %! assert(jac{2}(1,:,:),3*y'.^2,1e-12) |
|
|
277 %! assert(jac{1}(2,:,:),2*x',1e-12) |
|
|
278 %! assert(jac{3}(3,:,:),ones(size(z')),1e-12) |
