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 #include "G4BSplineSurface.hh"
00037 #include "G4BezierSurface.hh"
00038 #include "G4ControlPoints.hh"
00039 #include "G4BoundingBox3D.hh"
00040
00041 G4BSplineSurface::G4BSplineSurface()
00042 : ord(0), k_index(0), param(0.), Rational(0)
00043 {
00044 distance = kInfinity;
00045 dir=ROW;
00046 first_hit = Hit = (G4UVHit*)0;
00047 order[0] = 0; order[1] = 0;
00048 ctl_points = (G4ControlPoints*)0;
00049 u_knots = v_knots = tmp_knots = (G4KnotVector*)0;
00050 }
00051
00052
00053 G4BSplineSurface::G4BSplineSurface(const char*, G4Ray&)
00054 : dir(0), ord(0), k_index(0), param(0.), Rational(0)
00055 {
00056 distance = kInfinity;
00057 first_hit = Hit = (G4UVHit*)0;
00058 order[0] = 0; order[1] = 0;
00059 ctl_points = (G4ControlPoints*)0;
00060 u_knots = v_knots = tmp_knots = (G4KnotVector*)0;
00061 }
00062
00063
00064 G4BSplineSurface::G4BSplineSurface(G4int u, G4int v, G4KnotVector& u_kv,
00065 G4KnotVector& v_kv, G4ControlPoints& cp)
00066 : dir(0), ord(0), k_index(0), param(0.), Rational(0)
00067 {
00068 first_hit = Hit = (G4UVHit*)0;
00069
00070 order[0] = u+1; order[1] = v+1;
00071
00072 u_knots = new G4KnotVector(u_kv);
00073 v_knots = new G4KnotVector(v_kv);
00074 tmp_knots = (G4KnotVector*)0;
00075
00076 ctl_points = new G4ControlPoints(cp);
00077 }
00078
00079
00080 G4BSplineSurface::~G4BSplineSurface()
00081 {
00082 delete u_knots;
00083 delete v_knots;
00084 delete ctl_points;
00085 G4UVHit* temphit=Hit;
00086 Hit = first_hit;
00087 while(Hit!=(G4UVHit*)0)
00088 {
00089 Hit=Hit->GetNext();
00090 delete temphit;
00091 temphit=Hit;
00092 }
00093
00094
00095 }
00096
00097
00098 G4int G4BSplineSurface::Intersect(const G4Ray& rayref)
00099 {
00100 Intersected = 1;
00101 FindIntersections(rayref);
00102 G4BezierSurface *bez_ptr;
00103 bezier_list.MoveToFirst();
00104 distance = kInfinity;
00105
00106 while( bezier_list.GetSurface() != (G4Surface*)0)
00107 {
00108 bez_ptr = (G4BezierSurface*)bezier_list.GetSurface();
00109
00110 if(bez_ptr->IsActive())
00111 {
00112 if(distance > bez_ptr->GetDistance())
00113 {
00114
00115 closest_hit = bez_ptr->AveragePoint();
00116 distance = bez_ptr->GetDistance();
00117 }
00118 else
00119 {
00120
00121 bez_ptr->SetActive(0);
00122
00123
00124
00125 }
00126 }
00127
00128 bezier_list.Step();
00129 }
00130
00131 bezier_list.MoveToFirst();
00132
00133 if(bezier_list.GetSize())
00134 return 1;
00135 else
00136 {
00137 active=0;
00138 return 0;
00139 }
00140 }
00141
00142
00143 G4Point3D G4BSplineSurface::FinalIntersection()
00144 {
00145
00146 G4BezierSurface* bez_ptr;
00147 while ( bezier_list.GetSize() > 0 &&
00148 bezier_list.GetSurface() != (G4Surface*)0)
00149 {
00150 bez_ptr = (G4BezierSurface*)bezier_list.GetSurface();
00151 G4int tmp = 0;
00152
00153
00154
00155
00156 tmp = bez_ptr->BIntersect( bezier_list);
00157
00158 if(!tmp)
00159 {
00160 bezier_list.RemoveSurface(bez_ptr);
00161 if(bezier_list.GetSurface() != (G4Surface*)0)
00162 bezier_list.GetSurface()->SetActive(1);
00163 }
00164 else
00165 if(tmp==1)
00166 {
00167 active=1;
00168
00169 AddHit(bez_ptr->GetU(), bez_ptr->GetV());
00170
00171
00172 bezier_list.EmptyList();
00173 }
00174 else
00175 if(tmp==2)
00176 {
00177
00178
00179
00180
00181
00182
00183 bezier_list.MoveToFirst();
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193 if(bezier_list.GetSurface() != bezier_list.GetLastSurface())
00194 while (bezier_list.GetNext() != bezier_list.GetLastSurface())
00195 bezier_list.Step();
00196
00197 G4BezierSurface* tmpbz = (G4BezierSurface*) bezier_list.GetSurface();
00198 tmpbz->CalcBBox();
00199
00200
00201
00202 G4int result=0;
00203 if(tmpbz->bbox->GetTestResult())
00204 {
00205
00206 while(!result)
00207 result = tmpbz->ClipBothDirs();
00208 }
00209 else
00210 {
00211 bezier_list.RemoveSurface(tmpbz);
00212 }
00213
00214 tmpbz = (G4BezierSurface*) bezier_list.GetLastSurface();
00215 tmpbz->CalcBBox();
00216
00217
00218
00219 if(tmpbz->bbox->GetTestResult())
00220 {
00221 result = 0;
00222 while(!result)
00223 result = tmpbz->ClipBothDirs();
00224 }
00225 else
00226 {
00227 bezier_list.RemoveSurface(tmpbz);
00228 }
00229
00230 bezier_list.RemoveSurface(bez_ptr);
00231 bezier_list.MoveToFirst();
00232 }
00233
00234 bezier_list.Step();
00235 }
00236
00237 Hit = first_hit;
00238 G4Point3D result;
00239 if(Hit == (G4UVHit*)0)
00240 active = 0;
00241 else
00242 {
00243 while(Hit != (G4UVHit*)0)
00244 {
00245
00246
00247
00248 result = BSEvaluate();
00249
00250 Hit = Hit->GetNext();
00251 }
00252
00253 Hit = first_hit;
00254 }
00255
00256 return result;
00257 }
00258
00259
00260 void G4BSplineSurface::CalcBBox()
00261 {
00262
00263
00264
00265
00266
00267
00268 G4Point3D box_min = G4Point3D( PINFINITY);
00269 G4Point3D box_max = G4Point3D(-PINFINITY);
00270
00271
00272
00273
00274 for(register int a = ctl_points->GetRows()-1; a>=0;a--)
00275 for(register int b = ctl_points->GetCols()-1; b>=0;b--)
00276 {
00277 G4Point3D tmp = ctl_points->Get3D(a,b);
00278 if((box_min.x()) > (tmp.x())) box_min.setX(tmp.x());
00279 if((box_min.y()) > (tmp.y())) box_min.setY(tmp.y());
00280 if((box_min.z()) > (tmp.z())) box_min.setZ(tmp.z());
00281 if((box_max.x()) < (tmp.x())) box_max.setX(tmp.x());
00282 if((box_max.y()) < (tmp.y())) box_max.setY(tmp.y());
00283 if((box_max.z()) < (tmp.z())) box_max.setZ(tmp.z());
00284 }
00285 bbox = new G4BoundingBox3D( box_min, box_max);
00286 }
00287
00288
00289 G4ProjectedSurface* G4BSplineSurface::CopyToProjectedSurface
00290 (const G4Ray& rayref)
00291 {
00292 G4ProjectedSurface* proj_srf = new G4ProjectedSurface() ;
00293 proj_srf->PutOrder(0,GetOrder(0));
00294 proj_srf->PutOrder(1,GetOrder(1));
00295 proj_srf->dir = dir;
00296
00297 proj_srf->u_knots = new G4KnotVector(*u_knots);
00298 proj_srf->v_knots = new G4KnotVector(*v_knots);
00299 proj_srf->ctl_points = new G4ControlPoints
00300 (2, ctl_points->GetRows(), ctl_points->GetCols());
00301
00302 const G4Plane& plane1 = rayref.GetPlane(1);
00303 const G4Plane& plane2 = rayref.GetPlane(2);
00304 ProjectNURBSurfaceTo2D(plane1, plane2, proj_srf);
00305
00306 return proj_srf;
00307 }
00308
00309
00310 void G4BSplineSurface::FindIntersections(const G4Ray& rayref)
00311 {
00312
00313 G4ProjectedSurface* proj_srf = CopyToProjectedSurface(rayref);
00314
00315
00316 projected_list.AddSurface(proj_srf);
00317
00318
00319 while(projected_list.GetSize() > 0)
00320 {
00321
00322 proj_srf = (G4ProjectedSurface*)projected_list.GetSurface();
00323
00324
00325 proj_srf->CalcBBox();
00326
00327
00328
00329
00330 if(proj_srf->bbox->GetTestResult())
00331
00332 proj_srf->ConvertToBezier(projected_list, bezier_list);
00333
00334
00335 projected_list.RemoveSurface(proj_srf);
00336 }
00337
00338
00339 G4BezierSurface* bez_ptr;
00340 distance = kInfinity;
00341
00342 while(bezier_list.GetSurface() != (G4Surface*)0)
00343 {
00344 bez_ptr = (G4BezierSurface*)bezier_list.GetSurface();
00345
00346
00347 AddHit(bez_ptr->UAverage(), bez_ptr->VAverage());
00348
00349
00350
00351
00352
00353
00354 bez_ptr->SetAveragePoint(BSEvaluate());
00355
00356
00357 bez_ptr->CalcDistance(rayref.GetStart());
00358
00359
00360 if(bez_ptr->GetDistance() < distance) distance = bez_ptr->GetDistance();
00361
00362
00363 if (first_hit == Hit) first_hit = (G4UVHit*)0;
00364 delete Hit;
00365 Hit = (G4UVHit*)0;
00366
00367
00368 bezier_list.Step();
00369 }
00370
00371 bezier_list.MoveToFirst();
00372 if(bezier_list.GetSize() == 0)
00373 {
00374 active=0;
00375 return;
00376 }
00377
00378
00379 const G4Point3D& Pt = rayref.GetStart();
00380 const G4Vector3D& Dir = rayref.GetDir();
00381 G4Point3D TestPoint = G4Point3D( (0.00001*Dir) + Pt );
00382 G4BezierSurface* Bsrf = (G4BezierSurface*)bezier_list.GetSurface(0);
00383
00384 G4Point3D AveragePoint = Bsrf->AveragePoint();
00385 G4double TestDistance = TestPoint.distance2(AveragePoint);
00386
00387 if(TestDistance > distance)
00388
00389 active=0;
00390 }
00391
00392
00393 void G4BSplineSurface::AddHit(G4double u, G4double v)
00394 {
00395 if(Hit == (G4UVHit*)0)
00396 {
00397 first_hit = new G4UVHit(u,v);
00398 first_hit->SetNext(0);
00399 Hit = first_hit;
00400 }
00401 else
00402 {
00403 Hit->SetNext(new G4UVHit(u,v));
00404 Hit = Hit->GetNext();
00405 Hit->SetNext(0);
00406 }
00407 }
00408
00409
00410 void G4BSplineSurface::ProjectNURBSurfaceTo2D
00411 (const G4Plane& plane1, const G4Plane& plane2,
00412 register G4ProjectedSurface* proj_srf)
00413 {
00414
00415
00416
00417
00418
00419
00420 G4PointRat tmp = ctl_points->GetRat(0,0);
00421 G4int rational = tmp.GetType();
00422 G4Point3D psrfcoords;
00423 register int rows = ctl_points->GetRows();
00424 register int cols = ctl_points->GetCols();
00425
00426 for (register int i=0; i< rows; i++)
00427 for(register int j=0; j < cols;j++)
00428 {
00429 if ( rational==4 )
00430 {
00431 G4PointRat& srfcoords = ctl_points->GetRat(i, j);
00432
00433
00434
00435
00436
00437
00438 psrfcoords.setX(( srfcoords.x() * plane1.a
00439 +srfcoords.y() * plane1.b
00440 +srfcoords.z() * plane1.c
00441 -srfcoords.w() * plane1.d));
00442 psrfcoords.setY(( srfcoords.x() * plane2.a
00443 +srfcoords.y() * plane2.b
00444 +srfcoords.z() * plane2.c
00445 -srfcoords.w() * plane2.d));
00446
00447 proj_srf->ctl_points->put(i,j,psrfcoords);
00448 }
00449 else
00450 {
00451 G4Point3D srfcoords = ctl_points->Get3D(i, j);
00452
00453 psrfcoords.setX(( srfcoords.x() * plane1.a
00454 +srfcoords.y() * plane1.b
00455 +srfcoords.z() * plane1.c
00456 - plane1.d));
00457
00458 psrfcoords.setY(( srfcoords.x() * plane2.a
00459 +srfcoords.y() * plane2.b
00460 +srfcoords.z() * plane2.c
00461 - plane2.d));
00462
00463 proj_srf->ctl_points->put(i,j,psrfcoords);
00464 }
00465 }
00466 }
00467
00468
00469
00470
00471
00472 G4PointRat& G4BSplineSurface::InternalEvalCrv(G4int i, G4ControlPoints* crv)
00473 {
00474 if ( ord <= 1 )
00475 return crv->GetRat(i, k_index);
00476
00477 register int j = k_index;
00478
00479 while ( j > (k_index - ord + 1))
00480 {
00481 register G4double k1, k2;
00482
00483 k1 = tmp_knots->GetKnot((j + ord - 1));
00484 k2 = tmp_knots->GetKnot(j);
00485
00486 if ((std::abs(k1 - k2)) > kCarTolerance )
00487 {
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504 }
00505
00506 j--;
00507 }
00508
00509 ord = ord-1;
00510 return InternalEvalCrv(0, crv);
00511 }
00512
00513
00514 G4Point3D G4BSplineSurface::BSEvaluate()
00515 {
00516 register int i;
00517 register int row_size = ctl_points->GetRows();
00518 register G4ControlPoints *diff_curve;
00519 register G4ControlPoints* curves;
00520 G4Point3D result;
00521
00522
00523
00524
00525
00526 G4PointRat* tmp = &ctl_points->GetRat(0,0);
00527
00528 register int point_type = tmp->GetType();
00529 diff_curve = new G4ControlPoints(point_type, row_size, 1);
00530 k_index = u_knots->GetKnotIndex(Hit->GetU(), GetOrder(ROW) );
00531
00532 ord = GetOrder(ROW);
00533 if(k_index==-1)
00534 {
00535 delete diff_curve;
00536 active = 0;
00537 return result;
00538 }
00539
00540 curves=new G4ControlPoints(*ctl_points);
00541 tmp_knots = u_knots;
00542 param = Hit->GetU();
00543
00544 if(point_type == 4)
00545 {
00546 for ( i = 0; i < row_size; i++)
00547 {
00548 ord = GetOrder(ROW);
00549 G4PointRat rtr_pt = InternalEvalCrv(i, curves);
00550 diff_curve->put(0,i,rtr_pt);
00551 }
00552
00553 k_index = v_knots->GetKnotIndex( Hit->GetV(), GetOrder(COL) );
00554 if(k_index==-1)
00555 {
00556 delete diff_curve;
00557 delete curves;
00558 active = 0;
00559 return result;
00560 }
00561
00562 ord = GetOrder(COL);
00563 tmp_knots = v_knots;
00564 param = Hit->GetV();
00565
00566
00567
00568 G4PointRat rat_result(InternalEvalCrv(0, diff_curve));
00569
00570
00571
00572
00573 result.setX(rat_result.x()/rat_result.w());
00574 result.setY(rat_result.y()/rat_result.w());
00575 result.setZ(rat_result.z()/rat_result.w());
00576 }
00577 else
00578 if(point_type == 3)
00579 {
00580 for ( i = 0; i < row_size; i++)
00581 {
00582 ord = GetOrder(ROW);
00583
00584 G4Point3D rtr_pt = (InternalEvalCrv(i, curves)).pt();
00585 diff_curve->put(0,i,rtr_pt);
00586 }
00587
00588 k_index = v_knots->GetKnotIndex( Hit->GetV(), GetOrder(COL) );
00589 if(k_index==-1)
00590 {
00591 delete diff_curve;
00592 delete curves;
00593 active = 0;
00594 return result;
00595 }
00596
00597 ord = GetOrder(COL);
00598 tmp_knots = v_knots;
00599 param = Hit->GetV();
00600
00601
00602 result = (InternalEvalCrv(0, diff_curve)).pt();
00603 }
00604
00605 delete diff_curve;
00606 delete curves;
00607 closest_hit = result;
00608 return result;
00609 }
00610
00611
00612 G4Point3D G4BSplineSurface::Evaluation(const G4Ray&)
00613 {
00614
00615 G4UVHit* temphit=Hit;
00616 while(Hit!=(G4UVHit*)0)
00617 {
00618 Hit=Hit->GetNext();
00619 delete temphit;
00620 temphit=Hit;
00621 }
00622
00623
00624
00625
00626 closest_hit = FinalIntersection();
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637 return closest_hit;
00638 }
00639
00640
00641 G4double G4BSplineSurface::ClosestDistanceToPoint(const G4Point3D& Pt)
00642 {
00643 G4double PointDistance=0;
00644 PointDistance = ctl_points->ClosestDistanceToPoint(Pt);
00645 return PointDistance;
00646 }