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