G4BSplineSurface.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 // GEANT 4 class source file
00031 //
00032 // G4BSplineSurface.cc
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   // delete temphit;// remove last
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         // Put data from closest bezier to b-spline data struct
00115         closest_hit = bez_ptr->AveragePoint();
00116         distance = bez_ptr->GetDistance();
00117       }
00118       else
00119       {
00120         // Set other beziers as inactive
00121         bez_ptr->SetActive(0);
00122             
00123         // Remove beziers that are not closest
00124         //  bezier_list.RemoveSurface(bez_ptr);
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   // Compute the real intersection point.
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     // L. Broglia
00154     // Modify G4BezierSurface intersection function name 
00155     // tmp = bez_ptr->Intersect( bezier_list);
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         // Hit found
00169         AddHit(bez_ptr->GetU(), bez_ptr->GetV());
00170 
00171         // Delete beziers
00172         bezier_list.EmptyList();
00173       }
00174       else
00175         if(tmp==2)
00176         {       
00177           // The bezier was split so the last
00178           // two surfaces in the List should
00179           // be bbox tested and if passed
00180           // clipped in both dirs.
00181                 
00182           // Move to first
00183           bezier_list.MoveToFirst();
00184           // Find the second last.
00185 // What!?  Casting a G4Surface* to a G4SurfaceList*  !?!?!? - GC
00186 //
00187 //        if(bezier_list.index != bezier_list.last)
00188 //          while ( ((G4SurfaceList*)bezier_list.index)->next !=
00189 //                  bezier_list.last)  bezier_list.Step();
00190 //
00191 // Try the following instead (if that's the wished behavior)...
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 // L. Broglia           tmpbz->bbox->Test();
00201 
00202           G4int result=0;
00203           if(tmpbz->bbox->GetTestResult())
00204           {
00205             // Clip
00206             while(!result)
00207               result = tmpbz->ClipBothDirs();  
00208           }
00209           else
00210           {
00211             bezier_list.RemoveSurface(tmpbz);
00212           }
00213           // Second surface
00214           tmpbz = (G4BezierSurface*) bezier_list.GetLastSurface();
00215           tmpbz->CalcBBox();
00216 
00217 // L. Broglia           tmpbz->bbox->Test();
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   }//While....
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       // L. Broglia
00246       // Modify function name
00247       // result = Evaluate();
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   // Finds the bounds of the b-spline surface iow
00264   // calculates the bounds for a bounding box
00265   // to the surface. The bounding box is used
00266   // for a preliminary check of intersection.
00267   
00268   G4Point3D box_min = G4Point3D( PINFINITY);
00269   G4Point3D box_max = G4Point3D(-PINFINITY);        
00270   
00271   // Loop to search the whole control point mesh
00272   // for the minimum and maximum values for x, y and z.
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   // Do the projection to 2D
00313   G4ProjectedSurface* proj_srf = CopyToProjectedSurface(rayref);
00314 
00315   // Put surface in projected List
00316   projected_list.AddSurface(proj_srf);
00317   
00318   // Loop through List of projected surfaces
00319   while(projected_list.GetSize() > 0)
00320   {
00321     // Get first in List
00322     proj_srf = (G4ProjectedSurface*)projected_list.GetSurface();
00323     
00324     // Create the bounding box for the projected surface.
00325     proj_srf->CalcBBox();
00326     
00327 // L. Broglia   proj_srf->bbox->Test();
00328         
00329     // Check bbox test result is ok
00330     if(proj_srf->bbox->GetTestResult())
00331       // Convert the projected surface to a bezier. Split if necessary.
00332       proj_srf->ConvertToBezier(projected_list, bezier_list);
00333 
00334     // Remove projected surface
00335     projected_list.RemoveSurface(proj_srf);
00336   }
00337   
00338   // Loop through the bezier List
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     // Add a temporary Hit
00347     AddHit(bez_ptr->UAverage(), bez_ptr->VAverage());
00348         
00349     // Evaluate Hit
00350 
00351     // L. Broglia
00352     // Modify function name
00353     // bez_ptr->SetAveragePoint(Evaluate());
00354     bez_ptr->SetAveragePoint(BSEvaluate());
00355     
00356     // Calculate distance to ray origin
00357     bez_ptr->CalcDistance(rayref.GetStart());
00358 
00359     // Put closest to b_splines distance value
00360     if(bez_ptr->GetDistance() < distance) distance = bez_ptr->GetDistance();  
00361     
00362     // Remove the temporary Hit
00363     if (first_hit == Hit) first_hit = (G4UVHit*)0;
00364     delete Hit;
00365     Hit = (G4UVHit*)0;
00366     
00367     // Move to next in the List
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   // Check that approx Hit is in direction of ray
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     // Hit behind ray starting point, no intersection.
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   // Projects the nurb surface so that the z-axis = ray. 
00415   
00416   /* L. Broglia
00417   G4Point* tmp = (G4Point*)&ctl_points->get(0,0);
00418   */
00419 
00420   G4PointRat tmp = ctl_points->GetRat(0,0);
00421   G4int rational = tmp.GetType();// Get the type of control point
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 ) // 4 coordinates
00430       {
00431         G4PointRat& srfcoords = ctl_points->GetRat(i, j);
00432                 
00433 // L. Broglia   
00434 // Changes for new G4PointRat
00435 
00436         // Calculate the x- and y-coordinates for the new 
00437         // 2-D surface.  
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  // 3 coordinates
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 /* L. Broglia
00469    Changes for new G4PointRat 
00470 G4Point& G4BSplineSurface::InternalEvalCrv(G4int i, G4ControlPoints* crv)*/
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       /* L. Broglia             
00489       register G4PointRat* pts1 = &crv->get(i,j-1);
00490       register G4PointRat* pts2 = &crv->get(i,j  );         
00491       if(pts1->GetType()==3)
00492       {
00493         crv->CalcValues(k1, param, *(G4Point3D*)pts1, k2, *(G4Point3D*)pts2);
00494         crv->put(0, j, *(G4Point3D*)pts2);                           
00495       }
00496       else
00497       {
00498         crv->CalcValues(k1, param, *(G4PointRat*)pts1, k2, *(G4PointRat*)pts2);
00499         crv->put(0, j, *(G4PointRat*)pts2);           
00500       }
00501       register G4PointRat* pts1 = &crv->GetRat(i,j-1);
00502       register G4PointRat* pts2 = &crv->GetRat(i,j  );
00503       */
00504     }           
00505 
00506     j--;
00507   }     
00508 
00509   ord = ord-1;
00510   return InternalEvalCrv(0, crv); // Recursion
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   /* L. Broglia  
00523   G4Point* tmp = (G4Point*)&ctl_points->get(0,0);
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     // Evaluate the diff_curve...
00567     // G4PointRat rat_result = (G4PointRat&) InternalEvalCrv(0, diff_curve);
00568     G4PointRat rat_result(InternalEvalCrv(0, diff_curve));
00569     
00570     // Calc the 3D values.
00571     // L. Broglia
00572     // Changes for new G4PointRat
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         // G4Point3D rtr_pt  = (G4Point3D&) InternalEvalCrv(i, curves);
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       // Evaluate the diff_curve...
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   // Delete old UVhits
00615   G4UVHit* temphit=Hit;
00616   while(Hit!=(G4UVHit*)0)
00617   {
00618     Hit=Hit->GetNext();
00619     delete temphit;
00620     temphit=Hit;
00621   }
00622   
00623   // delete temphit;
00624   
00625   // Get the real Hit point
00626   closest_hit = FinalIntersection();
00627   
00628   // The following part (commented out) is old bullshit
00629   // Check that Hit is not in a void i.e. InnerBoundary.
00630   //    for(int a=0; a<NumberOfInnerBoundaries;a++)
00631   //      if(InnerBoundary[a]->Inside(closest_hit, rayref))
00632   //    {
00633   //      Active(0);
00634   //      Distance(kInfinity);
00635   //      return closest_hit;
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 }

Generated on Mon May 27 17:47:48 2013 for Geant4 by  doxygen 1.4.7