1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // |
---|
27 | // $Id: G4NURBS.cc,v 1.9 2006/06/29 19:06:42 gunter Exp $ |
---|
28 | // GEANT4 tag $Name: geant4-09-03 $ |
---|
29 | // |
---|
30 | // |
---|
31 | // Olivier Crumeyrolle 12 September 1996 |
---|
32 | |
---|
33 | // G4NURBS.cc |
---|
34 | // Implementation of class G4NURBS |
---|
35 | // OC 100796 |
---|
36 | |
---|
37 | |
---|
38 | #include "G4NURBS.hh" |
---|
39 | |
---|
40 | // G4NURBS.hh includes globals.hh which includes a lot of others |
---|
41 | // so no more includes required here |
---|
42 | |
---|
43 | //////////////////////////////////////////////////////////////////////// |
---|
44 | // Here start the real world. Please, check your armored jacket. // |
---|
45 | //////////////////////////////////////////////////////////////////////// |
---|
46 | |
---|
47 | std::ostream & operator << (std::ostream & inout_outStream, |
---|
48 | const G4NURBS & in_kNurb) |
---|
49 | { |
---|
50 | inout_outStream |
---|
51 | // the magic could be changed for good reasons only |
---|
52 | << "##ojc{NURBS}def[1.01.96.7] Just a magic. Could be added to /etc/magic" |
---|
53 | << "\n# NURBS Definition File (human and computer readable format)" |
---|
54 | << "\n# :" << in_kNurb.Whoami() |
---|
55 | << "\n# U order\tV order : " |
---|
56 | << '\n' << in_kNurb.GetUorder() << "\t\t" << in_kNurb.GetVorder(); |
---|
57 | // number of knots and knots themselves for U and V |
---|
58 | for (G4NURBS::t_direction dir = G4NURBS::U; dir < G4NURBS::NofD; |
---|
59 | /*(*(G4int *)(&dir))++*/ dir=(G4NURBS::t_direction)(((G4int)(dir))+1) ) |
---|
60 | { |
---|
61 | inout_outStream |
---|
62 | << "\n# Number of knots along " << G4NURBS::Tochar(dir) |
---|
63 | << '\n' << in_kNurb.GetnbrKnots(dir) |
---|
64 | << "\n# " << G4NURBS::Tochar(dir) << " knots vector (as a column)"; |
---|
65 | { // begin knots iteration |
---|
66 | G4double oneKnot; |
---|
67 | G4NURBS::KnotsIterator knotI(in_kNurb,dir); |
---|
68 | G4bool otherKnots; |
---|
69 | do |
---|
70 | { |
---|
71 | otherKnots = knotI.pick(&oneKnot); |
---|
72 | inout_outStream << "\n\t\t" << oneKnot; |
---|
73 | } |
---|
74 | while (otherKnots); |
---|
75 | } // end of knots iteration |
---|
76 | } // end of direction loop |
---|
77 | |
---|
78 | // number of control points in U and V direction |
---|
79 | // and controlpoints |
---|
80 | inout_outStream |
---|
81 | << "\n# Number of control points along U and V" |
---|
82 | << '\n' << in_kNurb.GetUnbrCtrlPts() |
---|
83 | << " " << in_kNurb.GetVnbrCtrlPts() |
---|
84 | << "\n# Control Points (one by line, U increasing first)"; |
---|
85 | { // begin of control points iteration |
---|
86 | G4NURBS::t_doubleCtrlPt oneCP; |
---|
87 | G4NURBS::CtrlPtsIterator cpI(in_kNurb); |
---|
88 | G4bool otherCPs; |
---|
89 | do |
---|
90 | { |
---|
91 | otherCPs = cpI.pick(&oneCP); |
---|
92 | inout_outStream |
---|
93 | << "\n\t" << oneCP[G4NURBS::X] |
---|
94 | << "\t" << oneCP[G4NURBS::Y] |
---|
95 | << "\t" << oneCP[G4NURBS::Z] |
---|
96 | << "\t" << oneCP[G4NURBS::W]; |
---|
97 | } |
---|
98 | while (otherCPs); |
---|
99 | } // end of control point iteration |
---|
100 | |
---|
101 | inout_outStream << "\n# That's all!" |
---|
102 | << G4endl; // endl do an \n and a flush |
---|
103 | return inout_outStream; |
---|
104 | } |
---|
105 | |
---|
106 | // the CC compiler issue some "maybe no value returned" |
---|
107 | // but everything is ok |
---|
108 | |
---|
109 | G4float G4NURBS::GetfloatKnot(t_direction in_dir, t_indKnot in_index) const |
---|
110 | { |
---|
111 | in_dir = (t_direction)(in_dir & DMask); |
---|
112 | if ( in_index < m[in_dir].nbrKnots ) |
---|
113 | return ((G4float)(m[in_dir].pKnots[in_index])); |
---|
114 | else |
---|
115 | { |
---|
116 | G4cerr << "\nERROR: G4NURBS::GetfloatKnot: index out of range\n" |
---|
117 | << "\n\t in_dir : " << G4int(in_dir) |
---|
118 | << ", in_index : " << G4int(in_index) |
---|
119 | << "m[in_dir].nbrKnots : " << m[in_dir].nbrKnots << G4endl; |
---|
120 | return ((G4float)m[in_dir].pKnots[m[in_dir].nbrKnots-1]); |
---|
121 | } |
---|
122 | } |
---|
123 | |
---|
124 | G4double G4NURBS::GetdoubleKnot(t_direction in_dir, t_indKnot in_index) const |
---|
125 | { |
---|
126 | in_dir = (t_direction)(in_dir & DMask); |
---|
127 | if ( in_index < m[in_dir].nbrKnots ) |
---|
128 | return (G4double)(m[in_dir].pKnots[in_index]); |
---|
129 | else |
---|
130 | { |
---|
131 | G4cerr << "\nERROR: G4NURBS::GetdoubleKnot: index out of range" |
---|
132 | << "\n\t in_dir : " << G4int(in_dir) |
---|
133 | << ", in_index : " << G4int(in_index) |
---|
134 | << "m[in_dir].nbrKnots : " << m[in_dir].nbrKnots |
---|
135 | << G4endl; |
---|
136 | return (G4double)(m[in_dir].pKnots[m[in_dir].nbrKnots-1]); |
---|
137 | } |
---|
138 | } |
---|
139 | |
---|
140 | G4NURBS::t_floatCtrlPt* |
---|
141 | G4NURBS::GetfloatCtrlPt(t_indCtrlPt in_onedimindex) const |
---|
142 | { |
---|
143 | if (in_onedimindex < mtotnbrCtrlPts) |
---|
144 | return TofloatCtrlPt(mpCtrlPts[in_onedimindex]); |
---|
145 | else |
---|
146 | { |
---|
147 | G4cerr << "\nERROR: G4NURBS::GetfloatCtrlPt: index out of range" |
---|
148 | << "\n\t in_onedimindex : " << in_onedimindex |
---|
149 | << " , mtotnbrCtrlPts : " << mtotnbrCtrlPts << G4endl; |
---|
150 | return TofloatCtrlPt(mpCtrlPts[mtotnbrCtrlPts-1]); |
---|
151 | } |
---|
152 | } |
---|
153 | |
---|
154 | G4NURBS::t_floatCtrlPt* |
---|
155 | G4NURBS::GetfloatCtrlPt(t_inddCtrlPt in_Uindex, t_inddCtrlPt in_Vindex) const |
---|
156 | { |
---|
157 | if ( (in_Uindex < m[U].nbrCtrlPts) && (in_Vindex < m[V].nbrCtrlPts) ) |
---|
158 | return TofloatCtrlPt(mpCtrlPts[To1d(in_Uindex, in_Vindex)]); |
---|
159 | else |
---|
160 | { |
---|
161 | G4cerr << "\nERROR: G4NURBS::GetfloatCtrlPt: index(s) out of range" |
---|
162 | << "\n\t in_Uindex : " << in_Uindex |
---|
163 | << " , in_Vindex : " << in_Vindex |
---|
164 | << " , UnbrCtrlPts : " << m[U].nbrCtrlPts |
---|
165 | << " , VnbrCtrlPts : " << m[V].nbrCtrlPts << G4endl; |
---|
166 | return TofloatCtrlPt(mpCtrlPts[mtotnbrCtrlPts-1]); |
---|
167 | } |
---|
168 | } |
---|
169 | |
---|
170 | G4NURBS::t_doubleCtrlPt* |
---|
171 | G4NURBS::GetdoubleCtrlPt(t_indCtrlPt in_onedimindex) const |
---|
172 | { |
---|
173 | if ( in_onedimindex < mtotnbrCtrlPts ) |
---|
174 | return TodoubleCtrlPt(mpCtrlPts[in_onedimindex]); |
---|
175 | else |
---|
176 | { |
---|
177 | G4cerr << "\nERROR: G4NURBS::getdoubleCtrlPts: index out of range" |
---|
178 | << "\n\t in_onedimindex : " << in_onedimindex |
---|
179 | << " , mtotnbrCtrlPts : " << mtotnbrCtrlPts << G4endl; |
---|
180 | return TodoubleCtrlPt(mpCtrlPts[mtotnbrCtrlPts-1]); |
---|
181 | } |
---|
182 | } |
---|
183 | |
---|
184 | G4NURBS::t_doubleCtrlPt* |
---|
185 | G4NURBS::GetdoubleCtrlPt(t_inddCtrlPt in_Uindex, t_inddCtrlPt in_Vindex) const |
---|
186 | { |
---|
187 | if ( (in_Uindex < m[U].nbrCtrlPts) && (in_Vindex < m[V].nbrCtrlPts) ) |
---|
188 | return TodoubleCtrlPt(mpCtrlPts[To1d(in_Uindex, in_Vindex)]); |
---|
189 | else |
---|
190 | { |
---|
191 | G4cerr << "\nERROR: G4NURBS::GetdoubleCtrlPt: index(s) out of range" |
---|
192 | << "\n\t in_Uindex : " << in_Uindex |
---|
193 | << " , in_Vindex : " << in_Vindex |
---|
194 | << " , UnbrCtrlPts : " << m[U].nbrCtrlPts |
---|
195 | << " , VnbrCtrlPts : " << m[V].nbrCtrlPts << G4endl; |
---|
196 | return TodoubleCtrlPt(mpCtrlPts[mtotnbrCtrlPts-1]); |
---|
197 | } |
---|
198 | } |
---|
199 | |
---|
200 | // Total copy |
---|
201 | G4float * G4NURBS::GetfloatAllKnots(t_direction in_dir) const |
---|
202 | { |
---|
203 | in_dir = (t_direction)(in_dir & DMask); |
---|
204 | G4float * p = new G4float [m[in_dir].nbrKnots]; |
---|
205 | for (t_indKnot i = 0; i < m[in_dir].nbrKnots; i++) |
---|
206 | p[i] = (G4float)m[in_dir].pKnots[i]; |
---|
207 | return p; |
---|
208 | } |
---|
209 | |
---|
210 | G4double * G4NURBS::GetdoubleAllKnots(t_direction in_dir) const |
---|
211 | { |
---|
212 | in_dir = (t_direction)(in_dir & DMask); |
---|
213 | G4double * p = new G4double [m[in_dir].nbrKnots]; |
---|
214 | for (t_indKnot i = 0; i < m[in_dir].nbrKnots; i++) |
---|
215 | p[i] = (G4double)m[in_dir].pKnots[i]; |
---|
216 | return p; |
---|
217 | } |
---|
218 | |
---|
219 | G4float * G4NURBS::GetfloatAllCtrlPts() const |
---|
220 | { |
---|
221 | G4float * p = new G4float [mtotnbrCtrlPts*NofC]; |
---|
222 | for (t_indKnot i = 0; i < mtotnbrCtrlPts*NofC; i++) |
---|
223 | p[i] = (G4float)(((t_Coord *)mpCtrlPts)[i]); |
---|
224 | return p; |
---|
225 | } |
---|
226 | |
---|
227 | G4double * G4NURBS::GetdoubleAllCtrlPts() const |
---|
228 | { |
---|
229 | G4double * p = new G4double [mtotnbrCtrlPts*NofC]; |
---|
230 | for (t_indKnot i = 0; i < mtotnbrCtrlPts*NofC; i++) |
---|
231 | p[i] = (G4double)(((t_Coord *)mpCtrlPts)[i]); |
---|
232 | return p; |
---|
233 | } |
---|
234 | |
---|
235 | // Iterators |
---|
236 | |
---|
237 | G4NURBS::KnotsIterator::KnotsIterator(const G4NURBS & in_rNurb, |
---|
238 | G4NURBS::t_direction in_dir, |
---|
239 | t_indKnot in_startIndex) |
---|
240 | : kmdir((G4NURBS::t_direction)(in_dir & G4NURBS::DMask)), |
---|
241 | kmpMax(in_rNurb.m[kmdir].pKnots + in_rNurb.m[kmdir].nbrKnots) |
---|
242 | { |
---|
243 | if (in_startIndex < in_rNurb.m[kmdir].nbrKnots) |
---|
244 | mp = in_rNurb.m[kmdir].pKnots + in_startIndex; |
---|
245 | else |
---|
246 | { |
---|
247 | G4cerr << "\nERROR: G4NURBS::KnotsIterator: in_startIndex out of range" |
---|
248 | << "\n\tin_startIndex : " << in_startIndex |
---|
249 | << ", nbr of knots : " << in_rNurb.m[kmdir].nbrKnots |
---|
250 | << "\n\t mp set to NULL, calls to picking functions will fail" |
---|
251 | << G4endl; |
---|
252 | mp = 0; |
---|
253 | } |
---|
254 | } |
---|
255 | |
---|
256 | G4bool G4NURBS::KnotsIterator::pick(G4double * inout_pDbl) |
---|
257 | { |
---|
258 | (*inout_pDbl) = (G4double)(*mp); |
---|
259 | return (G4bool)((++mp)<kmpMax); |
---|
260 | } |
---|
261 | |
---|
262 | G4bool G4NURBS::KnotsIterator::pick(G4float * inout_pFlt) |
---|
263 | { |
---|
264 | (*inout_pFlt) = (G4float)(*mp); |
---|
265 | return (G4bool)((++mp)<kmpMax); |
---|
266 | } |
---|
267 | |
---|
268 | G4NURBS::CtrlPtsCoordsIterator::CtrlPtsCoordsIterator(const G4NURBS & in_rNurb, |
---|
269 | t_indCtrlPt in_startCtrlPtIndex) |
---|
270 | : kmpMax((const t_Coord *)(in_rNurb.mpCtrlPts + in_rNurb.mtotnbrCtrlPts)) |
---|
271 | { |
---|
272 | if (in_startCtrlPtIndex < in_rNurb.mtotnbrCtrlPts ) |
---|
273 | mp = (const t_Coord *)(in_rNurb.mpCtrlPts + in_startCtrlPtIndex); |
---|
274 | else |
---|
275 | { |
---|
276 | G4cerr << "\nERROR: G4NURBS::CtrlPtsCoordsIterator: " |
---|
277 | << "in_startCtrlPtIndex out of range" |
---|
278 | << "\n\tin_startCtrlPtIndex : " << in_startCtrlPtIndex |
---|
279 | << ", nbr of CtrlPts : " << in_rNurb.mtotnbrCtrlPts |
---|
280 | << "\n\t mp set to NULL, calls to picking functions will fail" |
---|
281 | << G4endl; |
---|
282 | mp = 0; |
---|
283 | } |
---|
284 | } |
---|
285 | |
---|
286 | G4bool G4NURBS::CtrlPtsCoordsIterator::pick(G4double * inout_pDbl) |
---|
287 | { |
---|
288 | (*inout_pDbl) = (G4double)((*mp)); |
---|
289 | return (G4bool)((++mp)<kmpMax); |
---|
290 | } |
---|
291 | |
---|
292 | G4bool G4NURBS::CtrlPtsCoordsIterator::pick(G4float * inout_pFlt) |
---|
293 | { |
---|
294 | (*inout_pFlt) = (G4float)((*mp)); |
---|
295 | return (G4bool)((++mp)<kmpMax); |
---|
296 | } |
---|
297 | |
---|
298 | G4NURBS::CtrlPtsIterator::CtrlPtsIterator(const G4NURBS & in_rNurb, |
---|
299 | t_indCtrlPt in_startIndex) |
---|
300 | : kmpMax(in_rNurb.mpCtrlPts + in_rNurb.mtotnbrCtrlPts) |
---|
301 | { |
---|
302 | if (in_startIndex < in_rNurb.mtotnbrCtrlPts ) |
---|
303 | mp = (in_rNurb.mpCtrlPts + in_startIndex); |
---|
304 | else |
---|
305 | { |
---|
306 | G4cerr << "\nERROR: G4NURBS::CtrlPtsIterator: in_startIndex out of range" |
---|
307 | << "\n\tin_startIndex : " << in_startIndex |
---|
308 | << ", nbr of CtrlPts : " << in_rNurb.mtotnbrCtrlPts |
---|
309 | << "\n\t mp set to NULL, calls to picking functions will fail" |
---|
310 | << G4endl; |
---|
311 | mp = 0; |
---|
312 | } |
---|
313 | } |
---|
314 | |
---|
315 | G4bool G4NURBS::CtrlPtsIterator::pick(t_doubleCtrlPt * inout_pDblCtrlPt) |
---|
316 | { |
---|
317 | for (t_indCoord i = G4NURBS::X; i < G4NURBS::NofC; i++) |
---|
318 | (*inout_pDblCtrlPt)[i] = (G4double)((*mp)[i]); |
---|
319 | return (G4bool)((++mp)<kmpMax); |
---|
320 | } |
---|
321 | |
---|
322 | G4bool G4NURBS::CtrlPtsIterator::pick(t_floatCtrlPt * inout_pFltCtrlPt) |
---|
323 | { |
---|
324 | for (t_indCoord i = G4NURBS::X; i < G4NURBS::NofC; i++) |
---|
325 | (*inout_pFltCtrlPt)[i] = (G4float)((*mp)[i]); |
---|
326 | return (G4bool)((++mp)<kmpMax); |
---|
327 | } |
---|
328 | |
---|
329 | //////////////////////////////////////////////////////////////////////// |
---|
330 | // Building functions |
---|
331 | |
---|
332 | G4bool G4NURBS::MakeKnotVector(t_Dir & io_d, t_KnotVectorGenFlag in_KVGFlag) |
---|
333 | { |
---|
334 | G4bool isgood = (io_d.order + io_d.nbrCtrlPts == io_d.nbrKnots) |
---|
335 | && (io_d.pKnots == 0); |
---|
336 | if ( isgood ) |
---|
337 | { |
---|
338 | io_d.pKnots = new t_Knot [io_d.nbrKnots]; |
---|
339 | if (in_KVGFlag != UserDefined) |
---|
340 | { // let's do the knots |
---|
341 | t_indKnot indKnot = 0; |
---|
342 | t_index nbrCentralDistinctKnots = io_d.nbrCtrlPts-io_d.order; |
---|
343 | if ( (nbrCentralDistinctKnots % in_KVGFlag) == 0) |
---|
344 | { |
---|
345 | nbrCentralDistinctKnots /= in_KVGFlag; |
---|
346 | // first and last knots repeated 'order' Times |
---|
347 | for (t_index i=0; i < io_d.order; indKnot++,i++) |
---|
348 | { |
---|
349 | io_d.pKnots[indKnot] = 0; |
---|
350 | io_d.pKnots[indKnot+io_d.nbrCtrlPts] = 1; |
---|
351 | } |
---|
352 | |
---|
353 | t_Knot stepKnot = 1.0/(t_Knot)(nbrCentralDistinctKnots+1); |
---|
354 | t_Knot valKnot = stepKnot; |
---|
355 | |
---|
356 | // central knots |
---|
357 | for (t_indKnot j=0; j<nbrCentralDistinctKnots; valKnot+=stepKnot, j++) |
---|
358 | { |
---|
359 | for (t_indKnot k=0; k<t_indKnot(in_KVGFlag); indKnot++, k++) |
---|
360 | io_d.pKnots[indKnot] = valKnot; |
---|
361 | } |
---|
362 | } |
---|
363 | else isgood = false; |
---|
364 | } // end of knots making |
---|
365 | } |
---|
366 | return isgood; |
---|
367 | } |
---|
368 | |
---|
369 | |
---|
370 | std::ostream & operator << (std::ostream & io_ostr, |
---|
371 | G4NURBS::t_KnotVectorGenFlag in_f) |
---|
372 | { |
---|
373 | switch (in_f) |
---|
374 | { |
---|
375 | case G4NURBS::UserDefined: io_ostr << "UserDefined"; break; |
---|
376 | case G4NURBS::Regular: io_ostr << "Regular"; break; |
---|
377 | case G4NURBS::RegularRep: io_ostr << "RegularRep"; break; |
---|
378 | default: io_ostr << (G4int)in_f; |
---|
379 | } |
---|
380 | return io_ostr; |
---|
381 | } |
---|
382 | |
---|
383 | //////////////////////////////////////////////////////////////////////// |
---|
384 | // Constructors and co |
---|
385 | |
---|
386 | void G4NURBS::Conscheck() const |
---|
387 | { |
---|
388 | G4int dummy; |
---|
389 | t_direction dir; |
---|
390 | for (dummy=0; (dummy?(dir=V):(dir=U)),(dummy < NofD); dummy++) |
---|
391 | { |
---|
392 | if (m[dir].order<=0) |
---|
393 | { |
---|
394 | G4cerr << "\nFATAL ERROR: G4NURBS::G4NURBS: The order in the " |
---|
395 | << G4NURBS::Tochar(dir) |
---|
396 | << " direction must be >= 1" << G4endl; |
---|
397 | G4Exception("ERROR - G4NURBS::Conscheck()"); |
---|
398 | } |
---|
399 | if (m[dir].nbrCtrlPts<=0) |
---|
400 | { |
---|
401 | G4cerr << "\nFATAL ERROR: G4NURBS::G4NURBS: The number of control points " |
---|
402 | << G4NURBS::Tochar(dir) |
---|
403 | << " direction must be >= 1" << G4endl; |
---|
404 | G4Exception("ERROR - G4NURBS::Conscheck()"); |
---|
405 | } |
---|
406 | } // end of dummy |
---|
407 | } |
---|
408 | |
---|
409 | G4NURBS::G4NURBS ( t_order in_Uorder, t_order in_Vorder, |
---|
410 | t_inddCtrlPt in_UnbrCtrlPts, t_inddCtrlPt in_VnbrCtrlPts, |
---|
411 | t_CtrlPt * in_pCtrlPts, |
---|
412 | t_Knot * in_pUKnots, t_Knot * in_pVKnots, |
---|
413 | t_CheckFlag in_CheckFlag ) |
---|
414 | { |
---|
415 | m[U].order=in_Uorder; m[V].order=in_Vorder; |
---|
416 | m[U].nbrCtrlPts=in_UnbrCtrlPts; m[V].nbrCtrlPts=in_VnbrCtrlPts; |
---|
417 | |
---|
418 | mtotnbrCtrlPts = m[U].nbrCtrlPts * m[V].nbrCtrlPts; |
---|
419 | m[U].nbrKnots = m[U].order + m[U].nbrCtrlPts; |
---|
420 | m[V].nbrKnots = m[V].order + m[V].nbrCtrlPts; |
---|
421 | |
---|
422 | if (in_CheckFlag) |
---|
423 | Conscheck(); |
---|
424 | |
---|
425 | // CtrlPts |
---|
426 | if (! (mpCtrlPts = in_pCtrlPts) ) |
---|
427 | { |
---|
428 | G4cerr << "\nFATAL ERROR: G4NURBS::G4NURBS: " |
---|
429 | << "A NURBS MUST HAVE CONTROL POINTS!\n" |
---|
430 | << "\teven if they are defined later, the array must be allocated." |
---|
431 | << G4endl; |
---|
432 | G4Exception("ERROR - G4NURBS::G4NURBS()"); |
---|
433 | } |
---|
434 | //mnbralias = 0; |
---|
435 | |
---|
436 | // Knots |
---|
437 | t_direction dir; |
---|
438 | G4int dummy; |
---|
439 | for (dummy=0; (dummy?(dir=V):(dir=U)),(dummy < NofD); dummy++) |
---|
440 | { |
---|
441 | if ( !(m[dir].pKnots = (dummy?in_pVKnots:in_pUKnots)) ) |
---|
442 | { // make some regular knots between 0 & 1 |
---|
443 | if(!MakeKnotVector(m[dir], Regular)) |
---|
444 | { |
---|
445 | G4cerr << "\nFATAL ERROR: G4NURBS::G4NURBS: " |
---|
446 | << "Unable to make a Regular knot vector along " |
---|
447 | << G4NURBS::Tochar(dir) |
---|
448 | << " direction." |
---|
449 | << G4endl; |
---|
450 | G4Exception("ERROR - G4NURBS::G4NURBS()"); |
---|
451 | } |
---|
452 | //m[dir].nbralias = 0; |
---|
453 | } // end of knots-making |
---|
454 | } // end for dummy |
---|
455 | } // end of G4NURBS::G4NURBS |
---|
456 | |
---|
457 | // second constructor |
---|
458 | |
---|
459 | G4NURBS::G4NURBS( t_order in_Uorder, t_order in_Vorder, |
---|
460 | t_inddCtrlPt in_UnbrCtrlPts, t_inddCtrlPt in_VnbrCtrlPts, |
---|
461 | t_KnotVectorGenFlag in_UKVGFlag, |
---|
462 | t_KnotVectorGenFlag in_VKVGFlag, |
---|
463 | t_CheckFlag in_CheckFlag ) |
---|
464 | { |
---|
465 | m[U].order=in_Uorder; m[V].order=in_Vorder; |
---|
466 | m[U].nbrCtrlPts=in_UnbrCtrlPts; m[V].nbrCtrlPts=in_VnbrCtrlPts; |
---|
467 | |
---|
468 | mtotnbrCtrlPts = m[U].nbrCtrlPts * m[V].nbrCtrlPts; |
---|
469 | m[U].nbrKnots = m[U].order + m[U].nbrCtrlPts; |
---|
470 | m[V].nbrKnots = m[V].order + m[V].nbrCtrlPts; |
---|
471 | |
---|
472 | if (in_CheckFlag) |
---|
473 | Conscheck(); |
---|
474 | |
---|
475 | // Allocate CtrlPts |
---|
476 | mpCtrlPts = new t_CtrlPt [mtotnbrCtrlPts]; |
---|
477 | //mnbralias = 0; |
---|
478 | |
---|
479 | // Knots |
---|
480 | t_direction dir; |
---|
481 | G4int dummy; |
---|
482 | for (dummy=0; (dummy?(dir=V):(dir=U)),(dummy < NofD); dummy++) |
---|
483 | { |
---|
484 | t_KnotVectorGenFlag flag = (dummy?in_VKVGFlag:in_UKVGFlag); |
---|
485 | m[dir].pKnots = 0; // (allocation under our control) |
---|
486 | if ( flag && !MakeKnotVector(m[dir], flag) ) |
---|
487 | { |
---|
488 | G4cerr << "\nFATAL ERROR: G4NURBS::G4NURBS: " |
---|
489 | << "Unable to make knot vector along " |
---|
490 | << G4NURBS::Tochar(dir) |
---|
491 | << " direction. (" << m[dir].nbrKnots |
---|
492 | << " knots requested for a " |
---|
493 | << flag |
---|
494 | << " knots vector)" |
---|
495 | << G4endl; |
---|
496 | G4Exception("ERROR - G4NURBS::G4NURBS()"); |
---|
497 | } |
---|
498 | //m[dir].nbralias = 0; |
---|
499 | } |
---|
500 | } |
---|
501 | |
---|
502 | G4NURBS::G4NURBS(const G4NURBS & in_krNurb) |
---|
503 | : G4Visible(in_krNurb) |
---|
504 | { |
---|
505 | // we assume the in nurbs is ok |
---|
506 | |
---|
507 | // the number of CtrlPts can be copied straightly |
---|
508 | mtotnbrCtrlPts = in_krNurb.mtotnbrCtrlPts; |
---|
509 | |
---|
510 | // the main datas |
---|
511 | |
---|
512 | // but as m is an array of t_Dir and as t_Dir |
---|
513 | // is just a structure and not a class with a copy cons |
---|
514 | // whe need to duplicate the knots |
---|
515 | t_direction dir; |
---|
516 | G4int dummy; |
---|
517 | for (dummy=0; (dummy?(dir=V):(dir=U)),(dummy < NofD); dummy++) |
---|
518 | { |
---|
519 | // first we do a 'stupid' copy of m[dir] |
---|
520 | m[dir] = in_krNurb.m[dir]; |
---|
521 | // but as m is an array of t_Dir and as t_Dir |
---|
522 | // is just a structure and not a class with a copy cons |
---|
523 | // whe need to duplicate the knots |
---|
524 | m[dir].pKnots = new G4double [m[dir].nbrKnots]; |
---|
525 | // we copy the knots with memcpy. This function should be the fastest |
---|
526 | memcpy(m[dir].pKnots, in_krNurb.m[dir].pKnots, |
---|
527 | m[dir].nbrKnots * sizeof(G4double)); |
---|
528 | } // end of dummy loop |
---|
529 | |
---|
530 | // the control points |
---|
531 | // once again we need to do the copy |
---|
532 | mpCtrlPts = new t_CtrlPt [mtotnbrCtrlPts]; |
---|
533 | memcpy(mpCtrlPts, in_krNurb.mpCtrlPts, mtotnbrCtrlPts*sizeof(t_CtrlPt)); |
---|
534 | |
---|
535 | // and as it's very strange to copy a nurbs in G4 |
---|
536 | // we issue a warning : |
---|
537 | G4cerr << "\nWARNING: G4NURBS::G4NURBS(const G4NURBS &) used" << G4endl; |
---|
538 | } |
---|
539 | |
---|
540 | G4NURBS::~G4NURBS() |
---|
541 | { |
---|
542 | // we must free the two knots vector |
---|
543 | t_direction dir; |
---|
544 | G4int dummy; |
---|
545 | for (dummy=0; (dummy?(dir=V):(dir=U)),(dummy < NofD); dummy++) |
---|
546 | { |
---|
547 | if (m[dir].pKnots) |
---|
548 | delete m[dir].pKnots; // [m[dir].nbrKnots] if t_Knot become a class |
---|
549 | m[dir].pKnots = 0; |
---|
550 | } |
---|
551 | // now we free the CtrlPts array |
---|
552 | if (mpCtrlPts) |
---|
553 | delete [] mpCtrlPts; // [mtotnbrCtrlPts] if t_CtrlPt become a class |
---|
554 | mpCtrlPts = 0; |
---|
555 | } |
---|
556 | |
---|
557 | /************************************************************************ |
---|
558 | * * |
---|
559 | * Return the current knot the parameter u is less than or equal to. * |
---|
560 | * Find this "breakpoint" allows the evaluation routines to concentrate * |
---|
561 | * on only those control points actually effecting the curve around u.] * |
---|
562 | * * |
---|
563 | * m is the number of points on the curve (or surface direction) * |
---|
564 | * k is the order of the curve (or surface direction) * |
---|
565 | * kv is the knot vector ([0..m+k-1]) to find the break point in. * |
---|
566 | * * |
---|
567 | ************************************************************************/ |
---|
568 | static G4int FindBreakPoint(G4double u, const Float *kv, G4int m, G4int k) |
---|
569 | { |
---|
570 | G4int i; |
---|
571 | if (u == kv[m+1]) return m; /* Special case for closed interval */ |
---|
572 | i = m + k; |
---|
573 | while ((u < kv[i]) && (i > 0)) i--; |
---|
574 | return(i); |
---|
575 | } |
---|
576 | |
---|
577 | /************************************************************************ |
---|
578 | * * |
---|
579 | * Compute Bi,k(u), for i = 0..k. * |
---|
580 | * u the parameter of the spline to find the basis functions for* |
---|
581 | * brkPoint the start of the knot interval ("segment") * |
---|
582 | * kv the knot vector * |
---|
583 | * k the order of the curve * |
---|
584 | * bvals the array of returned basis values. * |
---|
585 | * * |
---|
586 | * (From Bartels, Beatty & Barsky, p.387) * |
---|
587 | * * |
---|
588 | ************************************************************************/ |
---|
589 | static void BasisFunctions(G4double u, G4int brkPoint, |
---|
590 | const Float *kv, G4int k, G4double *bvals) |
---|
591 | { |
---|
592 | G4int r, s, i; |
---|
593 | G4double omega; |
---|
594 | |
---|
595 | bvals[0] = 1.0; |
---|
596 | for (r=2; r <= k; r++) |
---|
597 | { |
---|
598 | i = brkPoint - r + 1; |
---|
599 | bvals[r-1] = 0.0; |
---|
600 | for (s=r-2; s >= 0; s--) |
---|
601 | { |
---|
602 | i++; |
---|
603 | if (i < 0) |
---|
604 | { |
---|
605 | omega = 0.0; |
---|
606 | } |
---|
607 | else |
---|
608 | { |
---|
609 | omega = (u - kv[i]) / (kv[i+r-1] - kv[i]); |
---|
610 | } |
---|
611 | bvals[s+1] = bvals[s+1] + (1.0-omega) * bvals[s]; |
---|
612 | bvals[s] = omega * bvals[s]; |
---|
613 | } |
---|
614 | } |
---|
615 | } |
---|
616 | |
---|
617 | /************************************************************************ |
---|
618 | * * |
---|
619 | * Compute derivatives of the basis functions Bi,k(u)' * |
---|
620 | * * |
---|
621 | ************************************************************************/ |
---|
622 | static void BasisDerivatives(G4double u, G4int brkPoint, |
---|
623 | const Float *kv, G4int k, G4double *dvals) |
---|
624 | { |
---|
625 | G4int s, i; |
---|
626 | G4double omega, knotScale; |
---|
627 | |
---|
628 | BasisFunctions(u, brkPoint, kv, k-1, dvals); |
---|
629 | |
---|
630 | dvals[k-1] = 0.0; /* BasisFunctions misses this */ |
---|
631 | |
---|
632 | knotScale = kv[brkPoint+1] - kv[brkPoint]; |
---|
633 | |
---|
634 | i = brkPoint - k + 1; |
---|
635 | for (s=k-2; s >= 0; s--) |
---|
636 | { |
---|
637 | i++; |
---|
638 | omega = knotScale * ((G4double)(k-1)) / (kv[i+k-1] - kv[i]); |
---|
639 | dvals[s+1] += -omega * dvals[s]; |
---|
640 | dvals[s] *= omega; |
---|
641 | } |
---|
642 | } |
---|
643 | |
---|
644 | /*********************************************************************** |
---|
645 | * * |
---|
646 | * Calculate a point p on NurbSurface n at a specific u, v * |
---|
647 | * using the tensor product. * |
---|
648 | * * |
---|
649 | * Note the valid parameter range for u and v is * |
---|
650 | * (kvU[orderU] <= u < kvU[numU), (kvV[orderV] <= v < kvV[numV]) * |
---|
651 | * * |
---|
652 | ***********************************************************************/ |
---|
653 | void G4NURBS::CalcPoint(G4double u, G4double v, G4Point3D &p, |
---|
654 | G4Vector3D &utan, G4Vector3D &vtan) const |
---|
655 | { |
---|
656 | #define MAXORDER 50 |
---|
657 | struct Point4 |
---|
658 | { |
---|
659 | G4double x, y, z, w; |
---|
660 | }; |
---|
661 | |
---|
662 | G4int i, j, ri, rj; |
---|
663 | G4int ubrkPoint, ufirst; |
---|
664 | G4double bu[MAXORDER], buprime[MAXORDER]; |
---|
665 | G4int vbrkPoint, vfirst; |
---|
666 | G4double bv[MAXORDER], bvprime[MAXORDER]; |
---|
667 | Point4 r, rutan, rvtan; |
---|
668 | |
---|
669 | r.x = 0.0; r.y = 0.0; r.z = 0.0; r.w = 0.0; |
---|
670 | rutan = r; rvtan = r; |
---|
671 | |
---|
672 | G4int numU = GetUnbrCtrlPts(); |
---|
673 | G4int numV = GetVnbrCtrlPts(); |
---|
674 | G4int orderU = GetUorder(); |
---|
675 | G4int orderV = GetVorder(); |
---|
676 | |
---|
677 | /* Evaluate non-uniform basis functions (and derivatives) */ |
---|
678 | |
---|
679 | ubrkPoint = FindBreakPoint(u, m[U].pKnots, numU-1, orderU); |
---|
680 | ufirst = ubrkPoint - orderU + 1; |
---|
681 | BasisFunctions (u, ubrkPoint, m[U].pKnots, orderU, bu); |
---|
682 | BasisDerivatives(u, ubrkPoint, m[U].pKnots, orderU, buprime); |
---|
683 | |
---|
684 | vbrkPoint = FindBreakPoint(v, m[V].pKnots, numV-1, orderV); |
---|
685 | vfirst = vbrkPoint - orderV + 1; |
---|
686 | BasisFunctions (v, vbrkPoint, m[V].pKnots, orderV, bv); |
---|
687 | BasisDerivatives(v, vbrkPoint, m[V].pKnots, orderV, bvprime); |
---|
688 | |
---|
689 | /* Weight control points against the basis functions */ |
---|
690 | |
---|
691 | t_doubleCtrlPt *cpoint; |
---|
692 | Point4 cp; |
---|
693 | G4double tmp; |
---|
694 | |
---|
695 | for (i=0; i<orderV; i++) |
---|
696 | { |
---|
697 | for (j=0; j<orderU; j++) |
---|
698 | { |
---|
699 | ri = orderV - 1 - i; |
---|
700 | rj = orderU - 1 - j; |
---|
701 | |
---|
702 | tmp = bu[rj] * bv[ri]; |
---|
703 | cpoint = GetdoubleCtrlPt(j+ufirst, i+vfirst); |
---|
704 | cp.x = *cpoint[G4NURBS::X]; |
---|
705 | cp.y = *cpoint[G4NURBS::Y]; |
---|
706 | cp.z = *cpoint[G4NURBS::Z]; |
---|
707 | cp.w = *cpoint[G4NURBS::W]; |
---|
708 | r.x += cp.x * tmp; |
---|
709 | r.y += cp.y * tmp; |
---|
710 | r.z += cp.z * tmp; |
---|
711 | r.w += cp.w * tmp; |
---|
712 | |
---|
713 | tmp = buprime[rj] * bv[ri]; |
---|
714 | rutan.x += cp.x * tmp; |
---|
715 | rutan.y += cp.y * tmp; |
---|
716 | rutan.z += cp.z * tmp; |
---|
717 | rutan.w += cp.w * tmp; |
---|
718 | |
---|
719 | tmp = bu[rj] * bvprime[ri]; |
---|
720 | rvtan.x += cp.x * tmp; |
---|
721 | rvtan.y += cp.y * tmp; |
---|
722 | rvtan.z += cp.z * tmp; |
---|
723 | rvtan.w += cp.w * tmp; |
---|
724 | } |
---|
725 | } |
---|
726 | |
---|
727 | /* Project tangents, using the quotient rule for differentiation */ |
---|
728 | |
---|
729 | G4double wsqrdiv = 1.0 / (r.w * r.w); |
---|
730 | |
---|
731 | utan.setX((r.w * rutan.x - rutan.w * r.x) * wsqrdiv); |
---|
732 | utan.setY((r.w * rutan.y - rutan.w * r.y) * wsqrdiv); |
---|
733 | utan.setZ((r.w * rutan.z - rutan.w * r.z) * wsqrdiv); |
---|
734 | |
---|
735 | vtan.setX((r.w * rvtan.x - rvtan.w * r.x) * wsqrdiv); |
---|
736 | vtan.setY((r.w * rvtan.y - rvtan.w * r.y) * wsqrdiv); |
---|
737 | vtan.setZ((r.w * rvtan.z - rvtan.w * r.z) * wsqrdiv); |
---|
738 | |
---|
739 | p.setX(r.x / r.w); |
---|
740 | p.setY(r.y / r.w); |
---|
741 | p.setZ(r.z / r.w); |
---|
742 | } |
---|