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