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