G4ClippablePolygon.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: G4ClippablePolygon.cc 67011 2013-01-29 16:17:41Z gcosmo $
00028 //
00029 // 
00030 // --------------------------------------------------------------------
00031 // GEANT 4 class source file
00032 //
00033 //
00034 // G4ClippablePolygon.cc
00035 //
00036 // Includes code from G4VSolid (P.Kent, V.Grichine, J.Allison)
00037 //
00038 // --------------------------------------------------------------------
00039 
00040 #include "G4ClippablePolygon.hh"
00041 
00042 #include "G4VoxelLimits.hh"
00043 #include "G4GeometryTolerance.hh"
00044 
00045 //
00046 // Constructor
00047 //
00048 G4ClippablePolygon::G4ClippablePolygon()
00049   : normal(0.,0.,0.)
00050 {
00051   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00052 }
00053 
00054 
00055 //
00056 // Destructor
00057 //
00058 G4ClippablePolygon::~G4ClippablePolygon()
00059 {
00060 }
00061 
00062 
00063 //
00064 // AddVertexInOrder
00065 //
00066 void G4ClippablePolygon::AddVertexInOrder( const G4ThreeVector vertex )
00067 {
00068   vertices.push_back( vertex );
00069 }
00070 
00071 
00072 //
00073 // ClearAllVertices
00074 //
00075 void G4ClippablePolygon::ClearAllVertices()
00076 {
00077   vertices.clear();
00078 }
00079 
00080 
00081 //
00082 // Clip
00083 //
00084 G4bool G4ClippablePolygon::Clip( const G4VoxelLimits &voxelLimit )
00085 {
00086   if (voxelLimit.IsLimited()) {
00087     ClipAlongOneAxis( voxelLimit, kXAxis );
00088     ClipAlongOneAxis( voxelLimit, kYAxis );
00089     ClipAlongOneAxis( voxelLimit, kZAxis );
00090   }
00091   
00092   return (vertices.size() > 0);
00093 }
00094 
00095 
00096 //
00097 // PartialClip
00098 //
00099 // Clip, while ignoring the indicated axis
00100 //
00101 G4bool G4ClippablePolygon::PartialClip( const G4VoxelLimits &voxelLimit,
00102                                         const EAxis IgnoreMe )
00103 {
00104   if (voxelLimit.IsLimited()) {
00105     if (IgnoreMe != kXAxis) ClipAlongOneAxis( voxelLimit, kXAxis );
00106     if (IgnoreMe != kYAxis) ClipAlongOneAxis( voxelLimit, kYAxis );
00107     if (IgnoreMe != kZAxis) ClipAlongOneAxis( voxelLimit, kZAxis );
00108   }
00109   
00110   return (vertices.size() > 0);
00111 }
00112 
00113 
00114 //
00115 // GetExtent
00116 //
00117 G4bool G4ClippablePolygon::GetExtent( const EAxis axis, 
00118                                             G4double &min,
00119                                             G4double &max ) const
00120 {
00121   //
00122   // Okay, how many entries do we have?
00123   //
00124   G4int noLeft = vertices.size();
00125   
00126   //
00127   // Return false if nothing is left
00128   //
00129   if (noLeft == 0) return false;
00130   
00131   //
00132   // Initialize min and max to our first vertex
00133   //
00134   min = max = vertices[0].operator()( axis );
00135   
00136   //
00137   // Compare to the rest
00138   //
00139   G4int i;
00140   for( i=1; i<noLeft; i++ )
00141   {
00142     G4double component = vertices[i].operator()( axis );
00143     if (component < min )
00144       min = component;
00145     else if (component > max )
00146       max = component;
00147   }
00148   
00149   return true;
00150 }
00151 
00152 
00153 //
00154 // GetMinPoint
00155 //
00156 // Returns pointer to minimum point along the specified axis.
00157 // Take care! Do not use pointer after destroying parent polygon.
00158 //
00159 const G4ThreeVector *G4ClippablePolygon::GetMinPoint( const EAxis axis ) const
00160 {
00161   G4int noLeft = vertices.size();
00162   if (noLeft==0)
00163     G4Exception("G4ClippablePolygon::GetMinPoint()",
00164                 "GeomSolids0002", FatalException, "Empty polygon.");
00165   
00166   const G4ThreeVector *answer = &(vertices[0]);
00167   G4double min = answer->operator()(axis);
00168 
00169   G4int i;
00170   for( i=1; i<noLeft; i++ )
00171   {
00172     G4double component = vertices[i].operator()( axis );
00173     if (component < min)
00174     {
00175       answer = &(vertices[i]);
00176       min = component;
00177     }
00178   }
00179   
00180   return answer;
00181 }
00182 
00183 
00184 //
00185 // GetMaxPoint
00186 //
00187 // Returns pointer to maximum point along the specified axis.
00188 // Take care! Do not use pointer after destroying parent polygon.
00189 //
00190 const G4ThreeVector *G4ClippablePolygon::GetMaxPoint( const EAxis axis ) const
00191 {
00192   G4int noLeft = vertices.size();
00193   if (noLeft==0)
00194     G4Exception("G4ClippablePolygon::GetMaxPoint()",
00195                 "GeomSolids0002", FatalException, "Empty polygon.");
00196   
00197   const G4ThreeVector *answer = &(vertices[0]);
00198   G4double max = answer->operator()(axis);
00199 
00200   G4int i;
00201   for( i=1; i<noLeft; i++ )
00202   {
00203     G4double component = vertices[i].operator()( axis );
00204     if (component > max)
00205     {
00206       answer = &(vertices[i]);
00207       max = component;
00208     }
00209   }
00210   
00211   return answer;
00212 }
00213     
00214 
00215 //
00216 // InFrontOf
00217 //
00218 // Decide if this polygon is in "front" of another when
00219 // viewed along the specified axis. For our purposes here,
00220 // it is sufficient to use the minimum extent of the
00221 // polygon along the axis to determine this.
00222 //
00223 // In case the minima of the two polygons are equal,
00224 // we use a more sophisticated test.
00225 //
00226 // Note that it is possible for the two following
00227 // statements to both return true or both return false:
00228 //         polygon1.InFrontOf(polygon2)
00229 //         polygon2.BehindOf(polygon1)
00230 //
00231 G4bool G4ClippablePolygon::InFrontOf( const G4ClippablePolygon &other,
00232                                             EAxis axis ) const
00233 {
00234   //
00235   // If things are empty, do something semi-sensible
00236   //
00237   G4int noLeft = vertices.size();
00238   if (noLeft==0) return false;
00239   
00240   if (other.Empty()) return true;
00241 
00242   //
00243   // Get minimum of other polygon
00244   //
00245   const G4ThreeVector *minPointOther = other.GetMinPoint( axis );
00246   const G4double minOther = minPointOther->operator()(axis);
00247   
00248   //
00249   // Get minimum of this polygon
00250   //
00251   const G4ThreeVector *minPoint = GetMinPoint( axis );
00252   const G4double min = minPoint->operator()(axis);
00253   
00254   //
00255   // Easy decision
00256   //
00257   if (min < minOther-kCarTolerance) return true;    // Clear winner
00258   
00259   if (minOther < min-kCarTolerance) return false;    // Clear loser
00260   
00261   //
00262   // We have a tie (this will not be all that rare since our
00263   // polygons are connected)
00264   //
00265   // Check to see if there is a vertex in the other polygon
00266   // that is behind this one (or vice versa)
00267   //
00268   G4bool answer;
00269   G4ThreeVector normalOther = other.GetNormal();
00270   
00271   if (std::fabs(normalOther(axis)) > std::fabs(normal(axis)))
00272   {
00273     G4double minP, maxP;
00274     GetPlanerExtent( *minPointOther, normalOther, minP, maxP );
00275     
00276     answer = (normalOther(axis) > 0) ? (minP < -kCarTolerance)
00277                                      : (maxP > +kCarTolerance);
00278   }
00279   else
00280   {
00281     G4double minP, maxP;
00282     other.GetPlanerExtent( *minPoint, normal, minP, maxP );
00283     
00284     answer = (normal(axis) > 0) ? (maxP > +kCarTolerance)
00285                                 : (minP < -kCarTolerance);
00286   }
00287   return answer;
00288 }
00289 
00290 //
00291 // BehindOf
00292 //
00293 // Decide if this polygon is behind another.
00294 // See notes in method "InFrontOf"
00295 //
00296 G4bool G4ClippablePolygon::BehindOf( const G4ClippablePolygon &other,
00297                                            EAxis axis ) const
00298 {
00299   //
00300   // If things are empty, do something semi-sensible
00301   //
00302   G4int noLeft = vertices.size();
00303   if (noLeft==0) return false;
00304   
00305   if (other.Empty()) return true;
00306 
00307   //
00308   // Get minimum of other polygon
00309   //
00310   const G4ThreeVector *maxPointOther = other.GetMaxPoint( axis );
00311   const G4double maxOther = maxPointOther->operator()(axis);
00312   
00313   //
00314   // Get minimum of this polygon
00315   //
00316   const G4ThreeVector *maxPoint = GetMaxPoint( axis );
00317   const G4double max = maxPoint->operator()(axis);
00318   
00319   //
00320   // Easy decision
00321   //
00322   if (max > maxOther+kCarTolerance) return true;    // Clear winner
00323   
00324   if (maxOther > max+kCarTolerance) return false;    // Clear loser
00325   
00326   //
00327   // We have a tie (this will not be all that rare since our
00328   // polygons are connected)
00329   //
00330   // Check to see if there is a vertex in the other polygon
00331   // that is in front of this one (or vice versa)
00332   //
00333   G4bool answer;
00334   G4ThreeVector normalOther = other.GetNormal();
00335   
00336   if (std::fabs(normalOther(axis)) > std::fabs(normal(axis)))
00337   {
00338     G4double minP, maxP;
00339     GetPlanerExtent( *maxPointOther, normalOther, minP, maxP );
00340     
00341     answer = (normalOther(axis) > 0) ? (maxP > +kCarTolerance)
00342                                      : (minP < -kCarTolerance);
00343   }
00344   else
00345   {
00346     G4double minP, maxP;
00347     other.GetPlanerExtent( *maxPoint, normal, minP, maxP );
00348     
00349     answer = (normal(axis) > 0) ? (minP < -kCarTolerance)
00350                                 : (maxP > +kCarTolerance);
00351   }
00352   return answer;
00353 }
00354 
00355 
00356 //
00357 // GetPlanerExtent
00358 //
00359 // Get min/max distance in or out of a plane
00360 //
00361 G4bool G4ClippablePolygon::GetPlanerExtent( const G4ThreeVector &pointOnPlane, 
00362                                             const G4ThreeVector &planeNormal,
00363                                                   G4double &min,
00364                                                   G4double &max ) const
00365 {
00366   //
00367   // Okay, how many entries do we have?
00368   //
00369   G4int noLeft = vertices.size();
00370   
00371   //
00372   // Return false if nothing is left
00373   //
00374   if (noLeft == 0) return false;
00375   
00376   //
00377   // Initialize min and max to our first vertex
00378   //
00379   min = max = planeNormal.dot(vertices[0]-pointOnPlane);
00380   
00381   //
00382   // Compare to the rest
00383   //
00384   G4int i;
00385   for( i=1; i<noLeft; i++ )
00386   {
00387     G4double component = planeNormal.dot(vertices[i] - pointOnPlane);
00388     if (component < min )
00389       min = component;
00390     else if (component > max )
00391       max = component;
00392   }
00393   
00394   return true;
00395 }
00396 
00397 
00398 //
00399 // Clip along just one axis, as specified in voxelLimit
00400 //
00401 void G4ClippablePolygon::ClipAlongOneAxis( const G4VoxelLimits &voxelLimit,
00402                                            const EAxis axis )
00403 {    
00404   if (!voxelLimit.IsLimited(axis)) return;
00405   
00406   G4ThreeVectorList tempPolygon;
00407 
00408   //
00409   // Build a "simple" voxelLimit that includes only the min extent
00410   // and apply this to our vertices, producing result in tempPolygon
00411   //
00412   G4VoxelLimits simpleLimit1;
00413   simpleLimit1.AddLimit( axis, voxelLimit.GetMinExtent(axis), kInfinity );
00414   ClipToSimpleLimits( vertices, tempPolygon, simpleLimit1 );
00415 
00416   //
00417   // If nothing is left from the above clip, we might as well return now
00418   // (but with an empty vertices)
00419   //
00420   if (tempPolygon.size() == 0)
00421   {
00422     vertices.clear();
00423     return;
00424   }
00425 
00426   //
00427   // Now do the same, but using a "simple" limit that includes only the max
00428   // extent. Apply this to out tempPolygon, producing result in vertices.
00429   //
00430   G4VoxelLimits simpleLimit2;
00431   simpleLimit2.AddLimit( axis, -kInfinity, voxelLimit.GetMaxExtent(axis) );
00432   ClipToSimpleLimits( tempPolygon, vertices, simpleLimit2 );
00433 
00434   //
00435   // If nothing is left, return now
00436   //
00437   if (vertices.size() == 0) return;
00438 }
00439 
00440 
00441 //
00442 // pVoxelLimits must be only limited along one axis, and either the maximum
00443 // along the axis must be +kInfinity, or the minimum -kInfinity
00444 //
00445 void G4ClippablePolygon::ClipToSimpleLimits( G4ThreeVectorList& pPolygon,
00446                                              G4ThreeVectorList& outputPolygon,
00447                                        const G4VoxelLimits& pVoxelLimit   )
00448 {
00449   G4int i;
00450   G4int noVertices=pPolygon.size();
00451   G4ThreeVector vEnd,vStart;
00452 
00453   outputPolygon.clear();
00454     
00455   for (i=0;i<noVertices;i++)
00456   {
00457     vStart=pPolygon[i];
00458     if (i==noVertices-1)
00459     {
00460       vEnd=pPolygon[0];
00461     }
00462     else
00463     {
00464       vEnd=pPolygon[i+1];
00465     }
00466 
00467     if (pVoxelLimit.Inside(vStart))
00468     {
00469       if (pVoxelLimit.Inside(vEnd))
00470       {
00471         // vStart and vEnd inside -> output end point
00472         //
00473         outputPolygon.push_back(vEnd);
00474       }
00475       else
00476       {
00477         // vStart inside, vEnd outside -> output crossing point
00478         //
00479         pVoxelLimit.ClipToLimits(vStart,vEnd);
00480         outputPolygon.push_back(vEnd);
00481       }
00482     }
00483     else
00484     {
00485       if (pVoxelLimit.Inside(vEnd))
00486       {
00487         // vStart outside, vEnd inside -> output inside section
00488         //
00489         pVoxelLimit.ClipToLimits(vStart,vEnd);
00490         outputPolygon.push_back(vStart);
00491         outputPolygon.push_back(vEnd);
00492       }
00493       else    // Both point outside -> no output
00494       {
00495       }
00496     }
00497   }
00498 }

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