[831] | 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 | // |
---|
[850] | 27 | // $Id: G4BSplineSurface.cc,v 1.15 2008/03/13 14:18:57 gcosmo Exp $ |
---|
[1228] | 28 | // GEANT4 tag $Name: geant4-09-03 $ |
---|
[831] | 29 | // |
---|
| 30 | // ---------------------------------------------------------------------- |
---|
| 31 | // GEANT 4 class source file |
---|
| 32 | // |
---|
| 33 | // G4BSplineSurface.cc |
---|
| 34 | // |
---|
| 35 | // ---------------------------------------------------------------------- |
---|
| 36 | |
---|
| 37 | #include "G4BSplineSurface.hh" |
---|
| 38 | #include "G4BezierSurface.hh" |
---|
| 39 | #include "G4ControlPoints.hh" |
---|
| 40 | #include "G4BoundingBox3D.hh" |
---|
| 41 | |
---|
| 42 | G4BSplineSurface::G4BSplineSurface() |
---|
| 43 | { |
---|
| 44 | distance = kInfinity; |
---|
| 45 | dir=ROW; |
---|
| 46 | first_hit = Hit = (G4UVHit*)0; |
---|
| 47 | ctl_points = (G4ControlPoints*)0; |
---|
| 48 | u_knots = v_knots = tmp_knots = (G4KnotVector*)0; |
---|
| 49 | } |
---|
| 50 | |
---|
| 51 | |
---|
| 52 | G4BSplineSurface::G4BSplineSurface(const char*, G4Ray&) |
---|
| 53 | { |
---|
| 54 | distance = kInfinity; |
---|
| 55 | first_hit = Hit = (G4UVHit*)0; |
---|
| 56 | ctl_points = (G4ControlPoints*)0; |
---|
| 57 | u_knots = v_knots = tmp_knots = (G4KnotVector*)0; |
---|
| 58 | } |
---|
| 59 | |
---|
| 60 | |
---|
| 61 | G4BSplineSurface::G4BSplineSurface(G4int u, G4int v, G4KnotVector& u_kv, |
---|
| 62 | G4KnotVector& v_kv, G4ControlPoints& cp) |
---|
| 63 | { |
---|
| 64 | first_hit = Hit = (G4UVHit*)0; |
---|
| 65 | |
---|
| 66 | order[0] = u+1; |
---|
| 67 | order[1] = v+1; |
---|
| 68 | |
---|
| 69 | u_knots = new G4KnotVector(u_kv); |
---|
| 70 | v_knots = new G4KnotVector(v_kv); |
---|
| 71 | tmp_knots = (G4KnotVector*)0; |
---|
| 72 | |
---|
| 73 | ctl_points = new G4ControlPoints(cp); |
---|
| 74 | } |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | G4BSplineSurface::~G4BSplineSurface() |
---|
| 78 | { |
---|
| 79 | delete u_knots; |
---|
| 80 | delete v_knots; |
---|
| 81 | delete ctl_points; |
---|
| 82 | G4UVHit* temphit=Hit; |
---|
| 83 | Hit = first_hit; |
---|
| 84 | while(Hit!=(G4UVHit*)0) |
---|
| 85 | { |
---|
| 86 | Hit=Hit->GetNext(); |
---|
| 87 | delete temphit; |
---|
| 88 | temphit=Hit; |
---|
| 89 | } |
---|
| 90 | // delete temphit;// remove last |
---|
| 91 | |
---|
| 92 | } |
---|
| 93 | |
---|
| 94 | |
---|
| 95 | G4int G4BSplineSurface::Intersect(const G4Ray& rayref) |
---|
| 96 | { |
---|
| 97 | Intersected = 1; |
---|
| 98 | FindIntersections(rayref); |
---|
| 99 | G4BezierSurface *bez_ptr; |
---|
| 100 | bezier_list.MoveToFirst(); |
---|
| 101 | distance = kInfinity; |
---|
| 102 | |
---|
| 103 | while( bezier_list.GetSurface() != (G4Surface*)0) |
---|
| 104 | { |
---|
| 105 | bez_ptr = (G4BezierSurface*)bezier_list.GetSurface(); |
---|
| 106 | |
---|
| 107 | if(bez_ptr->IsActive()) |
---|
| 108 | { |
---|
| 109 | if(distance > bez_ptr->GetDistance()) |
---|
| 110 | { |
---|
| 111 | // Put data from closest bezier to b-spline data struct |
---|
| 112 | closest_hit = bez_ptr->AveragePoint(); |
---|
| 113 | distance = bez_ptr->GetDistance(); |
---|
| 114 | } |
---|
| 115 | else |
---|
| 116 | { |
---|
| 117 | // Set other beziers as inactive |
---|
| 118 | bez_ptr->SetActive(0); |
---|
| 119 | |
---|
| 120 | // Remove beziers that are not closest |
---|
| 121 | // bezier_list.RemoveSurface(bez_ptr); |
---|
| 122 | } |
---|
| 123 | } |
---|
| 124 | |
---|
| 125 | bezier_list.Step(); |
---|
| 126 | } |
---|
| 127 | |
---|
| 128 | bezier_list.MoveToFirst(); |
---|
| 129 | |
---|
| 130 | if(bezier_list.GetSize()) |
---|
| 131 | return 1; |
---|
| 132 | else |
---|
| 133 | { |
---|
| 134 | active=0; |
---|
| 135 | return 0; |
---|
| 136 | } |
---|
| 137 | } |
---|
| 138 | |
---|
| 139 | |
---|
| 140 | G4Point3D G4BSplineSurface::FinalIntersection() |
---|
| 141 | { |
---|
| 142 | // Compute the real intersection point. |
---|
| 143 | G4BezierSurface* bez_ptr; |
---|
| 144 | while ( bezier_list.GetSize() > 0 && |
---|
| 145 | bezier_list.GetSurface() != (G4Surface*)0) |
---|
| 146 | { |
---|
| 147 | bez_ptr = (G4BezierSurface*)bezier_list.GetSurface(); |
---|
| 148 | int tmp = 0; |
---|
| 149 | |
---|
| 150 | // L. Broglia |
---|
| 151 | // Modify G4BezierSurface intersection function name |
---|
| 152 | // tmp = bez_ptr->Intersect( bezier_list); |
---|
| 153 | tmp = bez_ptr->BIntersect( bezier_list); |
---|
| 154 | |
---|
| 155 | if(!tmp) |
---|
| 156 | { |
---|
| 157 | bezier_list.RemoveSurface(bez_ptr); |
---|
| 158 | if(bezier_list.GetSurface() != (G4Surface*)0) |
---|
| 159 | bezier_list.GetSurface()->SetActive(1); |
---|
| 160 | } |
---|
| 161 | else |
---|
| 162 | if(tmp==1) |
---|
| 163 | { |
---|
| 164 | active=1; |
---|
| 165 | // Hit found |
---|
| 166 | AddHit(bez_ptr->GetU(), bez_ptr->GetV()); |
---|
| 167 | |
---|
| 168 | // Delete beziers |
---|
| 169 | bezier_list.EmptyList(); |
---|
| 170 | } |
---|
| 171 | else |
---|
| 172 | if(tmp==2) |
---|
| 173 | { |
---|
| 174 | // The bezier was split so the last |
---|
| 175 | // two surfaces in the List should |
---|
| 176 | // be bbox tested and if passed |
---|
| 177 | // clipped in both dirs. |
---|
| 178 | |
---|
| 179 | // Move to first |
---|
| 180 | bezier_list.MoveToFirst(); |
---|
| 181 | // Find the second last. |
---|
| 182 | // What!? Casting a G4Surface* to a G4SurfaceList* !?!?!? - GC |
---|
| 183 | // |
---|
| 184 | // if(bezier_list.index != bezier_list.last) |
---|
| 185 | // while ( ((G4SurfaceList*)bezier_list.index)->next != |
---|
| 186 | // bezier_list.last) bezier_list.Step(); |
---|
| 187 | // |
---|
| 188 | // Try the following instead (if that's the wished behavior)... |
---|
| 189 | // |
---|
| 190 | if(bezier_list.GetSurface() != bezier_list.GetLastSurface()) |
---|
| 191 | while (bezier_list.GetNext() != bezier_list.GetLastSurface()) |
---|
| 192 | bezier_list.Step(); |
---|
| 193 | |
---|
| 194 | G4BezierSurface* tmp = (G4BezierSurface*) bezier_list.GetSurface(); |
---|
| 195 | tmp->CalcBBox(); |
---|
| 196 | |
---|
| 197 | // L. Broglia tmp->bbox->Test(); |
---|
| 198 | |
---|
| 199 | int result=0; |
---|
| 200 | if(tmp->bbox->GetTestResult()) |
---|
| 201 | { |
---|
| 202 | // Clip |
---|
| 203 | while(!result) |
---|
| 204 | result = tmp->ClipBothDirs(); |
---|
| 205 | } |
---|
| 206 | else |
---|
| 207 | { |
---|
| 208 | bezier_list.RemoveSurface(tmp); |
---|
| 209 | } |
---|
| 210 | // Second surface |
---|
| 211 | tmp = (G4BezierSurface*) bezier_list.GetLastSurface(); |
---|
| 212 | tmp->CalcBBox(); |
---|
| 213 | |
---|
| 214 | // L. Broglia tmp->bbox->Test(); |
---|
| 215 | |
---|
| 216 | if(tmp->bbox->GetTestResult()) |
---|
| 217 | { |
---|
| 218 | result = 0; |
---|
| 219 | while(!result) |
---|
| 220 | result = tmp->ClipBothDirs(); |
---|
| 221 | } |
---|
| 222 | else |
---|
| 223 | { |
---|
| 224 | bezier_list.RemoveSurface(tmp); |
---|
| 225 | } |
---|
| 226 | |
---|
| 227 | bezier_list.RemoveSurface(bez_ptr); |
---|
| 228 | bezier_list.MoveToFirst(); |
---|
| 229 | } |
---|
| 230 | |
---|
| 231 | bezier_list.Step(); |
---|
| 232 | }//While.... |
---|
| 233 | |
---|
| 234 | Hit = first_hit; |
---|
| 235 | G4Point3D result; |
---|
| 236 | if(Hit == (G4UVHit*)0) |
---|
| 237 | active = 0; |
---|
| 238 | else |
---|
| 239 | { |
---|
| 240 | while(Hit != (G4UVHit*)0) |
---|
| 241 | { |
---|
| 242 | // L. Broglia |
---|
| 243 | // Modify function name |
---|
| 244 | // result = Evaluate(); |
---|
| 245 | result = BSEvaluate(); |
---|
| 246 | |
---|
| 247 | Hit = Hit->GetNext(); |
---|
| 248 | } |
---|
| 249 | |
---|
| 250 | Hit = first_hit; |
---|
| 251 | } |
---|
| 252 | |
---|
| 253 | return result; |
---|
| 254 | } |
---|
| 255 | |
---|
| 256 | |
---|
| 257 | void G4BSplineSurface::CalcBBox() |
---|
| 258 | { |
---|
| 259 | |
---|
| 260 | // Finds the bounds of the b-spline surface iow |
---|
| 261 | // calculates the bounds for a bounding box |
---|
| 262 | // to the surface. The bounding box is used |
---|
| 263 | // for a preliminary check of intersection. |
---|
| 264 | |
---|
| 265 | G4Point3D box_min = G4Point3D( PINFINITY); |
---|
| 266 | G4Point3D box_max = G4Point3D(-PINFINITY); |
---|
| 267 | |
---|
| 268 | // Loop to search the whole control point mesh |
---|
| 269 | // for the minimum and maximum values for x, y and z. |
---|
| 270 | |
---|
| 271 | for(register int a = ctl_points->GetRows()-1; a>=0;a--) |
---|
| 272 | for(register int b = ctl_points->GetCols()-1; b>=0;b--) |
---|
| 273 | { |
---|
| 274 | G4Point3D tmp = ctl_points->Get3D(a,b); |
---|
| 275 | if((box_min.x()) > (tmp.x())) box_min.setX(tmp.x()); |
---|
| 276 | if((box_min.y()) > (tmp.y())) box_min.setY(tmp.y()); |
---|
| 277 | if((box_min.z()) > (tmp.z())) box_min.setZ(tmp.z()); |
---|
| 278 | if((box_max.x()) < (tmp.x())) box_max.setX(tmp.x()); |
---|
| 279 | if((box_max.y()) < (tmp.y())) box_max.setY(tmp.y()); |
---|
| 280 | if((box_max.z()) < (tmp.z())) box_max.setZ(tmp.z()); |
---|
| 281 | } |
---|
| 282 | bbox = new G4BoundingBox3D( box_min, box_max); |
---|
| 283 | } |
---|
| 284 | |
---|
| 285 | |
---|
| 286 | G4ProjectedSurface* G4BSplineSurface::CopyToProjectedSurface |
---|
| 287 | (const G4Ray& rayref) |
---|
| 288 | { |
---|
| 289 | G4ProjectedSurface* proj_srf = new G4ProjectedSurface() ; |
---|
| 290 | proj_srf->PutOrder(0,GetOrder(0)); |
---|
| 291 | proj_srf->PutOrder(1,GetOrder(1)); |
---|
| 292 | proj_srf->dir = dir; |
---|
| 293 | |
---|
| 294 | proj_srf->u_knots = new G4KnotVector(*u_knots); |
---|
| 295 | proj_srf->v_knots = new G4KnotVector(*v_knots); |
---|
| 296 | proj_srf->ctl_points = new G4ControlPoints |
---|
| 297 | (2, ctl_points->GetRows(), ctl_points->GetCols()); |
---|
| 298 | |
---|
| 299 | const G4Plane& plane1 = rayref.GetPlane(1); |
---|
| 300 | const G4Plane& plane2 = rayref.GetPlane(2); |
---|
| 301 | ProjectNURBSurfaceTo2D(plane1, plane2, proj_srf); |
---|
| 302 | |
---|
| 303 | return proj_srf; |
---|
| 304 | } |
---|
| 305 | |
---|
| 306 | |
---|
| 307 | void G4BSplineSurface::FindIntersections(const G4Ray& rayref) |
---|
| 308 | { |
---|
| 309 | // Do the projection to 2D |
---|
| 310 | G4ProjectedSurface* proj_srf = CopyToProjectedSurface(rayref); |
---|
| 311 | |
---|
| 312 | // Put surface in projected List |
---|
| 313 | projected_list.AddSurface(proj_srf); |
---|
| 314 | |
---|
| 315 | // Loop through List of projected surfaces |
---|
| 316 | while(projected_list.GetSize() > 0) |
---|
| 317 | { |
---|
| 318 | // Get first in List |
---|
| 319 | proj_srf = (G4ProjectedSurface*)projected_list.GetSurface(); |
---|
| 320 | |
---|
| 321 | // Create the bounding box for the projected surface. |
---|
| 322 | proj_srf->CalcBBox(); |
---|
| 323 | |
---|
| 324 | // L. Broglia proj_srf->bbox->Test(); |
---|
| 325 | |
---|
| 326 | // Check bbox test result is ok |
---|
| 327 | if(proj_srf->bbox->GetTestResult()) |
---|
| 328 | // Convert the projected surface to a bezier. Split if necessary. |
---|
| 329 | proj_srf->ConvertToBezier(projected_list, bezier_list); |
---|
| 330 | |
---|
| 331 | // Remove projected surface |
---|
| 332 | projected_list.RemoveSurface(proj_srf); |
---|
| 333 | } |
---|
| 334 | |
---|
| 335 | // Loop through the bezier List |
---|
| 336 | G4BezierSurface* bez_ptr; |
---|
| 337 | distance = kInfinity; |
---|
| 338 | |
---|
| 339 | while(bezier_list.GetSurface() != (G4Surface*)0) |
---|
| 340 | { |
---|
| 341 | bez_ptr = (G4BezierSurface*)bezier_list.GetSurface(); |
---|
| 342 | |
---|
| 343 | // Add a temporary Hit |
---|
| 344 | AddHit(bez_ptr->UAverage(), bez_ptr->VAverage()); |
---|
| 345 | |
---|
| 346 | // Evaluate Hit |
---|
| 347 | |
---|
| 348 | // L. Broglia |
---|
| 349 | // Modify function name |
---|
| 350 | // bez_ptr->SetAveragePoint(Evaluate()); |
---|
| 351 | bez_ptr->SetAveragePoint(BSEvaluate()); |
---|
| 352 | |
---|
| 353 | // Calculate distance to ray origin |
---|
| 354 | bez_ptr->CalcDistance(rayref.GetStart()); |
---|
| 355 | |
---|
| 356 | // Put closest to b_splines distance value |
---|
| 357 | if(bez_ptr->GetDistance() < distance) distance = bez_ptr->GetDistance(); |
---|
| 358 | |
---|
| 359 | // Remove the temporary Hit |
---|
| 360 | if (first_hit == Hit) first_hit = (G4UVHit*)0; |
---|
| 361 | delete Hit; |
---|
| 362 | Hit = (G4UVHit*)0; |
---|
| 363 | |
---|
| 364 | // Move to next in the List |
---|
| 365 | bezier_list.Step(); |
---|
| 366 | } |
---|
| 367 | |
---|
| 368 | bezier_list.MoveToFirst(); |
---|
| 369 | if(bezier_list.GetSize() == 0) |
---|
| 370 | { |
---|
| 371 | active=0; |
---|
| 372 | return; |
---|
| 373 | } |
---|
| 374 | |
---|
| 375 | // Check that approx Hit is in direction of ray |
---|
| 376 | const G4Point3D& Pt = rayref.GetStart(); |
---|
| 377 | const G4Vector3D& Dir = rayref.GetDir(); |
---|
| 378 | G4Point3D TestPoint = G4Point3D( (0.00001*Dir) + Pt ); |
---|
| 379 | G4BezierSurface* Bsrf = (G4BezierSurface*)bezier_list.GetSurface(0); |
---|
| 380 | |
---|
| 381 | G4Point3D AveragePoint = Bsrf->AveragePoint(); |
---|
| 382 | G4double TestDistance = TestPoint.distance2(AveragePoint); |
---|
| 383 | |
---|
| 384 | if(TestDistance > distance) |
---|
| 385 | // Hit behind ray starting point, no intersection. |
---|
| 386 | active=0; |
---|
| 387 | } |
---|
| 388 | |
---|
| 389 | |
---|
| 390 | void G4BSplineSurface::AddHit(G4double u, G4double v) |
---|
| 391 | { |
---|
| 392 | if(Hit == (G4UVHit*)0) |
---|
| 393 | { |
---|
| 394 | first_hit = new G4UVHit(u,v); |
---|
| 395 | first_hit->SetNext(0); |
---|
| 396 | Hit = first_hit; |
---|
| 397 | } |
---|
| 398 | else |
---|
| 399 | { |
---|
| 400 | Hit->SetNext(new G4UVHit(u,v)); |
---|
| 401 | Hit = Hit->GetNext(); |
---|
| 402 | Hit->SetNext(0); |
---|
| 403 | } |
---|
| 404 | } |
---|
| 405 | |
---|
| 406 | |
---|
| 407 | void G4BSplineSurface::ProjectNURBSurfaceTo2D |
---|
| 408 | (const G4Plane& plane1, const G4Plane& plane2, |
---|
| 409 | register G4ProjectedSurface* proj_srf) |
---|
| 410 | { |
---|
| 411 | // Projects the nurb surface so that the z-axis = ray. |
---|
| 412 | |
---|
| 413 | /* L. Broglia |
---|
| 414 | G4Point* tmp = (G4Point*)&ctl_points->get(0,0); |
---|
| 415 | */ |
---|
| 416 | |
---|
| 417 | G4PointRat tmp = ctl_points->GetRat(0,0); |
---|
| 418 | int rational = tmp.GetType();// Get the type of control point |
---|
| 419 | G4Point3D psrfcoords; |
---|
| 420 | register int rows = ctl_points->GetRows(); |
---|
| 421 | register int cols = ctl_points->GetCols(); |
---|
| 422 | |
---|
| 423 | for (register int i=0; i< rows; i++) |
---|
| 424 | for(register int j=0; j < cols;j++) |
---|
| 425 | { |
---|
| 426 | if ( rational==4 ) // 4 coordinates |
---|
| 427 | { |
---|
| 428 | G4PointRat& srfcoords = ctl_points->GetRat(i, j); |
---|
| 429 | |
---|
| 430 | // L. Broglia |
---|
| 431 | // Changes for new G4PointRat |
---|
| 432 | |
---|
| 433 | // Calculate the x- and y-coordinates for the new |
---|
| 434 | // 2-D surface. |
---|
| 435 | psrfcoords.setX(( srfcoords.x() * plane1.a |
---|
| 436 | +srfcoords.y() * plane1.b |
---|
| 437 | +srfcoords.z() * plane1.c |
---|
| 438 | -srfcoords.w() * plane1.d)); |
---|
| 439 | psrfcoords.setY(( srfcoords.x() * plane2.a |
---|
| 440 | +srfcoords.y() * plane2.b |
---|
| 441 | +srfcoords.z() * plane2.c |
---|
| 442 | -srfcoords.w() * plane2.d)); |
---|
| 443 | |
---|
| 444 | proj_srf->ctl_points->put(i,j,psrfcoords); |
---|
| 445 | } |
---|
| 446 | else // 3 coordinates |
---|
| 447 | { |
---|
| 448 | G4Point3D srfcoords = ctl_points->Get3D(i, j); |
---|
| 449 | |
---|
| 450 | psrfcoords.setX(( srfcoords.x() * plane1.a |
---|
| 451 | +srfcoords.y() * plane1.b |
---|
| 452 | +srfcoords.z() * plane1.c |
---|
| 453 | - plane1.d)); |
---|
| 454 | |
---|
| 455 | psrfcoords.setY(( srfcoords.x() * plane2.a |
---|
| 456 | +srfcoords.y() * plane2.b |
---|
| 457 | +srfcoords.z() * plane2.c |
---|
| 458 | - plane2.d)); |
---|
| 459 | |
---|
| 460 | proj_srf->ctl_points->put(i,j,psrfcoords); |
---|
| 461 | } |
---|
| 462 | } |
---|
| 463 | } |
---|
| 464 | |
---|
| 465 | /* L. Broglia |
---|
| 466 | Changes for new G4PointRat |
---|
| 467 | G4Point& G4BSplineSurface::InternalEvalCrv(int i, G4ControlPoints* crv)*/ |
---|
| 468 | |
---|
| 469 | G4PointRat& G4BSplineSurface::InternalEvalCrv(int i, G4ControlPoints* crv) |
---|
| 470 | { |
---|
| 471 | if ( ord <= 1 ) |
---|
| 472 | return crv->GetRat(i, k_index); |
---|
| 473 | |
---|
| 474 | register int j = k_index; |
---|
| 475 | |
---|
| 476 | while ( j > (k_index - ord + 1)) |
---|
| 477 | { |
---|
| 478 | register G4double k1, k2; |
---|
| 479 | |
---|
| 480 | k1 = tmp_knots->GetKnot((j + ord - 1)); |
---|
| 481 | k2 = tmp_knots->GetKnot(j); |
---|
| 482 | |
---|
| 483 | if ((std::abs(k1 - k2)) > kCarTolerance ) |
---|
| 484 | { |
---|
| 485 | /* L. Broglia |
---|
| 486 | register G4PointRat* pts1 = &crv->get(i,j-1); |
---|
| 487 | register G4PointRat* pts2 = &crv->get(i,j ); |
---|
| 488 | if(pts1->GetType()==3) |
---|
| 489 | { |
---|
| 490 | crv->CalcValues(k1, param, *(G4Point3D*)pts1, k2, *(G4Point3D*)pts2); |
---|
| 491 | crv->put(0, j, *(G4Point3D*)pts2); |
---|
| 492 | } |
---|
| 493 | else |
---|
| 494 | { |
---|
| 495 | crv->CalcValues(k1, param, *(G4PointRat*)pts1, k2, *(G4PointRat*)pts2); |
---|
| 496 | crv->put(0, j, *(G4PointRat*)pts2); |
---|
| 497 | } |
---|
| 498 | register G4PointRat* pts1 = &crv->GetRat(i,j-1); |
---|
| 499 | register G4PointRat* pts2 = &crv->GetRat(i,j ); |
---|
| 500 | */ |
---|
| 501 | } |
---|
| 502 | |
---|
| 503 | j--; |
---|
| 504 | } |
---|
| 505 | |
---|
| 506 | ord = ord-1; |
---|
| 507 | return InternalEvalCrv(0, crv); // Recursion |
---|
| 508 | } |
---|
| 509 | |
---|
| 510 | |
---|
| 511 | G4Point3D G4BSplineSurface::BSEvaluate() |
---|
| 512 | { |
---|
| 513 | register int i; |
---|
| 514 | register int row_size = ctl_points->GetRows(); |
---|
| 515 | register G4ControlPoints *diff_curve; |
---|
| 516 | register G4ControlPoints* curves; |
---|
| 517 | G4Point3D result; |
---|
| 518 | |
---|
| 519 | /* L. Broglia |
---|
| 520 | G4Point* tmp = (G4Point*)&ctl_points->get(0,0); |
---|
| 521 | */ |
---|
| 522 | |
---|
| 523 | G4PointRat* tmp = &ctl_points->GetRat(0,0); |
---|
| 524 | |
---|
| 525 | register int point_type = tmp->GetType(); |
---|
| 526 | diff_curve = new G4ControlPoints(point_type, row_size, 1); |
---|
| 527 | k_index = u_knots->GetKnotIndex(Hit->GetU(), GetOrder(ROW) ); |
---|
| 528 | |
---|
| 529 | ord = GetOrder(ROW); |
---|
| 530 | if(k_index==-1) |
---|
| 531 | { |
---|
| 532 | delete diff_curve; |
---|
| 533 | active = 0; |
---|
| 534 | return result; |
---|
| 535 | } |
---|
| 536 | |
---|
| 537 | curves=new G4ControlPoints(*ctl_points); |
---|
| 538 | tmp_knots = u_knots; |
---|
| 539 | param = Hit->GetU(); |
---|
| 540 | |
---|
| 541 | if(point_type == 4) |
---|
| 542 | { |
---|
| 543 | for ( i = 0; i < row_size; i++) |
---|
| 544 | { |
---|
| 545 | ord = GetOrder(ROW); |
---|
| 546 | G4PointRat rtr_pt = InternalEvalCrv(i, curves); |
---|
| 547 | diff_curve->put(0,i,rtr_pt); |
---|
| 548 | } |
---|
| 549 | |
---|
| 550 | k_index = v_knots->GetKnotIndex( Hit->GetV(), GetOrder(COL) ); |
---|
| 551 | if(k_index==-1) |
---|
| 552 | { |
---|
| 553 | delete diff_curve; |
---|
| 554 | delete curves; |
---|
| 555 | active = 0; |
---|
| 556 | return result; |
---|
| 557 | } |
---|
| 558 | |
---|
| 559 | ord = GetOrder(COL); |
---|
| 560 | tmp_knots = v_knots; |
---|
| 561 | param = Hit->GetV(); |
---|
| 562 | |
---|
| 563 | // Evaluate the diff_curve... |
---|
| 564 | // G4PointRat rat_result = (G4PointRat&) InternalEvalCrv(0, diff_curve); |
---|
| 565 | G4PointRat rat_result(InternalEvalCrv(0, diff_curve)); |
---|
| 566 | |
---|
| 567 | // Calc the 3D values. |
---|
| 568 | // L. Broglia |
---|
| 569 | // Changes for new G4PointRat |
---|
| 570 | result.setX(rat_result.x()/rat_result.w()); |
---|
| 571 | result.setY(rat_result.y()/rat_result.w()); |
---|
| 572 | result.setZ(rat_result.z()/rat_result.w()); |
---|
| 573 | } |
---|
| 574 | else |
---|
| 575 | if(point_type == 3) |
---|
| 576 | { |
---|
| 577 | for ( i = 0; i < row_size; i++) |
---|
| 578 | { |
---|
| 579 | ord = GetOrder(ROW); |
---|
| 580 | // G4Point3D rtr_pt = (G4Point3D&) InternalEvalCrv(i, curves); |
---|
| 581 | G4Point3D rtr_pt = (InternalEvalCrv(i, curves)).pt(); |
---|
| 582 | diff_curve->put(0,i,rtr_pt); |
---|
| 583 | } |
---|
| 584 | |
---|
| 585 | k_index = v_knots->GetKnotIndex( Hit->GetV(), GetOrder(COL) ); |
---|
| 586 | if(k_index==-1) |
---|
| 587 | { |
---|
| 588 | delete diff_curve; |
---|
| 589 | delete curves; |
---|
| 590 | active = 0; |
---|
| 591 | return result; |
---|
| 592 | } |
---|
| 593 | |
---|
| 594 | ord = GetOrder(COL); |
---|
| 595 | tmp_knots = v_knots; |
---|
| 596 | param = Hit->GetV(); |
---|
| 597 | |
---|
| 598 | // Evaluate the diff_curve... |
---|
| 599 | result = (InternalEvalCrv(0, diff_curve)).pt(); |
---|
| 600 | } |
---|
| 601 | |
---|
| 602 | delete diff_curve; |
---|
| 603 | delete curves; |
---|
| 604 | closest_hit = result; |
---|
| 605 | return result; |
---|
| 606 | } |
---|
| 607 | |
---|
| 608 | |
---|
| 609 | G4Point3D G4BSplineSurface::Evaluation(const G4Ray&) |
---|
| 610 | { |
---|
| 611 | // Delete old UVhits |
---|
| 612 | G4UVHit* temphit=Hit; |
---|
| 613 | while(Hit!=(G4UVHit*)0) |
---|
| 614 | { |
---|
| 615 | Hit=Hit->GetNext(); |
---|
| 616 | delete temphit; |
---|
| 617 | temphit=Hit; |
---|
| 618 | } |
---|
| 619 | |
---|
| 620 | // delete temphit; |
---|
| 621 | |
---|
| 622 | // Get the real Hit point |
---|
| 623 | closest_hit = FinalIntersection(); |
---|
| 624 | |
---|
| 625 | // The following part (commented out) is old bullshit |
---|
| 626 | // Check that Hit is not in a void i.e. InnerBoundary. |
---|
| 627 | // for(int a=0; a<NumberOfInnerBoundaries;a++) |
---|
| 628 | // if(InnerBoundary[a]->Inside(closest_hit, rayref)) |
---|
| 629 | // { |
---|
| 630 | // Active(0); |
---|
| 631 | // Distance(kInfinity); |
---|
| 632 | // return closest_hit; |
---|
| 633 | // } |
---|
| 634 | return closest_hit; |
---|
| 635 | } |
---|
| 636 | |
---|
| 637 | |
---|
| 638 | G4double G4BSplineSurface::ClosestDistanceToPoint(const G4Point3D& Pt) |
---|
| 639 | { |
---|
| 640 | G4double PointDistance=0; |
---|
| 641 | PointDistance = ctl_points->ClosestDistanceToPoint(Pt); |
---|
| 642 | return PointDistance; |
---|
| 643 | } |
---|