Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/management/src/G4VoxelLimits.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // Class G4VoxelLimits implementation
 27 //
 28 // 13.07.95, P.Kent - Initial version
 29 // --------------------------------------------------------------------
 30 
 31 #include "G4VoxelLimits.hh"
 32 
 33 #include "G4ios.hh"
 34 
 35 ///////////////////////////////////////////////////////////////////////////
 36 //
 37 // Further restrict limits
 38 // No checks for illegal restrictions
 39 //
 40 void G4VoxelLimits::AddLimit( const EAxis pAxis, 
 41                               const G4double pMin,
 42                               const G4double pMax )
 43 {
 44   if ( pAxis == kXAxis )
 45   {
 46     if ( pMin > fxAxisMin ) fxAxisMin = pMin ;    
 47     if ( pMax < fxAxisMax ) fxAxisMax = pMax ;    
 48   }
 49   else if ( pAxis == kYAxis )
 50   {
 51     if ( pMin > fyAxisMin ) fyAxisMin = pMin ;    
 52     if ( pMax < fyAxisMax ) fyAxisMax = pMax ;
 53   }
 54   else
 55   { 
 56     assert( pAxis == kZAxis ) ;
 57 
 58     if ( pMin > fzAxisMin ) fzAxisMin = pMin ;
 59     if ( pMax < fzAxisMax ) fzAxisMax = pMax ;
 60   }
 61 }
 62 
 63 ///////////////////////////////////////////////////////////////////////////
 64 //
 65 // ClipToLimits
 66 //
 67 // Clip the line segment pStart->pEnd to the volume described by the
 68 // current limits. Return true if the line remains after clipping,
 69 // else false, and leave the vectors in an undefined state.
 70 //
 71 // Process:
 72 //
 73 // Use Cohen-Sutherland clipping in 3D
 74 // [Fundamentals of Interactive Computer Graphics,Foley & Van Dam]
 75 //
 76 G4bool G4VoxelLimits::ClipToLimits( G4ThreeVector& pStart,
 77                                     G4ThreeVector& pEnd      ) const
 78 {
 79   G4int sCode, eCode ;
 80   G4bool remainsAfterClip ;
 81     
 82   // Determine if line is trivially inside (both outcodes==0) or outside
 83   // (logical AND of outcodes !=0)
 84 
 85   sCode = OutCode(pStart) ;
 86   eCode = OutCode(pEnd)   ;
 87 
 88   if ( (sCode & eCode) != 0 )
 89   {
 90     // Trivially outside, no intersection with region
 91 
 92     remainsAfterClip = false;
 93   }
 94   else if ( sCode == 0 && eCode == 0 )
 95   {
 96     // Trivially inside, no intersections
 97 
 98     remainsAfterClip = true ;
 99   }
100   else
101   {
102     // Line segment *may* cut volume boundaries
103     // At most, one end point is inside
104 
105     G4double x1, y1, z1, x2, y2, z2 ;
106 
107     x1 = pStart.x() ;
108     y1 = pStart.y() ;
109     z1 = pStart.z() ;
110 
111     x2 = pEnd.x() ;
112     y2 = pEnd.y() ;
113     z2 = pEnd.z() ;
114 
115     while ( sCode != eCode )  // Loop checking, 06.08.2015, G.Cosmo
116     {
117       // Copy vectors to work variables x1-z1,x2-z2
118       // Ensure x1-z1 lies outside volume, swapping vectors and outcodes
119       // if necessary
120 
121       if ( sCode != 0 )
122       {
123         if ( (sCode & 0x01) != 0 )  // Clip against fxAxisMin
124         {
125           z1 += (fxAxisMin-x1)*(z2-z1)/(x2-x1);
126           y1 += (fxAxisMin-x1)*(y2-y1)/(x2-x1);
127           x1  = fxAxisMin;
128         }
129         else if ( (sCode & 0x02) != 0 ) // Clip against fxAxisMax
130         {
131           z1 += (fxAxisMax-x1)*(z2-z1)/(x2-x1);
132           y1 += (fxAxisMax-x1)*(y2-y1)/(x2-x1);
133           x1  = fxAxisMax ;
134         }
135         else if ( (sCode & 0x04) != 0 )  // Clip against fyAxisMin
136         {
137           x1 += (fyAxisMin-y1)*(x2-x1)/(y2-y1);
138           z1 += (fyAxisMin-y1)*(z2-z1)/(y2-y1);
139           y1  = fyAxisMin;
140         }
141         else if ( (sCode & 0x08) != 0 )  // Clip against fyAxisMax
142         {
143           x1 += (fyAxisMax-y1)*(x2-x1)/(y2-y1);
144           z1 += (fyAxisMax-y1)*(z2-z1)/(y2-y1);
145           y1  = fyAxisMax;
146         }
147         else if ( (sCode & 0x10) != 0 )  // Clip against fzAxisMin
148         {
149           x1 += (fzAxisMin-z1)*(x2-x1)/(z2-z1);
150           y1 += (fzAxisMin-z1)*(y2-y1)/(z2-z1);
151           z1  = fzAxisMin;
152         }
153         else if ( (sCode & 0x20) != 0 )  // Clip against fzAxisMax
154         {
155           x1 += (fzAxisMax-z1)*(x2-x1)/(z2-z1);
156           y1 += (fzAxisMax-z1)*(y2-y1)/(z2-z1);
157           z1  = fzAxisMax;
158         }
159       }
160       if ( eCode != 0 )  // Clip 2nd end: repeat of 1st, but 1<>2
161       {
162         if ( (eCode & 0x01) != 0 )  // Clip against fxAxisMin
163         {
164           z2 += (fxAxisMin-x2)*(z1-z2)/(x1-x2);
165           y2 += (fxAxisMin-x2)*(y1-y2)/(x1-x2);
166           x2  = fxAxisMin;
167         }
168         else if ( (eCode & 0x02) != 0 )  // Clip against fxAxisMax
169         {
170           z2 += (fxAxisMax-x2)*(z1-z2)/(x1-x2);
171           y2 += (fxAxisMax-x2)*(y1-y2)/(x1-x2);
172           x2  = fxAxisMax;
173         }
174         else if ( (eCode & 0x04) != 0 )  // Clip against fyAxisMin
175         {
176           x2 += (fyAxisMin-y2)*(x1-x2)/(y1-y2);
177           z2 += (fyAxisMin-y2)*(z1-z2)/(y1-y2);
178           y2  = fyAxisMin;
179         }
180         else if ((eCode&0x08) != 0)  // Clip against fyAxisMax
181         {
182           x2 += (fyAxisMax-y2)*(x1-x2)/(y1-y2);
183           z2 += (fyAxisMax-y2)*(z1-z2)/(y1-y2);
184           y2  = fyAxisMax;
185         }
186         else if ( (eCode & 0x10) != 0 )  // Clip against fzAxisMin
187         {
188           x2 += (fzAxisMin-z2)*(x1-x2)/(z1-z2);
189           y2 += (fzAxisMin-z2)*(y1-y2)/(z1-z2);
190           z2  = fzAxisMin;
191         }
192         else if ( (eCode & 0x20) != 0 )  // Clip against fzAxisMax
193         {
194           x2 += (fzAxisMax-z2)*(x1-x2)/(z1-z2);
195           y2 += (fzAxisMax-z2)*(y1-y2)/(z1-z2);
196           z2  = fzAxisMax;
197         }
198       }
199       pStart = G4ThreeVector(x1,y1,z1);
200       pEnd   = G4ThreeVector(x2,y2,z2);
201       sCode  = OutCode(pStart);
202       eCode  = OutCode(pEnd);
203     }
204     remainsAfterClip = sCode == 0 && eCode == 0;
205   }
206   return remainsAfterClip;
207 }
208 
209 ////////////////////////////////////////////////////////////////////////////
210 //
211 // Calculate the `outcode' for the specified vector:
212 // The following bits are set:
213 //   0      pVec.x()<fxAxisMin && IsXLimited()
214 //   1      pVec.x()>fxAxisMax && IsXLimited()
215 //   2      pVec.y()<fyAxisMin && IsYLimited()
216 //   3      pVec.y()>fyAxisMax && IsYLimited()
217 //   4      pVec.z()<fzAxisMin && IsZLimited()
218 //   5      pVec.z()>fzAxisMax && IsZLimited()
219 //
220 G4int G4VoxelLimits::OutCode( const G4ThreeVector& pVec ) const
221 {
222   G4int code = 0 ;                // The outcode
223 
224   if ( IsXLimited() )
225   {
226     if ( pVec.x() < fxAxisMin ) code |= 0x01 ;
227     if ( pVec.x() > fxAxisMax ) code |= 0x02 ;
228   }
229   if ( IsYLimited() )
230   {
231     if ( pVec.y() < fyAxisMin ) code |= 0x04 ;
232     if ( pVec.y() > fyAxisMax ) code |= 0x08 ;
233   }
234   if (IsZLimited())
235   {
236     if ( pVec.z() < fzAxisMin ) code |= 0x10 ;
237     if ( pVec.z() > fzAxisMax ) code |= 0x20 ;
238   }
239   return code;
240 }
241 
242 ///////////////////////////////////////////////////////////////////////////////
243 
244 std::ostream& operator << (std::ostream& os, const G4VoxelLimits& pLim)
245 {
246     os << "{";
247     if (pLim.IsXLimited())
248         {
249             os << "(" << pLim.GetMinXExtent() 
250                << "," << pLim.GetMaxXExtent() << ") ";
251         }
252     else
253         {
254             os << "(-,-) ";
255         }
256     if (pLim.IsYLimited())
257         {
258             os << "(" << pLim.GetMinYExtent() 
259                << "," << pLim.GetMaxYExtent() << ") ";
260         }
261     else
262         {
263             os << "(-,-) ";
264         }
265     if (pLim.IsZLimited())
266         {
267             os << "(" << pLim.GetMinZExtent()
268                << "," << pLim.GetMaxZExtent() << ")";
269         }
270     else
271         {
272             os << "(-,-)";
273         }
274     os << "}";
275     return os;
276 }
277