G4NURBS.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // 
00030 // Olivier Crumeyrolle  12 September 1996
00031 
00032 // G4NURBS.cc
00033 // Implementation of class G4NURBS
00034 // OC 100796
00035 
00036 
00037 #include "G4NURBS.hh"
00038 
00039 // G4NURBS.hh includes globals.hh which includes a lot of others
00040 // so no more includes required here
00041 
00043 //    Here start the real world. Please, check your armored jacket.   //
00045 
00046 std::ostream & operator << (std::ostream & inout_outStream,
00047                               const G4NURBS & in_kNurb)
00048 {
00049   inout_outStream
00050   // the magic could be changed for good reasons only
00051       << "##ojc{NURBS}def[1.01.96.7]   Just a magic. Could be added to /etc/magic"
00052       << "\n# NURBS Definition File (human and computer readable format)"
00053       << "\n# :" << in_kNurb.Whoami()
00054       << "\n# U order\tV order : " 
00055       << '\n' << in_kNurb.GetUorder() << "\t\t" << in_kNurb.GetVorder();
00056   // number of knots and knots themselves for U and V
00057   for (G4NURBS::t_direction dir = G4NURBS::U; dir < G4NURBS::NofD;
00058        /*(*(G4int *)(&dir))++*/ dir=(G4NURBS::t_direction)(((G4int)(dir))+1) )
00059   {
00060     inout_outStream
00061       << "\n# Number of knots along " << G4NURBS::Tochar(dir)
00062       << '\n' << in_kNurb.GetnbrKnots(dir)
00063       << "\n# " << G4NURBS::Tochar(dir) << " knots vector (as a column)";
00064     {  // begin knots iteration
00065        G4double oneKnot;
00066        G4NURBS::KnotsIterator knotI(in_kNurb,dir);
00067        G4bool otherKnots;
00068        do 
00069        {
00070          otherKnots = knotI.pick(&oneKnot);
00071          inout_outStream << "\n\t\t" << oneKnot;
00072        }
00073        while (otherKnots);
00074     }  // end of knots iteration
00075   }  // end of direction loop
00076 
00077   // number of control points in U and V direction
00078   // and controlpoints
00079   inout_outStream
00080       << "\n# Number of control points along U and V"
00081       << '\n' << in_kNurb.GetUnbrCtrlPts() 
00082       << "   " << in_kNurb.GetVnbrCtrlPts()
00083       << "\n# Control Points (one by line, U increasing first)";
00084   { // begin of control points iteration
00085     G4NURBS::t_doubleCtrlPt   oneCP;
00086     G4NURBS::CtrlPtsIterator  cpI(in_kNurb);
00087     G4bool otherCPs;
00088     do
00089     {
00090       otherCPs = cpI.pick(&oneCP);
00091       inout_outStream 
00092         << "\n\t" << oneCP[G4NURBS::X]
00093         << "\t" << oneCP[G4NURBS::Y]
00094         << "\t" << oneCP[G4NURBS::Z]
00095         << "\t" << oneCP[G4NURBS::W];
00096     }
00097     while (otherCPs);
00098   } // end of control point iteration
00099 
00100   inout_outStream << "\n# That's all!"
00101                   << G4endl;  // endl do an \n and a flush
00102   return inout_outStream;
00103 }
00104 
00105 // the CC compiler issue some "maybe no value returned"
00106 // but everything is ok
00107 
00108 G4float G4NURBS::GetfloatKnot(t_direction in_dir, t_indKnot in_index) const
00109 {
00110   in_dir = (t_direction)(in_dir & DMask);
00111   if ( in_index < m[in_dir].nbrKnots )
00112     return ((G4float)(m[in_dir].pKnots[in_index])); 
00113   else
00114   {
00115     G4cerr << "\nERROR: G4NURBS::GetfloatKnot: index out of range\n"
00116            << "\n\t in_dir : " << G4int(in_dir)
00117            << ", in_index : " << G4int(in_index)
00118            << "m[in_dir].nbrKnots : " << m[in_dir].nbrKnots << G4endl;
00119     return ((G4float)m[in_dir].pKnots[m[in_dir].nbrKnots-1]); 
00120   }
00121 }
00122  
00123 G4double G4NURBS::GetdoubleKnot(t_direction in_dir, t_indKnot in_index) const
00124 {
00125   in_dir = (t_direction)(in_dir & DMask);
00126   if ( in_index < m[in_dir].nbrKnots )
00127     return (G4double)(m[in_dir].pKnots[in_index]);
00128   else
00129   {
00130     G4cerr << "\nERROR: G4NURBS::GetdoubleKnot: index out of range"
00131            << "\n\t in_dir : " << G4int(in_dir)
00132            << ", in_index : " << G4int(in_index)
00133            << "m[in_dir].nbrKnots : " << m[in_dir].nbrKnots
00134            << G4endl;
00135     return (G4double)(m[in_dir].pKnots[m[in_dir].nbrKnots-1]); 
00136   }
00137 }
00138 
00139 G4NURBS::t_floatCtrlPt*
00140 G4NURBS::GetfloatCtrlPt(t_indCtrlPt in_onedimindex) const
00141 {
00142   if (in_onedimindex < mtotnbrCtrlPts)
00143     return TofloatCtrlPt(mpCtrlPts[in_onedimindex]);       
00144   else
00145   {
00146     G4cerr << "\nERROR: G4NURBS::GetfloatCtrlPt: index out of range"
00147            << "\n\t in_onedimindex : " << in_onedimindex
00148            << " , mtotnbrCtrlPts : " << mtotnbrCtrlPts << G4endl;
00149     return TofloatCtrlPt(mpCtrlPts[mtotnbrCtrlPts-1]);  
00150   }    
00151 }
00152 
00153 G4NURBS::t_floatCtrlPt*
00154 G4NURBS::GetfloatCtrlPt(t_inddCtrlPt in_Uindex, t_inddCtrlPt in_Vindex) const
00155 {
00156   if ( (in_Uindex < m[U].nbrCtrlPts) && (in_Vindex < m[V].nbrCtrlPts) )
00157     return TofloatCtrlPt(mpCtrlPts[To1d(in_Uindex, in_Vindex)]);
00158   else
00159   {
00160     G4cerr << "\nERROR: G4NURBS::GetfloatCtrlPt: index(s) out of range"
00161            << "\n\t in_Uindex : " << in_Uindex
00162            << " , in_Vindex : " << in_Vindex
00163            << " , UnbrCtrlPts : " << m[U].nbrCtrlPts
00164            << " , VnbrCtrlPts : " << m[V].nbrCtrlPts << G4endl;
00165     return TofloatCtrlPt(mpCtrlPts[mtotnbrCtrlPts-1]);  
00166   }   
00167 }
00168 
00169 G4NURBS::t_doubleCtrlPt*
00170 G4NURBS::GetdoubleCtrlPt(t_indCtrlPt in_onedimindex) const
00171 {
00172   if ( in_onedimindex < mtotnbrCtrlPts )
00173     return TodoubleCtrlPt(mpCtrlPts[in_onedimindex]);
00174   else
00175   {
00176     G4cerr << "\nERROR: G4NURBS::getdoubleCtrlPts: index out of range"
00177            << "\n\t in_onedimindex : " << in_onedimindex
00178            << " , mtotnbrCtrlPts : " << mtotnbrCtrlPts << G4endl;
00179     return TodoubleCtrlPt(mpCtrlPts[mtotnbrCtrlPts-1]);  
00180   }   
00181 }
00182 
00183 G4NURBS::t_doubleCtrlPt*
00184 G4NURBS::GetdoubleCtrlPt(t_inddCtrlPt in_Uindex, t_inddCtrlPt in_Vindex) const
00185 {
00186   if ( (in_Uindex < m[U].nbrCtrlPts) && (in_Vindex < m[V].nbrCtrlPts) )
00187     return TodoubleCtrlPt(mpCtrlPts[To1d(in_Uindex, in_Vindex)]);
00188   else
00189   {
00190     G4cerr << "\nERROR: G4NURBS::GetdoubleCtrlPt: index(s) out of range"
00191            << "\n\t in_Uindex : " << in_Uindex
00192            << " , in_Vindex : " << in_Vindex
00193            << " , UnbrCtrlPts : " << m[U].nbrCtrlPts
00194            << " , VnbrCtrlPts : " << m[V].nbrCtrlPts << G4endl;
00195     return TodoubleCtrlPt(mpCtrlPts[mtotnbrCtrlPts-1]);  
00196   }  
00197 }
00198  
00199 // Total copy
00200 G4float * G4NURBS::GetfloatAllKnots(t_direction in_dir) const
00201 {
00202   in_dir = (t_direction)(in_dir & DMask);
00203   G4float * p = new G4float [m[in_dir].nbrKnots];
00204   for (t_indKnot i = 0; i < m[in_dir].nbrKnots; i++)
00205     p[i] = (G4float)m[in_dir].pKnots[i];
00206   return p;
00207 }
00208 
00209 G4double * G4NURBS::GetdoubleAllKnots(t_direction in_dir) const
00210 {
00211   in_dir = (t_direction)(in_dir & DMask);
00212   G4double * p = new G4double [m[in_dir].nbrKnots];
00213   for (t_indKnot i = 0; i < m[in_dir].nbrKnots; i++)
00214     p[i] = (G4double)m[in_dir].pKnots[i];
00215   return p;
00216 }
00217 
00218 G4float * G4NURBS::GetfloatAllCtrlPts() const
00219 {
00220   G4float * p = new G4float [mtotnbrCtrlPts*NofC];
00221   for (t_indKnot i = 0; i < mtotnbrCtrlPts*NofC; i++)
00222     p[i] = (G4float)(((t_Coord *)mpCtrlPts)[i]);
00223   return p;
00224 }
00225 
00226 G4double * G4NURBS::GetdoubleAllCtrlPts() const
00227 {
00228   G4double * p = new G4double [mtotnbrCtrlPts*NofC];
00229   for (t_indKnot i = 0; i < mtotnbrCtrlPts*NofC; i++)
00230     p[i] = (G4double)(((t_Coord *)mpCtrlPts)[i]);
00231   return p;
00232 }
00233 
00234 // Iterators
00235 
00236 G4NURBS::KnotsIterator::KnotsIterator(const G4NURBS & in_rNurb,
00237                                       G4NURBS::t_direction in_dir,
00238                                       t_indKnot in_startIndex)
00239  : kmdir((G4NURBS::t_direction)(in_dir &  G4NURBS::DMask)),
00240    kmpMax(in_rNurb.m[kmdir].pKnots + in_rNurb.m[kmdir].nbrKnots)
00241 {
00242   if (in_startIndex < in_rNurb.m[kmdir].nbrKnots)
00243     mp = in_rNurb.m[kmdir].pKnots + in_startIndex;
00244   else
00245   {
00246     G4cerr   << "\nERROR: G4NURBS::KnotsIterator: in_startIndex out of range"
00247              << "\n\tin_startIndex : " << in_startIndex
00248              << ", nbr of knots : " << in_rNurb.m[kmdir].nbrKnots
00249              << "\n\t mp set to NULL, calls to picking functions will fail"
00250              << G4endl;
00251     mp = 0;
00252   }
00253 }
00254 
00255 G4bool G4NURBS::KnotsIterator::pick(G4double * inout_pDbl)
00256 {
00257   (*inout_pDbl) = (G4double)(*mp);
00258   return (G4bool)((++mp)<kmpMax);
00259 }
00260 
00261 G4bool G4NURBS::KnotsIterator::pick(G4float * inout_pFlt)
00262 {
00263   (*inout_pFlt) = (G4float)(*mp);
00264   return (G4bool)((++mp)<kmpMax);
00265 }
00266 
00267 G4NURBS::CtrlPtsCoordsIterator::CtrlPtsCoordsIterator(const G4NURBS & in_rNurb,
00268                                                t_indCtrlPt in_startCtrlPtIndex)
00269  : kmpMax((const t_Coord *)(in_rNurb.mpCtrlPts + in_rNurb.mtotnbrCtrlPts))
00270 {
00271   if (in_startCtrlPtIndex < in_rNurb.mtotnbrCtrlPts )
00272     mp = (const t_Coord *)(in_rNurb.mpCtrlPts + in_startCtrlPtIndex);
00273   else
00274   {
00275     G4cerr << "\nERROR: G4NURBS::CtrlPtsCoordsIterator: "
00276            << "in_startCtrlPtIndex out of range"
00277            << "\n\tin_startCtrlPtIndex : " << in_startCtrlPtIndex
00278            << ", nbr of CtrlPts : " << in_rNurb.mtotnbrCtrlPts 
00279            << "\n\t mp set to NULL, calls to picking functions will fail"
00280            << G4endl;
00281     mp = 0;
00282   }
00283 }
00284 
00285 G4bool G4NURBS::CtrlPtsCoordsIterator::pick(G4double  * inout_pDbl)
00286 {
00287   (*inout_pDbl) = (G4double)((*mp));
00288   return (G4bool)((++mp)<kmpMax);
00289 }
00290 
00291 G4bool G4NURBS::CtrlPtsCoordsIterator::pick(G4float * inout_pFlt)
00292 {
00293   (*inout_pFlt) = (G4float)((*mp));
00294   return (G4bool)((++mp)<kmpMax);
00295 }
00296 
00297 G4NURBS::CtrlPtsIterator::CtrlPtsIterator(const G4NURBS & in_rNurb,
00298                                           t_indCtrlPt in_startIndex)
00299  : kmpMax(in_rNurb.mpCtrlPts + in_rNurb.mtotnbrCtrlPts)
00300 {
00301   if (in_startIndex < in_rNurb.mtotnbrCtrlPts )
00302     mp = (in_rNurb.mpCtrlPts + in_startIndex);
00303   else
00304   {
00305     G4cerr << "\nERROR: G4NURBS::CtrlPtsIterator: in_startIndex out of range"
00306            << "\n\tin_startIndex : " << in_startIndex
00307            << ", nbr of CtrlPts : " << in_rNurb.mtotnbrCtrlPts 
00308            << "\n\t mp set to NULL, calls to picking functions will fail"
00309            << G4endl;
00310     mp = 0;
00311   }
00312 }
00313 
00314 G4bool G4NURBS::CtrlPtsIterator::pick(t_doubleCtrlPt * inout_pDblCtrlPt)
00315 {
00316   for (t_indCoord i = G4NURBS::X; i < G4NURBS::NofC; i++)
00317     (*inout_pDblCtrlPt)[i] = (G4double)((*mp)[i]);
00318   return (G4bool)((++mp)<kmpMax);
00319 }
00320 
00321 G4bool G4NURBS::CtrlPtsIterator::pick(t_floatCtrlPt * inout_pFltCtrlPt)
00322 {
00323   for (t_indCoord i = G4NURBS::X; i < G4NURBS::NofC; i++)
00324     (*inout_pFltCtrlPt)[i] = (G4float)((*mp)[i]);
00325   return (G4bool)((++mp)<kmpMax);
00326 }
00327 
00329 // Building functions
00330 
00331 G4bool G4NURBS::MakeKnotVector(t_Dir & io_d, t_KnotVectorGenFlag in_KVGFlag)
00332 {
00333   G4bool isgood = (io_d.order + io_d.nbrCtrlPts == io_d.nbrKnots)
00334                   && (io_d.pKnots == 0);
00335   if ( isgood )
00336   {
00337     io_d.pKnots = new t_Knot [io_d.nbrKnots];
00338     if (in_KVGFlag != UserDefined)
00339     {  // let's do the knots
00340       t_indKnot indKnot = 0;
00341       t_index nbrCentralDistinctKnots = io_d.nbrCtrlPts-io_d.order;
00342       if ( (nbrCentralDistinctKnots % in_KVGFlag) == 0)
00343       {
00344         nbrCentralDistinctKnots /= in_KVGFlag; 
00345         // first and last knots repeated 'order' Times
00346         for (t_index i=0; i < io_d.order; indKnot++,i++)
00347         {
00348           io_d.pKnots[indKnot] = 0;
00349           io_d.pKnots[indKnot+io_d.nbrCtrlPts] = 1;
00350         }
00351 
00352         t_Knot stepKnot = 1.0/(t_Knot)(nbrCentralDistinctKnots+1);
00353         t_Knot valKnot = stepKnot;
00354 
00355         // central knots
00356         for (t_indKnot j=0; j<nbrCentralDistinctKnots; valKnot+=stepKnot, j++)
00357         {
00358           for (t_indKnot k=0; k<t_indKnot(in_KVGFlag); indKnot++, k++)
00359             io_d.pKnots[indKnot] = valKnot;
00360         }
00361       }
00362       else isgood = false;
00363     } // end of knots making
00364   }  
00365   return isgood;
00366 }
00367 
00368 
00369 std::ostream & operator << (std::ostream & io_ostr,
00370                               G4NURBS::t_KnotVectorGenFlag in_f)
00371 {
00372   switch (in_f) 
00373   {
00374     case G4NURBS::UserDefined: io_ostr << "UserDefined"; break;
00375     case G4NURBS::Regular:     io_ostr << "Regular"; break;
00376     case G4NURBS::RegularRep:  io_ostr << "RegularRep"; break;
00377     default:                   io_ostr << (G4int)in_f;
00378   }
00379   return io_ostr;
00380 }
00381 
00383 // Constructors and co
00384 
00385 void G4NURBS::Conscheck() const
00386 {
00387   G4int dummy;
00388   t_direction dir;
00389   for (dummy=0; (dummy?(dir=V):(dir=U)),(dummy < NofD); dummy++)
00390   { 
00391     if (m[dir].order<=0)
00392     { 
00393       G4ExceptionDescription ed;
00394       ed << "The order in the "
00395          << G4NURBS::Tochar(dir) 
00396          << " direction must be >= 1" << G4endl;
00397       G4Exception("G4NURBS::Conscheck()",
00398                   "greps9001", FatalException, ed);
00399     }
00400     if (m[dir].nbrCtrlPts<=0)
00401     {
00402       G4ExceptionDescription ed;
00403       ed << "The number of control points "
00404              << G4NURBS::Tochar(dir) 
00405              << " direction must be >= 1" << G4endl;
00406       G4Exception("G4NURBS::Conscheck()",
00407                   "greps9002", FatalException, ed);
00408     }
00409   }  // end of dummy
00410 }
00411 
00412 G4NURBS::G4NURBS ( t_order in_Uorder, t_order in_Vorder,
00413                    t_inddCtrlPt in_UnbrCtrlPts, t_inddCtrlPt in_VnbrCtrlPts,
00414                    t_CtrlPt * in_pCtrlPts,
00415                    t_Knot * in_pUKnots, t_Knot * in_pVKnots,
00416                    t_CheckFlag in_CheckFlag )
00417 {
00418   m[U].order=in_Uorder; m[V].order=in_Vorder;
00419   m[U].nbrCtrlPts=in_UnbrCtrlPts; m[V].nbrCtrlPts=in_VnbrCtrlPts;
00420 
00421   mtotnbrCtrlPts = m[U].nbrCtrlPts * m[V].nbrCtrlPts;
00422   m[U].nbrKnots = m[U].order + m[U].nbrCtrlPts;
00423   m[V].nbrKnots = m[V].order + m[V].nbrCtrlPts;
00424     
00425   if (in_CheckFlag)
00426     Conscheck();
00427 
00428   // CtrlPts
00429   if (! (mpCtrlPts = in_pCtrlPts) )
00430   {
00431     G4ExceptionDescription ed;
00432     ed << "A NURBS MUST HAVE CONTROL POINTS!\n"
00433        << "\teven if they are defined later, the array must be allocated."
00434        << G4endl;
00435     G4Exception("G4NURBS::G4NURBS()",
00436                 "greps9003", FatalException, ed);
00437   }
00438   //mnbralias = 0;
00439 
00440   // Knots
00441   t_direction dir;
00442   G4int dummy;
00443   for (dummy=0; (dummy?(dir=V):(dir=U)),(dummy < NofD); dummy++)
00444   {
00445     if ( !(m[dir].pKnots = (dummy?in_pVKnots:in_pUKnots)) )
00446     {  // make some regular knots between 0 & 1
00447       if(!MakeKnotVector(m[dir], Regular))
00448       {
00449         G4ExceptionDescription ed;
00450         ed << "Unable to make a Regular knot vector along "
00451            << G4NURBS::Tochar(dir)
00452            << " direction."
00453            << G4endl;
00454         G4Exception("G4NURBS::G4NURBS()",
00455                     "greps9004", FatalException, ed);
00456       }
00457       //m[dir].nbralias = 0;
00458     }  // end of knots-making
00459   }  // end for dummy
00460 } // end of G4NURBS::G4NURBS 
00461 
00462 // second constructor
00463 
00464 G4NURBS::G4NURBS( t_order in_Uorder, t_order in_Vorder,
00465                   t_inddCtrlPt in_UnbrCtrlPts, t_inddCtrlPt in_VnbrCtrlPts,
00466                   t_KnotVectorGenFlag in_UKVGFlag,
00467                   t_KnotVectorGenFlag in_VKVGFlag,
00468                   t_CheckFlag in_CheckFlag )
00469 {
00470   m[U].order=in_Uorder; m[V].order=in_Vorder;
00471   m[U].nbrCtrlPts=in_UnbrCtrlPts; m[V].nbrCtrlPts=in_VnbrCtrlPts;
00472 
00473   mtotnbrCtrlPts = m[U].nbrCtrlPts * m[V].nbrCtrlPts;
00474   m[U].nbrKnots = m[U].order + m[U].nbrCtrlPts;
00475   m[V].nbrKnots = m[V].order + m[V].nbrCtrlPts;
00476     
00477   if (in_CheckFlag)
00478     Conscheck();
00479 
00480   // Allocate CtrlPts
00481   mpCtrlPts = new t_CtrlPt [mtotnbrCtrlPts];
00482   //mnbralias = 0;
00483 
00484   // Knots
00485   t_direction dir;
00486   G4int dummy;
00487   for (dummy=0; (dummy?(dir=V):(dir=U)),(dummy < NofD); dummy++)
00488   {
00489     t_KnotVectorGenFlag flag = (dummy?in_VKVGFlag:in_UKVGFlag);
00490     m[dir].pKnots = 0;  // (allocation under our control)
00491     if ( flag != UserDefined && !MakeKnotVector(m[dir], flag) )
00492     {
00493       G4ExceptionDescription ed;
00494       ed << "Unable to make knot vector along "
00495          << G4NURBS::Tochar(dir)
00496          << " direction. (" << m[dir].nbrKnots 
00497          << " knots requested for a " 
00498          << flag 
00499          << " knots vector)"
00500          << G4endl;
00501       G4Exception("G4NURBS::G4NURBS()",
00502                   "greps9005", FatalException, ed);
00503     }
00504     //m[dir].nbralias = 0;
00505   }
00506 }
00507 
00508 G4NURBS::G4NURBS(const G4NURBS & in_krNurb)
00509   : G4Visible(in_krNurb)
00510 {
00511   // we assume the in nurbs is ok
00512 
00513   // the number of CtrlPts can be copied straightly
00514   mtotnbrCtrlPts = in_krNurb.mtotnbrCtrlPts;
00515 
00516   // the main datas
00517 
00518   // but as m is an array of t_Dir and as t_Dir
00519   // is just a structure and not a class with a copy cons
00520   // whe need to duplicate the knots
00521   t_direction dir;
00522   G4int dummy;
00523   for (dummy=0; (dummy?(dir=V):(dir=U)),(dummy < NofD); dummy++)
00524   {
00525     // first we do a 'stupid' copy of m[dir]
00526     m[dir] = in_krNurb.m[dir];
00527     // but as m is an array of t_Dir and as t_Dir
00528     // is just a structure and not a class with a copy cons
00529     // whe need to duplicate the knots
00530     m[dir].pKnots = new G4double [m[dir].nbrKnots];
00531     // we copy the knots with memcpy. This function should be the fastest
00532     memcpy(m[dir].pKnots, in_krNurb.m[dir].pKnots,
00533            m[dir].nbrKnots * sizeof(G4double));
00534   }  // end of dummy loop
00535 
00536   // the control points
00537   // once again we need to do the copy
00538   mpCtrlPts = new t_CtrlPt [mtotnbrCtrlPts];
00539   memcpy(mpCtrlPts, in_krNurb.mpCtrlPts, mtotnbrCtrlPts*sizeof(t_CtrlPt)); 
00540 
00541   // and as it's very strange to copy a nurbs in G4
00542   // we issue a warning :
00543   G4cerr << "\nWARNING: G4NURBS::G4NURBS(const G4NURBS &) used" << G4endl;
00544 }
00545 
00546 G4NURBS::~G4NURBS()
00547 {
00548   // we must free the two knots vector
00549   t_direction dir;
00550   G4int dummy;
00551   for (dummy=0; (dummy?(dir=V):(dir=U)),(dummy < NofD); dummy++)
00552   {
00553     if (m[dir].pKnots)
00554       delete m[dir].pKnots;  // [m[dir].nbrKnots] if t_Knot become a class
00555     m[dir].pKnots = 0;
00556   }
00557   // now we free the CtrlPts array
00558   if (mpCtrlPts)
00559     delete [] mpCtrlPts;    // [mtotnbrCtrlPts] if t_CtrlPt become a class
00560   mpCtrlPts = 0;
00561 }
00562 
00563 /************************************************************************
00564  *                                                                      *
00565  * Return the current knot the parameter u is less than or equal to.    *
00566  * Find this "breakpoint" allows the evaluation routines to concentrate *
00567  * on only those control points actually effecting the curve around u.] *
00568  *                                                                      *
00569  *  np  is the number of points on the curve (or surface direction)     *
00570  *  k   is the order of the curve (or surface direction)                *
00571  *  kv  is the knot vector ([0..np+k-1]) to find the break point in.    *
00572  *                                                                      *
00573  ************************************************************************/
00574 static G4int FindBreakPoint(G4double u, const Float *kv, G4int np, G4int k)
00575 {
00576   G4int i;
00577   if (u == kv[np+1]) return np;      /* Special case for closed interval */
00578   i = np + k;
00579   while ((u < kv[i]) && (i > 0)) i--;
00580   return(i);
00581 }
00582 
00583 /************************************************************************
00584  *                                                                      *
00585  * Compute Bi,k(u), for i = 0..k.                                       *
00586  *  u        the parameter of the spline to find the basis functions for*
00587  *  brkPoint the start of the knot interval ("segment")                 *
00588  *  kv       the knot vector                                            *
00589  *  k        the order of the curve                                     *
00590  *  bvals    the array of returned basis values.                        *
00591  *                                                                      *
00592  * (From Bartels, Beatty & Barsky, p.387)                               *
00593  *                                                                      *
00594  ************************************************************************/
00595 static void BasisFunctions(G4double u, G4int brkPoint,
00596                            const Float *kv, G4int k, G4double *bvals)
00597 {
00598   G4int r, q, i;
00599   G4double omega;
00600 
00601   bvals[0] = 1.0;
00602   for (r=2; r <= k; r++)
00603   {
00604     i = brkPoint - r + 1;
00605     bvals[r-1] = 0.0;
00606     for (q=r-2; q >= 0; q--)
00607     {
00608       i++;
00609       if (i < 0)
00610       {
00611         omega = 0.0;
00612       }
00613       else
00614       {
00615         omega = (u - kv[i]) / (kv[i+r-1] - kv[i]);
00616       }
00617       bvals[q+1] = bvals[q+1] + (1.0-omega) * bvals[q];
00618       bvals[q]   = omega * bvals[q];
00619     }
00620   }
00621 }
00622 
00623 /************************************************************************
00624  *                                                                      *
00625  * Compute derivatives of the basis functions Bi,k(u)'                  *
00626  *                                                                      *
00627  ************************************************************************/
00628 static void BasisDerivatives(G4double u, G4int brkPoint,
00629                              const Float *kv, G4int k, G4double *dvals)
00630 {
00631   G4int q, i;
00632   G4double omega, knotScale;
00633 
00634   BasisFunctions(u, brkPoint, kv, k-1, dvals);
00635 
00636   dvals[k-1] = 0.0;      /* BasisFunctions misses this */
00637 
00638   knotScale = kv[brkPoint+1] - kv[brkPoint];
00639 
00640   i = brkPoint - k + 1;
00641   for (q=k-2; q >= 0; q--)
00642   {
00643     i++;
00644     omega = knotScale * ((G4double)(k-1)) / (kv[i+k-1] - kv[i]);
00645     dvals[q+1] += -omega * dvals[q];
00646     dvals[q] *= omega;
00647   }
00648 }
00649 
00650 /***********************************************************************
00651  *                                                                     *
00652  * Calculate a point p on NurbSurface n at a specific u, v             *
00653  * using the tensor product.                                           *
00654  *                                                                     *
00655  * Note the valid parameter range for u and v is                       *
00656  * (kvU[orderU] <= u < kvU[numU), (kvV[orderV] <= v < kvV[numV])       *
00657  *                                                                     *
00658  ***********************************************************************/
00659 void G4NURBS::CalcPoint(G4double u, G4double v, G4Point3D &p,
00660                         G4Vector3D &utan, G4Vector3D &vtan) const
00661 {
00662 #define MAXORDER 50
00663   struct Point4
00664   {
00665     G4double x, y, z, w;
00666   };
00667 
00668   G4int i, j, ri, rj;
00669   G4int ubrkPoint, ufirst;
00670   G4double bu[MAXORDER], buprime[MAXORDER];
00671   G4int vbrkPoint, vfirst;
00672   G4double bv[MAXORDER], bvprime[MAXORDER];
00673   Point4 r, rutan, rvtan;
00674 
00675   r.x = 0.0; r.y = 0.0; r.z = 0.0; r.w = 0.0;
00676   rutan = r;   rvtan = r;
00677 
00678   G4int numU   = GetUnbrCtrlPts();
00679   G4int numV   = GetVnbrCtrlPts();
00680   G4int orderU = GetUorder();
00681   G4int orderV = GetVorder();
00682   
00683   /* Evaluate non-uniform basis functions (and derivatives) */
00684   
00685   ubrkPoint = FindBreakPoint(u, m[U].pKnots, numU-1, orderU);
00686   ufirst    = ubrkPoint - orderU + 1;
00687   BasisFunctions  (u, ubrkPoint, m[U].pKnots, orderU, bu);
00688   BasisDerivatives(u, ubrkPoint, m[U].pKnots, orderU, buprime);
00689 
00690   vbrkPoint = FindBreakPoint(v, m[V].pKnots, numV-1, orderV);
00691   vfirst    = vbrkPoint - orderV + 1;
00692   BasisFunctions  (v, vbrkPoint, m[V].pKnots, orderV, bv);
00693   BasisDerivatives(v, vbrkPoint, m[V].pKnots, orderV, bvprime);
00694 
00695   /* Weight control points against the basis functions */
00696 
00697   t_doubleCtrlPt *cpoint;
00698   Point4 cp;
00699   G4double tmp;
00700 
00701   for (i=0; i<orderV; i++)
00702   {
00703     for (j=0; j<orderU; j++)
00704     {
00705       ri = orderV - 1 - i;
00706       rj = orderU - 1 - j;
00707 
00708       tmp = bu[rj] * bv[ri];
00709 
00710       cpoint = GetdoubleCtrlPt(j+ufirst, i+vfirst);
00711       cp.x = *cpoint[G4NURBS::X];
00712       cp.y = *cpoint[G4NURBS::Y];
00713       cp.z = *cpoint[G4NURBS::Z];
00714       cp.w = *cpoint[G4NURBS::W];
00715       delete [] cpoint;
00716 
00717       r.x += cp.x * tmp;
00718       r.y += cp.y * tmp;
00719       r.z += cp.z * tmp;
00720       r.w += cp.w * tmp;
00721       
00722       tmp = buprime[rj] * bv[ri];
00723       rutan.x += cp.x * tmp;
00724       rutan.y += cp.y * tmp;
00725       rutan.z += cp.z * tmp;
00726       rutan.w += cp.w * tmp;
00727 
00728       tmp = bu[rj] * bvprime[ri];
00729       rvtan.x += cp.x * tmp;
00730       rvtan.y += cp.y * tmp;
00731       rvtan.z += cp.z * tmp;
00732       rvtan.w += cp.w * tmp;
00733     }
00734   }
00735 
00736   /* Project tangents, using the quotient rule for differentiation */
00737   
00738   G4double wsqrdiv = 1.0 / (r.w * r.w);
00739 
00740   utan.setX((r.w * rutan.x - rutan.w * r.x) * wsqrdiv);
00741   utan.setY((r.w * rutan.y - rutan.w * r.y) * wsqrdiv);
00742   utan.setZ((r.w * rutan.z - rutan.w * r.z) * wsqrdiv);
00743 
00744   vtan.setX((r.w * rvtan.x - rvtan.w * r.x) * wsqrdiv);
00745   vtan.setY((r.w * rvtan.y - rvtan.w * r.y) * wsqrdiv);
00746   vtan.setZ((r.w * rvtan.z - rvtan.w * r.z) * wsqrdiv);
00747 
00748   p.setX(r.x / r.w);
00749   p.setY(r.y / r.w);
00750   p.setZ(r.z / r.w);
00751 }

Generated on Mon May 27 17:49:07 2013 for Geant4 by  doxygen 1.4.7