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