00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include "G4NURBS.hh"
00038
00039
00040
00041
00043
00045
00046 std::ostream & operator << (std::ostream & inout_outStream,
00047 const G4NURBS & in_kNurb)
00048 {
00049 inout_outStream
00050
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
00057 for (G4NURBS::t_direction dir = G4NURBS::U; dir < G4NURBS::NofD;
00058 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 {
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 }
00075 }
00076
00077
00078
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 {
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 }
00099
00100 inout_outStream << "\n# That's all!"
00101 << G4endl;
00102 return inout_outStream;
00103 }
00104
00105
00106
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
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
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
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 {
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
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
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 }
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
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 }
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
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
00439
00440
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 {
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
00458 }
00459 }
00460 }
00461
00462
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
00481 mpCtrlPts = new t_CtrlPt [mtotnbrCtrlPts];
00482
00483
00484
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;
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
00505 }
00506 }
00507
00508 G4NURBS::G4NURBS(const G4NURBS & in_krNurb)
00509 : G4Visible(in_krNurb)
00510 {
00511
00512
00513
00514 mtotnbrCtrlPts = in_krNurb.mtotnbrCtrlPts;
00515
00516
00517
00518
00519
00520
00521 t_direction dir;
00522 G4int dummy;
00523 for (dummy=0; (dummy?(dir=V):(dir=U)),(dummy < NofD); dummy++)
00524 {
00525
00526 m[dir] = in_krNurb.m[dir];
00527
00528
00529
00530 m[dir].pKnots = new G4double [m[dir].nbrKnots];
00531
00532 memcpy(m[dir].pKnots, in_krNurb.m[dir].pKnots,
00533 m[dir].nbrKnots * sizeof(G4double));
00534 }
00535
00536
00537
00538 mpCtrlPts = new t_CtrlPt [mtotnbrCtrlPts];
00539 memcpy(mpCtrlPts, in_krNurb.mpCtrlPts, mtotnbrCtrlPts*sizeof(t_CtrlPt));
00540
00541
00542
00543 G4cerr << "\nWARNING: G4NURBS::G4NURBS(const G4NURBS &) used" << G4endl;
00544 }
00545
00546 G4NURBS::~G4NURBS()
00547 {
00548
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;
00555 m[dir].pKnots = 0;
00556 }
00557
00558 if (mpCtrlPts)
00559 delete [] mpCtrlPts;
00560 mpCtrlPts = 0;
00561 }
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
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;
00578 i = np + k;
00579 while ((u < kv[i]) && (i > 0)) i--;
00580 return(i);
00581 }
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
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
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;
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
00653
00654
00655
00656
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
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
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
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 }