source: trunk/source/graphics_reps/src/G4NURBS.cc@ 1235

Last change on this file since 1235 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 24.7 KB
RevLine 
[830]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 $
[1228]28// GEANT4 tag $Name: geant4-09-03 $
[830]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
47std::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
109G4float 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
124G4double 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
140G4NURBS::t_floatCtrlPt*
141G4NURBS::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
154G4NURBS::t_floatCtrlPt*
155G4NURBS::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
170G4NURBS::t_doubleCtrlPt*
171G4NURBS::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
184G4NURBS::t_doubleCtrlPt*
185G4NURBS::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
201G4float * 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
210G4double * 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
219G4float * 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
227G4double * 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
237G4NURBS::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
256G4bool G4NURBS::KnotsIterator::pick(G4double * inout_pDbl)
257{
258 (*inout_pDbl) = (G4double)(*mp);
259 return (G4bool)((++mp)<kmpMax);
260}
261
262G4bool G4NURBS::KnotsIterator::pick(G4float * inout_pFlt)
263{
264 (*inout_pFlt) = (G4float)(*mp);
265 return (G4bool)((++mp)<kmpMax);
266}
267
268G4NURBS::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
286G4bool G4NURBS::CtrlPtsCoordsIterator::pick(G4double * inout_pDbl)
287{
288 (*inout_pDbl) = (G4double)((*mp));
289 return (G4bool)((++mp)<kmpMax);
290}
291
292G4bool G4NURBS::CtrlPtsCoordsIterator::pick(G4float * inout_pFlt)
293{
294 (*inout_pFlt) = (G4float)((*mp));
295 return (G4bool)((++mp)<kmpMax);
296}
297
298G4NURBS::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
315G4bool 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
322G4bool 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
332G4bool 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
370std::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
386void 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
409G4NURBS::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
459G4NURBS::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
502G4NURBS::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
540G4NURBS::~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 ************************************************************************/
568static 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 ************************************************************************/
589static 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 ************************************************************************/
622static 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 ***********************************************************************/
653void 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}
Note: See TracBrowser for help on using the repository browser.