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

Last change on this file since 1346 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 24.7 KB
Line 
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-04-beta-01 $
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.