| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // Class G4VoxelLimits implementation <<
27 // 26 //
28 // 13.07.95, P.Kent - Initial version << 27 // $Id: G4VoxelLimits.cc,v 1.11 2006/06/29 18:34:11 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-08-03-patch-01 $
>> 29 //
>> 30 // class G4VoxelLimits
>> 31 //
>> 32 // Implementation
>> 33 //
>> 34 // History:
>> 35 //
>> 36 // 14.03.02 V. Grichine, cosmetics
>> 37 // 13.07.95 P.Kent Initial version
29 // ------------------------------------------- 38 // --------------------------------------------------------------------
30 39
31 #include "G4VoxelLimits.hh" 40 #include "G4VoxelLimits.hh"
32 41
33 #include "G4ios.hh" 42 #include "G4ios.hh"
34 43
35 ////////////////////////////////////////////// 44 ///////////////////////////////////////////////////////////////////////////
36 // 45 //
>> 46 // Empty constructor and destructor
>> 47 //
>> 48
>> 49 G4VoxelLimits::G4VoxelLimits()
>> 50 : fxAxisMin(-kInfinity),fxAxisMax(kInfinity),
>> 51 fyAxisMin(-kInfinity),fyAxisMax(kInfinity),
>> 52 fzAxisMin(-kInfinity),fzAxisMax(kInfinity)
>> 53 {
>> 54 }
>> 55
>> 56 G4VoxelLimits::~G4VoxelLimits()
>> 57 {
>> 58 }
>> 59
>> 60 ///////////////////////////////////////////////////////////////////////////
>> 61 //
37 // Further restrict limits 62 // Further restrict limits
38 // No checks for illegal restrictions 63 // No checks for illegal restrictions
39 // 64 //
>> 65
40 void G4VoxelLimits::AddLimit( const EAxis pAxi 66 void G4VoxelLimits::AddLimit( const EAxis pAxis,
41 const G4double p 67 const G4double pMin,
42 const G4double p 68 const G4double pMax )
43 { 69 {
44 if ( pAxis == kXAxis ) 70 if ( pAxis == kXAxis )
45 { 71 {
46 if ( pMin > fxAxisMin ) fxAxisMin = pMin ; 72 if ( pMin > fxAxisMin ) fxAxisMin = pMin ;
47 if ( pMax < fxAxisMax ) fxAxisMax = pMax ; 73 if ( pMax < fxAxisMax ) fxAxisMax = pMax ;
48 } 74 }
49 else if ( pAxis == kYAxis ) 75 else if ( pAxis == kYAxis )
50 { 76 {
51 if ( pMin > fyAxisMin ) fyAxisMin = pMin ; 77 if ( pMin > fyAxisMin ) fyAxisMin = pMin ;
52 if ( pMax < fyAxisMax ) fyAxisMax = pMax ; 78 if ( pMax < fyAxisMax ) fyAxisMax = pMax ;
53 } 79 }
54 else 80 else
55 { 81 {
56 assert( pAxis == kZAxis ) ; 82 assert( pAxis == kZAxis ) ;
57 83
58 if ( pMin > fzAxisMin ) fzAxisMin = pMin ; 84 if ( pMin > fzAxisMin ) fzAxisMin = pMin ;
59 if ( pMax < fzAxisMax ) fzAxisMax = pMax ; 85 if ( pMax < fzAxisMax ) fzAxisMax = pMax ;
60 } 86 }
61 } 87 }
62 88
63 ////////////////////////////////////////////// 89 ///////////////////////////////////////////////////////////////////////////
64 // 90 //
65 // ClipToLimits 91 // ClipToLimits
66 // 92 //
67 // Clip the line segment pStart->pEnd to the v 93 // Clip the line segment pStart->pEnd to the volume described by the
68 // current limits. Return true if the line rem 94 // current limits. Return true if the line remains after clipping,
69 // else false, and leave the vectors in an und 95 // else false, and leave the vectors in an undefined state.
70 // 96 //
71 // Process: 97 // Process:
72 // 98 //
73 // Use Cohen-Sutherland clipping in 3D 99 // Use Cohen-Sutherland clipping in 3D
74 // [Fundamentals of Interactive Computer Graph 100 // [Fundamentals of Interactive Computer Graphics,Foley & Van Dam]
75 // 101 //
>> 102
76 G4bool G4VoxelLimits::ClipToLimits( G4ThreeVec 103 G4bool G4VoxelLimits::ClipToLimits( G4ThreeVector& pStart,
77 G4ThreeVec 104 G4ThreeVector& pEnd ) const
78 { 105 {
79 G4int sCode, eCode ; 106 G4int sCode, eCode ;
80 G4bool remainsAfterClip ; 107 G4bool remainsAfterClip ;
81 108
82 // Determine if line is trivially inside (bo 109 // Determine if line is trivially inside (both outcodes==0) or outside
83 // (logical AND of outcodes !=0) 110 // (logical AND of outcodes !=0)
84 111
85 sCode = OutCode(pStart) ; 112 sCode = OutCode(pStart) ;
86 eCode = OutCode(pEnd) ; 113 eCode = OutCode(pEnd) ;
87 114
88 if ( (sCode & eCode) != 0 ) << 115 if ( sCode & eCode )
89 { 116 {
90 // Trivially outside, no intersection with 117 // Trivially outside, no intersection with region
91 118
92 remainsAfterClip = false; 119 remainsAfterClip = false;
93 } 120 }
94 else if ( sCode == 0 && eCode == 0 ) 121 else if ( sCode == 0 && eCode == 0 )
95 { 122 {
96 // Trivially inside, no intersections 123 // Trivially inside, no intersections
97 124
98 remainsAfterClip = true ; 125 remainsAfterClip = true ;
99 } 126 }
100 else 127 else
101 { 128 {
102 // Line segment *may* cut volume boundarie 129 // Line segment *may* cut volume boundaries
103 // At most, one end point is inside 130 // At most, one end point is inside
104 131
105 G4double x1, y1, z1, x2, y2, z2 ; 132 G4double x1, y1, z1, x2, y2, z2 ;
106 133
107 x1 = pStart.x() ; 134 x1 = pStart.x() ;
108 y1 = pStart.y() ; 135 y1 = pStart.y() ;
109 z1 = pStart.z() ; 136 z1 = pStart.z() ;
110 137
111 x2 = pEnd.x() ; 138 x2 = pEnd.x() ;
112 y2 = pEnd.y() ; 139 y2 = pEnd.y() ;
113 z2 = pEnd.z() ; 140 z2 = pEnd.z() ;
114 << 141 /*
115 while ( sCode != eCode ) // Loop checking << 142 if( std::abs(x1-x2) < kCarTolerance*kCarTolerance)
>> 143 {
>> 144 G4cout<<"x1 = "<<x1<<"\t"<<"x2 = "<<x2<<G4endl;
>> 145 }
>> 146 if( std::abs(y1-y2) < kCarTolerance*kCarTolerance)
>> 147 {
>> 148 G4cout<<"y1 = "<<y1<<"\t"<<"y2 = "<<y2<<G4endl;
>> 149 }
>> 150 if( std::abs(z1-z2) < kCarTolerance*kCarTolerance)
>> 151 {
>> 152 G4cout<<"z1 = "<<z1<<"\t"<<"z2 = "<<z2<<G4endl;
>> 153 }
>> 154 */
>> 155 while ( sCode != eCode )
116 { 156 {
117 // Copy vectors to work variables x1-z1, 157 // Copy vectors to work variables x1-z1,x2-z2
118 // Ensure x1-z1 lies outside volume, swa 158 // Ensure x1-z1 lies outside volume, swapping vectors and outcodes
119 // if necessary 159 // if necessary
120 160
121 if ( sCode != 0 ) << 161 if ( sCode )
122 { 162 {
123 if ( (sCode & 0x01) != 0 ) // Clip ag << 163 if ( sCode & 0x01 ) // Clip against fxAxisMin
124 { 164 {
125 z1 += (fxAxisMin-x1)*(z2-z1)/(x2-x1) 165 z1 += (fxAxisMin-x1)*(z2-z1)/(x2-x1);
126 y1 += (fxAxisMin-x1)*(y2-y1)/(x2-x1) 166 y1 += (fxAxisMin-x1)*(y2-y1)/(x2-x1);
127 x1 = fxAxisMin; 167 x1 = fxAxisMin;
128 } 168 }
129 else if ( (sCode & 0x02) != 0 ) // Cli << 169 else if ( sCode & 0x02 ) // Clip against fxAxisMax
130 { 170 {
131 z1 += (fxAxisMax-x1)*(z2-z1)/(x2-x1) 171 z1 += (fxAxisMax-x1)*(z2-z1)/(x2-x1);
132 y1 += (fxAxisMax-x1)*(y2-y1)/(x2-x1) 172 y1 += (fxAxisMax-x1)*(y2-y1)/(x2-x1);
133 x1 = fxAxisMax ; 173 x1 = fxAxisMax ;
134 } 174 }
135 else if ( (sCode & 0x04) != 0 ) // Cl << 175 else if ( sCode & 0x04 ) // Clip against fyAxisMin
136 { 176 {
137 x1 += (fyAxisMin-y1)*(x2-x1)/(y2-y1) 177 x1 += (fyAxisMin-y1)*(x2-x1)/(y2-y1);
138 z1 += (fyAxisMin-y1)*(z2-z1)/(y2-y1) 178 z1 += (fyAxisMin-y1)*(z2-z1)/(y2-y1);
139 y1 = fyAxisMin; 179 y1 = fyAxisMin;
140 } 180 }
141 else if ( (sCode & 0x08) != 0 ) // Cl << 181 else if ( sCode & 0x08 ) // Clip against fyAxisMax
142 { 182 {
143 x1 += (fyAxisMax-y1)*(x2-x1)/(y2-y1) 183 x1 += (fyAxisMax-y1)*(x2-x1)/(y2-y1);
144 z1 += (fyAxisMax-y1)*(z2-z1)/(y2-y1) 184 z1 += (fyAxisMax-y1)*(z2-z1)/(y2-y1);
145 y1 = fyAxisMax; 185 y1 = fyAxisMax;
146 } 186 }
147 else if ( (sCode & 0x10) != 0 ) // Cl << 187 else if ( sCode & 0x10 ) // Clip against fzAxisMin
148 { 188 {
149 x1 += (fzAxisMin-z1)*(x2-x1)/(z2-z1) 189 x1 += (fzAxisMin-z1)*(x2-x1)/(z2-z1);
150 y1 += (fzAxisMin-z1)*(y2-y1)/(z2-z1) 190 y1 += (fzAxisMin-z1)*(y2-y1)/(z2-z1);
151 z1 = fzAxisMin; 191 z1 = fzAxisMin;
152 } 192 }
153 else if ( (sCode & 0x20) != 0 ) // Cl << 193 else if ( sCode & 0x20 ) // Clip against fzAxisMax
154 { 194 {
155 x1 += (fzAxisMax-z1)*(x2-x1)/(z2-z1) 195 x1 += (fzAxisMax-z1)*(x2-x1)/(z2-z1);
156 y1 += (fzAxisMax-z1)*(y2-y1)/(z2-z1) 196 y1 += (fzAxisMax-z1)*(y2-y1)/(z2-z1);
157 z1 = fzAxisMax; 197 z1 = fzAxisMax;
158 } 198 }
159 } 199 }
160 if ( eCode != 0 ) // Clip 2nd end: repe << 200 if ( eCode ) // Clip 2nd end: repeat of 1st, but 1<>2
161 { 201 {
162 if ( (eCode & 0x01) != 0 ) // Clip ag << 202 if ( eCode & 0x01 ) // Clip against fxAxisMin
163 { 203 {
164 z2 += (fxAxisMin-x2)*(z1-z2)/(x1-x2) 204 z2 += (fxAxisMin-x2)*(z1-z2)/(x1-x2);
165 y2 += (fxAxisMin-x2)*(y1-y2)/(x1-x2) 205 y2 += (fxAxisMin-x2)*(y1-y2)/(x1-x2);
166 x2 = fxAxisMin; 206 x2 = fxAxisMin;
167 } 207 }
168 else if ( (eCode & 0x02) != 0 ) // Cl << 208 else if ( eCode & 0x02 ) // Clip against fxAxisMax
169 { 209 {
170 z2 += (fxAxisMax-x2)*(z1-z2)/(x1-x2) 210 z2 += (fxAxisMax-x2)*(z1-z2)/(x1-x2);
171 y2 += (fxAxisMax-x2)*(y1-y2)/(x1-x2) 211 y2 += (fxAxisMax-x2)*(y1-y2)/(x1-x2);
172 x2 = fxAxisMax; 212 x2 = fxAxisMax;
173 } 213 }
174 else if ( (eCode & 0x04) != 0 ) // Cl << 214 else if ( eCode & 0x04 ) // Clip against fyAxisMin
175 { 215 {
176 x2 += (fyAxisMin-y2)*(x1-x2)/(y1-y2) 216 x2 += (fyAxisMin-y2)*(x1-x2)/(y1-y2);
177 z2 += (fyAxisMin-y2)*(z1-z2)/(y1-y2) 217 z2 += (fyAxisMin-y2)*(z1-z2)/(y1-y2);
178 y2 = fyAxisMin; 218 y2 = fyAxisMin;
179 } 219 }
180 else if ((eCode&0x08) != 0) // Clip a << 220 else if (eCode&0x08) // Clip against fyAxisMax
181 { 221 {
182 x2 += (fyAxisMax-y2)*(x1-x2)/(y1-y2) 222 x2 += (fyAxisMax-y2)*(x1-x2)/(y1-y2);
183 z2 += (fyAxisMax-y2)*(z1-z2)/(y1-y2) 223 z2 += (fyAxisMax-y2)*(z1-z2)/(y1-y2);
184 y2 = fyAxisMax; 224 y2 = fyAxisMax;
185 } 225 }
186 else if ( (eCode & 0x10) != 0 ) // Cl << 226 else if ( eCode & 0x10 ) // Clip against fzAxisMin
187 { 227 {
188 x2 += (fzAxisMin-z2)*(x1-x2)/(z1-z2) 228 x2 += (fzAxisMin-z2)*(x1-x2)/(z1-z2);
189 y2 += (fzAxisMin-z2)*(y1-y2)/(z1-z2) 229 y2 += (fzAxisMin-z2)*(y1-y2)/(z1-z2);
190 z2 = fzAxisMin; 230 z2 = fzAxisMin;
191 } 231 }
192 else if ( (eCode & 0x20) != 0 ) // Cl << 232 else if ( eCode & 0x20 ) // Clip against fzAxisMax
193 { 233 {
194 x2 += (fzAxisMax-z2)*(x1-x2)/(z1-z2) 234 x2 += (fzAxisMax-z2)*(x1-x2)/(z1-z2);
195 y2 += (fzAxisMax-z2)*(y1-y2)/(z1-z2) 235 y2 += (fzAxisMax-z2)*(y1-y2)/(z1-z2);
196 z2 = fzAxisMax; 236 z2 = fzAxisMax;
197 } 237 }
198 } 238 }
>> 239 // G4endl; G4cout<<"x1 = "<<x1<<"\t"<<"x2 = "<<x2<<G4endl<<G4endl;
199 pStart = G4ThreeVector(x1,y1,z1); 240 pStart = G4ThreeVector(x1,y1,z1);
200 pEnd = G4ThreeVector(x2,y2,z2); 241 pEnd = G4ThreeVector(x2,y2,z2);
201 sCode = OutCode(pStart); 242 sCode = OutCode(pStart);
202 eCode = OutCode(pEnd); 243 eCode = OutCode(pEnd);
203 } 244 }
204 remainsAfterClip = sCode == 0 && eCode == << 245 if ( sCode == 0 && eCode == 0 ) remainsAfterClip = true;
>> 246 else remainsAfterClip = false;
205 } 247 }
206 return remainsAfterClip; 248 return remainsAfterClip;
207 } 249 }
208 250
209 ////////////////////////////////////////////// 251 ////////////////////////////////////////////////////////////////////////////
210 // 252 //
211 // Calculate the `outcode' for the specified v 253 // Calculate the `outcode' for the specified vector:
212 // The following bits are set: 254 // The following bits are set:
213 // 0 pVec.x()<fxAxisMin && IsXLimited() 255 // 0 pVec.x()<fxAxisMin && IsXLimited()
214 // 1 pVec.x()>fxAxisMax && IsXLimited() 256 // 1 pVec.x()>fxAxisMax && IsXLimited()
215 // 2 pVec.y()<fyAxisMin && IsYLimited() 257 // 2 pVec.y()<fyAxisMin && IsYLimited()
216 // 3 pVec.y()>fyAxisMax && IsYLimited() 258 // 3 pVec.y()>fyAxisMax && IsYLimited()
217 // 4 pVec.z()<fzAxisMin && IsZLimited() 259 // 4 pVec.z()<fzAxisMin && IsZLimited()
218 // 5 pVec.z()>fzAxisMax && IsZLimited() 260 // 5 pVec.z()>fzAxisMax && IsZLimited()
219 // 261 //
>> 262
220 G4int G4VoxelLimits::OutCode( const G4ThreeVec 263 G4int G4VoxelLimits::OutCode( const G4ThreeVector& pVec ) const
221 { 264 {
222 G4int code = 0 ; // The outco 265 G4int code = 0 ; // The outcode
223 266
224 if ( IsXLimited() ) 267 if ( IsXLimited() )
225 { 268 {
226 if ( pVec.x() < fxAxisMin ) code |= 0x01 ; 269 if ( pVec.x() < fxAxisMin ) code |= 0x01 ;
227 if ( pVec.x() > fxAxisMax ) code |= 0x02 ; 270 if ( pVec.x() > fxAxisMax ) code |= 0x02 ;
228 } 271 }
229 if ( IsYLimited() ) 272 if ( IsYLimited() )
230 { 273 {
231 if ( pVec.y() < fyAxisMin ) code |= 0x04 ; 274 if ( pVec.y() < fyAxisMin ) code |= 0x04 ;
232 if ( pVec.y() > fyAxisMax ) code |= 0x08 ; 275 if ( pVec.y() > fyAxisMax ) code |= 0x08 ;
233 } 276 }
234 if (IsZLimited()) 277 if (IsZLimited())
235 { 278 {
236 if ( pVec.z() < fzAxisMin ) code |= 0x10 ; 279 if ( pVec.z() < fzAxisMin ) code |= 0x10 ;
237 if ( pVec.z() > fzAxisMax ) code |= 0x20 ; 280 if ( pVec.z() > fzAxisMax ) code |= 0x20 ;
238 } 281 }
239 return code; 282 return code;
240 } 283 }
241 284
242 ////////////////////////////////////////////// 285 ///////////////////////////////////////////////////////////////////////////////
243 286
244 std::ostream& operator << (std::ostream& os, c 287 std::ostream& operator << (std::ostream& os, const G4VoxelLimits& pLim)
245 { 288 {
246 os << "{"; 289 os << "{";
247 if (pLim.IsXLimited()) 290 if (pLim.IsXLimited())
248 { 291 {
249 os << "(" << pLim.GetMinXExtent() 292 os << "(" << pLim.GetMinXExtent()
250 << "," << pLim.GetMaxXExtent() 293 << "," << pLim.GetMaxXExtent() << ") ";
251 } 294 }
252 else 295 else
253 { 296 {
254 os << "(-,-) "; 297 os << "(-,-) ";
255 } 298 }
256 if (pLim.IsYLimited()) 299 if (pLim.IsYLimited())
257 { 300 {
258 os << "(" << pLim.GetMinYExtent() 301 os << "(" << pLim.GetMinYExtent()
259 << "," << pLim.GetMaxYExtent() 302 << "," << pLim.GetMaxYExtent() << ") ";
260 } 303 }
261 else 304 else
262 { 305 {
263 os << "(-,-) "; 306 os << "(-,-) ";
264 } 307 }
265 if (pLim.IsZLimited()) 308 if (pLim.IsZLimited())
266 { 309 {
267 os << "(" << pLim.GetMinZExtent() 310 os << "(" << pLim.GetMinZExtent()
268 << "," << pLim.GetMaxZExtent() 311 << "," << pLim.GetMaxZExtent() << ")";
269 } 312 }
270 else 313 else
271 { 314 {
272 os << "(-,-)"; 315 os << "(-,-)";
273 } 316 }
274 os << "}"; 317 os << "}";
275 return os; 318 return os;
276 } 319 }
277 320