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 // Implementation of G4BoundingEnvelope 26 // Implementation of G4BoundingEnvelope 27 // 27 // 28 // 2016.05.25, E.Tcherniaev - initial version 28 // 2016.05.25, E.Tcherniaev - initial version 29 // ------------------------------------------- 29 // -------------------------------------------------------------------- 30 30 31 #include <cmath> 31 #include <cmath> 32 << 33 #include "globals.hh" 32 #include "globals.hh" >> 33 34 #include "G4BoundingEnvelope.hh" 34 #include "G4BoundingEnvelope.hh" 35 #include "G4GeometryTolerance.hh" 35 #include "G4GeometryTolerance.hh" 36 #include "G4Normal3D.hh" << 37 36 38 const G4double kCarTolerance = 37 const G4double kCarTolerance = 39 G4GeometryTolerance::GetInstance()->GetSurfa 38 G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 40 39 41 ////////////////////////////////////////////// 40 /////////////////////////////////////////////////////////////////////// 42 // 41 // 43 // Constructor from an axis aligned bounding b 42 // Constructor from an axis aligned bounding box 44 // 43 // 45 G4BoundingEnvelope::G4BoundingEnvelope(const G 44 G4BoundingEnvelope::G4BoundingEnvelope(const G4ThreeVector& pMin, 46 const G 45 const G4ThreeVector& pMax) 47 : fMin(pMin), fMax(pMax) 46 : fMin(pMin), fMax(pMax) 48 { 47 { 49 // Check correctness of bounding box 48 // Check correctness of bounding box 50 // 49 // 51 CheckBoundingBox(); 50 CheckBoundingBox(); 52 } 51 } 53 52 54 ////////////////////////////////////////////// 53 /////////////////////////////////////////////////////////////////////// 55 // 54 // 56 // Constructor from a sequence of polygons 55 // Constructor from a sequence of polygons 57 // 56 // 58 G4BoundingEnvelope:: 57 G4BoundingEnvelope:: 59 G4BoundingEnvelope(const std::vector<const G4T 58 G4BoundingEnvelope(const std::vector<const G4ThreeVectorList*>& polygons) 60 : fPolygons(&polygons) 59 : fPolygons(&polygons) 61 { 60 { 62 // Check correctness of polygons 61 // Check correctness of polygons 63 // 62 // 64 CheckBoundingPolygons(); 63 CheckBoundingPolygons(); 65 64 66 // Set bounding box 65 // Set bounding box 67 // 66 // 68 G4double xmin = kInfinity, ymin = kInfinit 67 G4double xmin = kInfinity, ymin = kInfinity, zmin = kInfinity; 69 G4double xmax = -kInfinity, ymax = -kInfinit 68 G4double xmax = -kInfinity, ymax = -kInfinity, zmax = -kInfinity; 70 for (const auto & polygon : *fPolygons) << 69 for (auto ibase = fPolygons->cbegin(); ibase != fPolygons->cend(); ++ibase) 71 { << 70 { 72 for (auto ipoint = polygon->cbegin(); ipoi << 71 for (auto ipoint = (*ibase)->cbegin(); ipoint != (*ibase)->cend(); ++ipoint) 73 { 72 { 74 G4double x = ipoint->x(); << 73 G4double x = ipoint->x(); 75 if (x < xmin) xmin = x; 74 if (x < xmin) xmin = x; 76 if (x > xmax) xmax = x; 75 if (x > xmax) xmax = x; 77 G4double y = ipoint->y(); << 76 G4double y = ipoint->y(); 78 if (y < ymin) ymin = y; 77 if (y < ymin) ymin = y; 79 if (y > ymax) ymax = y; 78 if (y > ymax) ymax = y; 80 G4double z = ipoint->z(); << 79 G4double z = ipoint->z(); 81 if (z < zmin) zmin = z; 80 if (z < zmin) zmin = z; 82 if (z > zmax) zmax = z; 81 if (z > zmax) zmax = z; 83 } 82 } 84 } 83 } 85 fMin.set(xmin,ymin,zmin); 84 fMin.set(xmin,ymin,zmin); 86 fMax.set(xmax,ymax,zmax); 85 fMax.set(xmax,ymax,zmax); 87 86 88 // Check correctness of bounding box 87 // Check correctness of bounding box 89 // 88 // 90 CheckBoundingBox(); 89 CheckBoundingBox(); 91 } 90 } 92 91 93 ////////////////////////////////////////////// 92 /////////////////////////////////////////////////////////////////////// 94 // 93 // 95 // Constructor from a bounding box and a seque 94 // Constructor from a bounding box and a sequence of polygons 96 // 95 // 97 G4BoundingEnvelope:: 96 G4BoundingEnvelope:: 98 G4BoundingEnvelope( const G4ThreeVector& pMin, 97 G4BoundingEnvelope( const G4ThreeVector& pMin, 99 const G4ThreeVector& pMax, 98 const G4ThreeVector& pMax, 100 const std::vector<const G4 99 const std::vector<const G4ThreeVectorList*>& polygons) 101 : fMin(pMin), fMax(pMax), fPolygons(&polygon 100 : fMin(pMin), fMax(pMax), fPolygons(&polygons) 102 { 101 { 103 // Check correctness of bounding box and pol 102 // Check correctness of bounding box and polygons 104 // << 103 // 105 CheckBoundingBox(); 104 CheckBoundingBox(); 106 CheckBoundingPolygons(); 105 CheckBoundingPolygons(); 107 } 106 } 108 107 109 ////////////////////////////////////////////// 108 /////////////////////////////////////////////////////////////////////// 110 // 109 // >> 110 // Destructor >> 111 // >> 112 G4BoundingEnvelope::~G4BoundingEnvelope() >> 113 { >> 114 } >> 115 >> 116 /////////////////////////////////////////////////////////////////////// >> 117 // 111 // Check correctness of the axis aligned bound 118 // Check correctness of the axis aligned bounding box 112 // 119 // 113 void G4BoundingEnvelope::CheckBoundingBox() 120 void G4BoundingEnvelope::CheckBoundingBox() 114 { 121 { 115 if (fMin.x() >= fMax.x() || fMin.y() >= fMax 122 if (fMin.x() >= fMax.x() || fMin.y() >= fMax.y() || fMin.z() >= fMax.z()) 116 { 123 { 117 std::ostringstream message; 124 std::ostringstream message; 118 message << "Badly defined bounding box (mi 125 message << "Badly defined bounding box (min >= max)!" 119 << "\npMin = " << fMin << 126 << "\npMin = " << fMin 120 << "\npMax = " << fMax; 127 << "\npMax = " << fMax; 121 G4Exception("G4BoundingEnvelope::CheckBoun 128 G4Exception("G4BoundingEnvelope::CheckBoundingBox()", 122 "GeomMgt0001", JustWarning, me 129 "GeomMgt0001", JustWarning, message); 123 } 130 } 124 } 131 } 125 132 126 ////////////////////////////////////////////// 133 /////////////////////////////////////////////////////////////////////// 127 // 134 // 128 // Check correctness of the sequence of boundi 135 // Check correctness of the sequence of bounding polygons. 129 // Firsf and last polygons may consist of a si 136 // Firsf and last polygons may consist of a single vertex 130 // 137 // 131 void G4BoundingEnvelope::CheckBoundingPolygons 138 void G4BoundingEnvelope::CheckBoundingPolygons() 132 { 139 { 133 std::size_t nbases = fPolygons->size(); << 140 G4int nbases = fPolygons->size(); 134 if (nbases < 2) 141 if (nbases < 2) 135 { 142 { 136 std::ostringstream message; 143 std::ostringstream message; 137 message << "Wrong number of polygons in th 144 message << "Wrong number of polygons in the sequence: " << nbases 138 << "\nShould be at least two!"; 145 << "\nShould be at least two!"; 139 G4Exception("G4BoundingEnvelope::CheckBoun << 146 G4Exception("G4BoundingEnvelope::CheckBoundingPolygons()", 140 "GeomMgt0001", FatalException, 147 "GeomMgt0001", FatalException, message); 141 return; 148 return; 142 } 149 } 143 150 144 std::size_t nsize = std::max((*fPolygons)[0 << 151 G4int nsize = std::max((*fPolygons)[0]->size(),(*fPolygons)[1]->size()); 145 if (nsize < 3) 152 if (nsize < 3) 146 { 153 { 147 std::ostringstream message; 154 std::ostringstream message; 148 message << "Badly constructed polygons!" 155 message << "Badly constructed polygons!" 149 << "\nNumber of polygons: " << nba 156 << "\nNumber of polygons: " << nbases 150 << "\nPolygon #0 size: " << (*fPol 157 << "\nPolygon #0 size: " << (*fPolygons)[0]->size() 151 << "\nPolygon #1 size: " << (*fPol 158 << "\nPolygon #1 size: " << (*fPolygons)[1]->size() 152 << "\n..."; 159 << "\n..."; 153 G4Exception("G4BoundingEnvelope::CheckBoun << 160 G4Exception("G4BoundingEnvelope::CheckBoundingPolygons()", 154 "GeomMgt0001", FatalException, 161 "GeomMgt0001", FatalException, message); 155 return; 162 return; 156 } 163 } 157 164 158 for (std::size_t k=0; k<nbases; ++k) << 165 for (G4int k=0; k<nbases; ++k) 159 { 166 { 160 std::size_t np = (*fPolygons)[k]->size(); << 167 G4int np = (*fPolygons)[k]->size(); 161 if (np == nsize) continue; 168 if (np == nsize) continue; 162 if (np == 1 && k==0) continue; 169 if (np == 1 && k==0) continue; 163 if (np == 1 && k==nbases-1) continue; 170 if (np == 1 && k==nbases-1) continue; 164 std::ostringstream message; 171 std::ostringstream message; 165 message << "Badly constructed polygons!" 172 message << "Badly constructed polygons!" 166 << "\nNumber of polygons: " << nba 173 << "\nNumber of polygons: " << nbases 167 << "\nPolygon #" << k << " size: " 174 << "\nPolygon #" << k << " size: " << np 168 << "\nexpected size: " << nsize; 175 << "\nexpected size: " << nsize; 169 G4Exception("G4BoundingEnvelope::SetBoundi << 176 G4Exception("G4BoundingEnvelope::SetBoundingPolygons()", 170 "GeomMgt0001", FatalException, 177 "GeomMgt0001", FatalException, message); 171 return; 178 return; 172 } << 179 } 173 } 180 } 174 181 175 ////////////////////////////////////////////// 182 /////////////////////////////////////////////////////////////////////// 176 // 183 // 177 // Quick comparison: bounding box vs voxel, it 184 // Quick comparison: bounding box vs voxel, it return true if further 178 // calculations are not needed 185 // calculations are not needed 179 // 186 // 180 G4bool 187 G4bool 181 G4BoundingEnvelope:: 188 G4BoundingEnvelope:: 182 BoundingBoxVsVoxelLimits(const EAxis pAxis, 189 BoundingBoxVsVoxelLimits(const EAxis pAxis, 183 const G4VoxelLimits& 190 const G4VoxelLimits& pVoxelLimits, 184 const G4Transform3D& 191 const G4Transform3D& pTransform3D, 185 G4double& pMin, G4dou 192 G4double& pMin, G4double& pMax) const 186 { 193 { 187 pMin = kInfinity; 194 pMin = kInfinity; 188 pMax = -kInfinity; 195 pMax = -kInfinity; 189 G4double xminlim = pVoxelLimits.GetMinXExten 196 G4double xminlim = pVoxelLimits.GetMinXExtent(); 190 G4double xmaxlim = pVoxelLimits.GetMaxXExten 197 G4double xmaxlim = pVoxelLimits.GetMaxXExtent(); 191 G4double yminlim = pVoxelLimits.GetMinYExten 198 G4double yminlim = pVoxelLimits.GetMinYExtent(); 192 G4double ymaxlim = pVoxelLimits.GetMaxYExten 199 G4double ymaxlim = pVoxelLimits.GetMaxYExtent(); 193 G4double zminlim = pVoxelLimits.GetMinZExten 200 G4double zminlim = pVoxelLimits.GetMinZExtent(); 194 G4double zmaxlim = pVoxelLimits.GetMaxZExten 201 G4double zmaxlim = pVoxelLimits.GetMaxZExtent(); 195 202 196 // Special case of pure translation 203 // Special case of pure translation 197 // 204 // 198 if (pTransform3D.xx()==1 && pTransform3D.yy( 205 if (pTransform3D.xx()==1 && pTransform3D.yy()==1 && pTransform3D.zz()==1) 199 { 206 { 200 G4double xmin = fMin.x() + pTransform3D.dx 207 G4double xmin = fMin.x() + pTransform3D.dx(); 201 G4double xmax = fMax.x() + pTransform3D.dx 208 G4double xmax = fMax.x() + pTransform3D.dx(); 202 G4double ymin = fMin.y() + pTransform3D.dy 209 G4double ymin = fMin.y() + pTransform3D.dy(); 203 G4double ymax = fMax.y() + pTransform3D.dy 210 G4double ymax = fMax.y() + pTransform3D.dy(); 204 G4double zmin = fMin.z() + pTransform3D.dz 211 G4double zmin = fMin.z() + pTransform3D.dz(); 205 G4double zmax = fMax.z() + pTransform3D.dz 212 G4double zmax = fMax.z() + pTransform3D.dz(); 206 213 207 if (xmin-kCarTolerance > xmaxlim) return t 214 if (xmin-kCarTolerance > xmaxlim) return true; 208 if (xmax+kCarTolerance < xminlim) return t 215 if (xmax+kCarTolerance < xminlim) return true; 209 if (ymin-kCarTolerance > ymaxlim) return t 216 if (ymin-kCarTolerance > ymaxlim) return true; 210 if (ymax+kCarTolerance < yminlim) return t 217 if (ymax+kCarTolerance < yminlim) return true; 211 if (zmin-kCarTolerance > zmaxlim) return t 218 if (zmin-kCarTolerance > zmaxlim) return true; 212 if (zmax+kCarTolerance < zminlim) return t 219 if (zmax+kCarTolerance < zminlim) return true; 213 220 214 if (xmin >= xminlim && xmax <= xmaxlim && 221 if (xmin >= xminlim && xmax <= xmaxlim && 215 ymin >= yminlim && ymax <= ymaxlim && 222 ymin >= yminlim && ymax <= ymaxlim && 216 zmin >= zminlim && zmax <= zmaxlim) 223 zmin >= zminlim && zmax <= zmaxlim) 217 { 224 { 218 if (pAxis == kXAxis) 225 if (pAxis == kXAxis) 219 { 226 { 220 pMin = (xmin-kCarTolerance < xminlim) 227 pMin = (xmin-kCarTolerance < xminlim) ? xminlim : xmin; 221 pMax = (xmax+kCarTolerance > xmaxlim) 228 pMax = (xmax+kCarTolerance > xmaxlim) ? xmaxlim : xmax; 222 } << 229 } 223 else if (pAxis == kYAxis) 230 else if (pAxis == kYAxis) 224 { 231 { 225 pMin = (ymin-kCarTolerance < yminlim) 232 pMin = (ymin-kCarTolerance < yminlim) ? yminlim : ymin; 226 pMax = (ymax+kCarTolerance > ymaxlim) 233 pMax = (ymax+kCarTolerance > ymaxlim) ? ymaxlim : ymax; 227 } 234 } 228 else if (pAxis == kZAxis) 235 else if (pAxis == kZAxis) 229 { 236 { 230 pMin = (zmin-kCarTolerance < zminlim) 237 pMin = (zmin-kCarTolerance < zminlim) ? zminlim : zmin; 231 pMax = (zmax+kCarTolerance > zmaxlim) 238 pMax = (zmax+kCarTolerance > zmaxlim) ? zmaxlim : zmax; 232 } 239 } 233 pMin -= kCarTolerance; 240 pMin -= kCarTolerance; 234 pMax += kCarTolerance; 241 pMax += kCarTolerance; 235 return true; 242 return true; 236 } 243 } 237 } 244 } 238 245 239 // Find max scale factor of the transformati << 246 // Find max scale factor of the transformation, set delta 240 // equal to kCarTolerance multiplied by the 247 // equal to kCarTolerance multiplied by the scale factor 241 // 248 // 242 G4double scale = FindScaleFactor(pTransform3 249 G4double scale = FindScaleFactor(pTransform3D); 243 G4double delta = kCarTolerance*scale; << 250 G4double delta = kCarTolerance*scale; 244 251 245 // Set the sphere surrounding the bounding b 252 // Set the sphere surrounding the bounding box 246 // 253 // 247 G4Point3D center = pTransform3D*G4Point3D(0. 254 G4Point3D center = pTransform3D*G4Point3D(0.5*(fMin+fMax)); 248 G4double radius = 0.5*scale*(fMax-fMin).mag 255 G4double radius = 0.5*scale*(fMax-fMin).mag() + delta; 249 256 250 // Check if the sphere surrounding the bound 257 // Check if the sphere surrounding the bounding box is outside 251 // the voxel limits 258 // the voxel limits 252 // 259 // 253 if (center.x()-radius > xmaxlim) return true 260 if (center.x()-radius > xmaxlim) return true; 254 if (center.y()-radius > ymaxlim) return true 261 if (center.y()-radius > ymaxlim) return true; 255 if (center.z()-radius > zmaxlim) return true 262 if (center.z()-radius > zmaxlim) return true; 256 if (center.x()+radius < xminlim) return true 263 if (center.x()+radius < xminlim) return true; 257 if (center.y()+radius < yminlim) return true 264 if (center.y()+radius < yminlim) return true; 258 if (center.z()+radius < zminlim) return true 265 if (center.z()+radius < zminlim) return true; 259 return false; 266 return false; 260 } 267 } 261 268 262 ////////////////////////////////////////////// 269 /////////////////////////////////////////////////////////////////////// 263 // 270 // 264 // Calculate extent of the specified bounding 271 // Calculate extent of the specified bounding envelope 265 // 272 // 266 G4bool 273 G4bool 267 G4BoundingEnvelope::CalculateExtent(const EAxi 274 G4BoundingEnvelope::CalculateExtent(const EAxis pAxis, 268 const G4Vo 275 const G4VoxelLimits& pVoxelLimits, 269 const G4Tr 276 const G4Transform3D& pTransform3D, 270 G4double& 277 G4double& pMin, G4double& pMax) const 271 { 278 { 272 pMin = kInfinity; 279 pMin = kInfinity; 273 pMax = -kInfinity; 280 pMax = -kInfinity; 274 G4double xminlim = pVoxelLimits.GetMinXExten 281 G4double xminlim = pVoxelLimits.GetMinXExtent(); 275 G4double xmaxlim = pVoxelLimits.GetMaxXExten 282 G4double xmaxlim = pVoxelLimits.GetMaxXExtent(); 276 G4double yminlim = pVoxelLimits.GetMinYExten 283 G4double yminlim = pVoxelLimits.GetMinYExtent(); 277 G4double ymaxlim = pVoxelLimits.GetMaxYExten 284 G4double ymaxlim = pVoxelLimits.GetMaxYExtent(); 278 G4double zminlim = pVoxelLimits.GetMinZExten 285 G4double zminlim = pVoxelLimits.GetMinZExtent(); 279 G4double zmaxlim = pVoxelLimits.GetMaxZExten 286 G4double zmaxlim = pVoxelLimits.GetMaxZExtent(); 280 287 281 // Special case of pure translation 288 // Special case of pure translation 282 // 289 // 283 if (pTransform3D.xx()==1 && pTransform3D.yy( 290 if (pTransform3D.xx()==1 && pTransform3D.yy()==1 && pTransform3D.zz()==1) 284 { 291 { 285 G4double xmin = fMin.x() + pTransform3D.dx 292 G4double xmin = fMin.x() + pTransform3D.dx(); 286 G4double xmax = fMax.x() + pTransform3D.dx 293 G4double xmax = fMax.x() + pTransform3D.dx(); 287 G4double ymin = fMin.y() + pTransform3D.dy 294 G4double ymin = fMin.y() + pTransform3D.dy(); 288 G4double ymax = fMax.y() + pTransform3D.dy 295 G4double ymax = fMax.y() + pTransform3D.dy(); 289 G4double zmin = fMin.z() + pTransform3D.dz 296 G4double zmin = fMin.z() + pTransform3D.dz(); 290 G4double zmax = fMax.z() + pTransform3D.dz 297 G4double zmax = fMax.z() + pTransform3D.dz(); 291 298 292 if (xmin-kCarTolerance > xmaxlim) return f 299 if (xmin-kCarTolerance > xmaxlim) return false; 293 if (xmax+kCarTolerance < xminlim) return f 300 if (xmax+kCarTolerance < xminlim) return false; 294 if (ymin-kCarTolerance > ymaxlim) return f 301 if (ymin-kCarTolerance > ymaxlim) return false; 295 if (ymax+kCarTolerance < yminlim) return f 302 if (ymax+kCarTolerance < yminlim) return false; 296 if (zmin-kCarTolerance > zmaxlim) return f 303 if (zmin-kCarTolerance > zmaxlim) return false; 297 if (zmax+kCarTolerance < zminlim) return f 304 if (zmax+kCarTolerance < zminlim) return false; 298 305 299 if (fPolygons == nullptr) 306 if (fPolygons == nullptr) 300 { 307 { 301 if (pAxis == kXAxis) 308 if (pAxis == kXAxis) 302 { 309 { 303 pMin = (xmin-kCarTolerance < xminlim) 310 pMin = (xmin-kCarTolerance < xminlim) ? xminlim : xmin; 304 pMax = (xmax+kCarTolerance > xmaxlim) 311 pMax = (xmax+kCarTolerance > xmaxlim) ? xmaxlim : xmax; 305 } << 312 } 306 else if (pAxis == kYAxis) 313 else if (pAxis == kYAxis) 307 { 314 { 308 pMin = (ymin-kCarTolerance < yminlim) 315 pMin = (ymin-kCarTolerance < yminlim) ? yminlim : ymin; 309 pMax = (ymax+kCarTolerance > ymaxlim) 316 pMax = (ymax+kCarTolerance > ymaxlim) ? ymaxlim : ymax; 310 } 317 } 311 else if (pAxis == kZAxis) 318 else if (pAxis == kZAxis) 312 { 319 { 313 pMin = (zmin-kCarTolerance < zminlim) 320 pMin = (zmin-kCarTolerance < zminlim) ? zminlim : zmin; 314 pMax = (zmax+kCarTolerance > zmaxlim) 321 pMax = (zmax+kCarTolerance > zmaxlim) ? zmaxlim : zmax; 315 } 322 } 316 pMin -= kCarTolerance; 323 pMin -= kCarTolerance; 317 pMax += kCarTolerance; 324 pMax += kCarTolerance; 318 return true; 325 return true; 319 } 326 } 320 } 327 } 321 328 322 // Find max scale factor of the transformati << 329 // Find max scale factor of the transformation, set delta 323 // equal to kCarTolerance multiplied by the 330 // equal to kCarTolerance multiplied by the scale factor 324 // 331 // 325 G4double scale = FindScaleFactor(pTransform3 332 G4double scale = FindScaleFactor(pTransform3D); 326 G4double delta = kCarTolerance*scale; << 333 G4double delta = kCarTolerance*scale; 327 334 328 // Set the sphere surrounding the bounding b 335 // Set the sphere surrounding the bounding box 329 // 336 // 330 G4Point3D center = pTransform3D*G4Point3D(0. 337 G4Point3D center = pTransform3D*G4Point3D(0.5*(fMin+fMax)); 331 G4double radius = 0.5*scale*(fMax-fMin).mag 338 G4double radius = 0.5*scale*(fMax-fMin).mag() + delta; 332 339 333 // Check if the sphere surrounding the bound 340 // Check if the sphere surrounding the bounding box is within 334 // the voxel limits, if so then transform on 341 // the voxel limits, if so then transform only one coordinate 335 // 342 // 336 if (center.x()-radius >= xminlim && center.x 343 if (center.x()-radius >= xminlim && center.x()+radius <= xmaxlim && 337 center.y()-radius >= yminlim && center.y 344 center.y()-radius >= yminlim && center.y()+radius <= ymaxlim && 338 center.z()-radius >= zminlim && center.z 345 center.z()-radius >= zminlim && center.z()+radius <= zmaxlim ) 339 { 346 { 340 G4double cx, cy, cz, cd; 347 G4double cx, cy, cz, cd; 341 if (pAxis == kXAxis) 348 if (pAxis == kXAxis) 342 { << 349 { 343 cx = pTransform3D.xx(); << 350 cx = pTransform3D.xx(); 344 cy = pTransform3D.xy(); 351 cy = pTransform3D.xy(); 345 cz = pTransform3D.xz(); 352 cz = pTransform3D.xz(); 346 cd = pTransform3D.dx(); 353 cd = pTransform3D.dx(); 347 } 354 } 348 else if (pAxis == kYAxis) 355 else if (pAxis == kYAxis) 349 { << 356 { 350 cx = pTransform3D.yx(); << 357 cx = pTransform3D.yx(); 351 cy = pTransform3D.yy(); 358 cy = pTransform3D.yy(); 352 cz = pTransform3D.yz(); 359 cz = pTransform3D.yz(); 353 cd = pTransform3D.dy(); 360 cd = pTransform3D.dy(); 354 } 361 } 355 else if (pAxis == kZAxis) 362 else if (pAxis == kZAxis) 356 { << 363 { 357 cx = pTransform3D.zx(); << 364 cx = pTransform3D.zx(); 358 cy = pTransform3D.zy(); 365 cy = pTransform3D.zy(); 359 cz = pTransform3D.zz(); 366 cz = pTransform3D.zz(); 360 cd = pTransform3D.dz(); 367 cd = pTransform3D.dz(); 361 } 368 } 362 else 369 else 363 { << 370 { 364 cx = cy = cz = cd = kInfinity; 371 cx = cy = cz = cd = kInfinity; 365 } 372 } 366 G4double emin = kInfinity, emax = -kInfini 373 G4double emin = kInfinity, emax = -kInfinity; 367 if (fPolygons == nullptr) 374 if (fPolygons == nullptr) 368 { 375 { 369 G4double coor; 376 G4double coor; 370 coor = cx*fMin.x() + cy*fMin.y() + cz*fM 377 coor = cx*fMin.x() + cy*fMin.y() + cz*fMin.z() + cd; 371 if (coor < emin) emin = coor; 378 if (coor < emin) emin = coor; 372 if (coor > emax) emax = coor; 379 if (coor > emax) emax = coor; 373 coor = cx*fMax.x() + cy*fMin.y() + cz*fM 380 coor = cx*fMax.x() + cy*fMin.y() + cz*fMin.z() + cd; 374 if (coor < emin) emin = coor; 381 if (coor < emin) emin = coor; 375 if (coor > emax) emax = coor; 382 if (coor > emax) emax = coor; 376 coor = cx*fMax.x() + cy*fMax.y() + cz*fM 383 coor = cx*fMax.x() + cy*fMax.y() + cz*fMin.z() + cd; 377 if (coor < emin) emin = coor; 384 if (coor < emin) emin = coor; 378 if (coor > emax) emax = coor; 385 if (coor > emax) emax = coor; 379 coor = cx*fMin.x() + cy*fMax.y() + cz*fM 386 coor = cx*fMin.x() + cy*fMax.y() + cz*fMin.z() + cd; 380 if (coor < emin) emin = coor; 387 if (coor < emin) emin = coor; 381 if (coor > emax) emax = coor; 388 if (coor > emax) emax = coor; 382 coor = cx*fMin.x() + cy*fMin.y() + cz*fM 389 coor = cx*fMin.x() + cy*fMin.y() + cz*fMax.z() + cd; 383 if (coor < emin) emin = coor; 390 if (coor < emin) emin = coor; 384 if (coor > emax) emax = coor; 391 if (coor > emax) emax = coor; 385 coor = cx*fMax.x() + cy*fMin.y() + cz*fM 392 coor = cx*fMax.x() + cy*fMin.y() + cz*fMax.z() + cd; 386 if (coor < emin) emin = coor; 393 if (coor < emin) emin = coor; 387 if (coor > emax) emax = coor; 394 if (coor > emax) emax = coor; 388 coor = cx*fMax.x() + cy*fMax.y() + cz*fM 395 coor = cx*fMax.x() + cy*fMax.y() + cz*fMax.z() + cd; 389 if (coor < emin) emin = coor; 396 if (coor < emin) emin = coor; 390 if (coor > emax) emax = coor; 397 if (coor > emax) emax = coor; 391 coor = cx*fMin.x() + cy*fMax.y() + cz*fM 398 coor = cx*fMin.x() + cy*fMax.y() + cz*fMax.z() + cd; 392 if (coor < emin) emin = coor; 399 if (coor < emin) emin = coor; 393 if (coor > emax) emax = coor; 400 if (coor > emax) emax = coor; 394 } 401 } 395 else 402 else 396 { 403 { 397 for (const auto & polygon : *fPolygons) << 404 for (auto ibase=fPolygons->cbegin(); ibase!=fPolygons->cend(); ++ibase) 398 { << 405 { 399 for (auto ipoint=polygon->cbegin(); ip << 406 for (auto ipoint=(*ibase)->cbegin(); ipoint!=(*ibase)->cend(); ++ipoint) 400 { 407 { 401 G4double coor = ipoint->x()*cx + ipo 408 G4double coor = ipoint->x()*cx + ipoint->y()*cy + ipoint->z()*cz + cd; 402 if (coor < emin) emin = coor; 409 if (coor < emin) emin = coor; 403 if (coor > emax) emax = coor; 410 if (coor > emax) emax = coor; 404 } 411 } 405 } 412 } 406 } 413 } 407 pMin = emin - delta; 414 pMin = emin - delta; 408 pMax = emax + delta; 415 pMax = emax + delta; 409 return true; 416 return true; 410 } << 417 } 411 418 412 // Check if the sphere surrounding the bound 419 // Check if the sphere surrounding the bounding box is outside 413 // the voxel limits 420 // the voxel limits 414 // 421 // 415 if (center.x()-radius > xmaxlim) return fals 422 if (center.x()-radius > xmaxlim) return false; 416 if (center.y()-radius > ymaxlim) return fals 423 if (center.y()-radius > ymaxlim) return false; 417 if (center.z()-radius > zmaxlim) return fals 424 if (center.z()-radius > zmaxlim) return false; 418 if (center.x()+radius < xminlim) return fals 425 if (center.x()+radius < xminlim) return false; 419 if (center.y()+radius < yminlim) return fals 426 if (center.y()+radius < yminlim) return false; 420 if (center.z()+radius < zminlim) return fals 427 if (center.z()+radius < zminlim) return false; 421 428 422 // Transform polygons << 429 // Allocate memory for transformed polygons >> 430 // >> 431 G4int nbases = (fPolygons == 0) ? 2 : fPolygons->size(); >> 432 std::vector<G4Polygon3D*> bases(nbases); >> 433 if (fPolygons == nullptr) >> 434 { >> 435 bases[0] = new G4Polygon3D(4); >> 436 bases[1] = new G4Polygon3D(4); >> 437 } >> 438 else >> 439 { >> 440 for (G4int i=0; i<nbases; ++i) >> 441 { >> 442 bases[i] = new G4Polygon3D((*fPolygons)[i]->size()); >> 443 } >> 444 } >> 445 >> 446 // Transform vertices 423 // 447 // 424 std::vector<G4Point3D> vertices; << 448 TransformVertices(pTransform3D, bases); 425 std::vector<std::pair<G4int, G4int>> bases; << 426 TransformVertices(pTransform3D, vertices, ba << 427 std::size_t nbases = bases.size(); << 428 449 429 // Create adjusted G4VoxelLimits box. New li << 450 // Create adjusted G4VoxelLimits box. New limits are extended by 430 // delta, kCarTolerance multiplied by max sc 451 // delta, kCarTolerance multiplied by max scale factor of 431 // the transformation 452 // the transformation 432 // 453 // 433 EAxis axes[] = { kXAxis, kYAxis, kZAxis }; << 454 EAxis axis[] = { kXAxis,kYAxis,kZAxis }; 434 G4VoxelLimits limits; // default is unlimite << 455 G4VoxelLimits limits; // default is unlimited 435 for (const auto & iAxis : axes) << 456 for (auto i=0; i<3; ++i) 436 { << 457 { 437 if (pVoxelLimits.IsLimited(iAxis)) << 458 if (pVoxelLimits.IsLimited(axis[i])) 438 { << 459 { 439 G4double emin = pVoxelLimits.GetMinExten << 460 G4double emin = pVoxelLimits.GetMinExtent(axis[i]) - delta; 440 G4double emax = pVoxelLimits.GetMaxExten << 461 G4double emax = pVoxelLimits.GetMaxExtent(axis[i]) + delta; 441 limits.AddLimit(iAxis, emin, emax); << 462 limits.AddLimit(axis[i], emin, emax); 442 } 463 } 443 } 464 } 444 465 445 // Main loop along the set of prisms 466 // Main loop along the set of prisms 446 // 467 // 447 G4Polygon3D baseA, baseB; << 448 G4Segment3D extent; 468 G4Segment3D extent; 449 extent.first = G4Point3D( kInfinity, kInfin 469 extent.first = G4Point3D( kInfinity, kInfinity, kInfinity); 450 extent.second = G4Point3D(-kInfinity,-kInfin 470 extent.second = G4Point3D(-kInfinity,-kInfinity,-kInfinity); 451 for (std::size_t k=0; k<nbases-1; ++k) << 471 for (G4int k=0; k<nbases-1; ++k) 452 { 472 { 453 baseA.resize(bases[k].second); << 454 for (G4int i = 0; i < bases[k].second; ++i << 455 baseA[i] = vertices[bases[k].first + i]; << 456 << 457 baseB.resize(bases[k+1].second); << 458 for (G4int i = 0; i < bases[k+1].second; + << 459 baseB[i] = vertices[bases[k+1].first + i << 460 << 461 // Find bounding box of current prism 473 // Find bounding box of current prism >> 474 G4Polygon3D* baseA = bases[k]; >> 475 G4Polygon3D* baseB = bases[k+1]; 462 G4Segment3D prismAABB; 476 G4Segment3D prismAABB; 463 GetPrismAABB(baseA, baseB, prismAABB); << 477 GetPrismAABB(*baseA, *baseB, prismAABB); 464 478 465 // Check if prismAABB is completely within 479 // Check if prismAABB is completely within the voxel limits 466 if (prismAABB.first.x() >= limits.GetMinXE 480 if (prismAABB.first.x() >= limits.GetMinXExtent() && 467 prismAABB.first.y() >= limits.GetMinYE 481 prismAABB.first.y() >= limits.GetMinYExtent() && 468 prismAABB.first.z() >= limits.GetMinZE 482 prismAABB.first.z() >= limits.GetMinZExtent() && 469 prismAABB.second.x()<= limits.GetMaxXE 483 prismAABB.second.x()<= limits.GetMaxXExtent() && 470 prismAABB.second.y()<= limits.GetMaxYE 484 prismAABB.second.y()<= limits.GetMaxYExtent() && 471 prismAABB.second.z()<= limits.GetMaxZE << 485 prismAABB.second.z()<= limits.GetMaxZExtent()) 472 { 486 { 473 if (extent.first.x() > prismAABB.first. 487 if (extent.first.x() > prismAABB.first.x()) 474 extent.first.setX( prismAABB.first.x() << 488 extent.first.setX( prismAABB.first.x() ); 475 if (extent.first.y() > prismAABB.first. 489 if (extent.first.y() > prismAABB.first.y()) 476 extent.first.setY( prismAABB.first.y() << 490 extent.first.setY( prismAABB.first.y() ); 477 if (extent.first.z() > prismAABB.first. 491 if (extent.first.z() > prismAABB.first.z()) 478 extent.first.setZ( prismAABB.first.z() << 492 extent.first.setZ( prismAABB.first.z() ); 479 if (extent.second.x() < prismAABB.second 493 if (extent.second.x() < prismAABB.second.x()) 480 extent.second.setX(prismAABB.second.x( << 494 extent.second.setX(prismAABB.second.x()); 481 if (extent.second.y() < prismAABB.second 495 if (extent.second.y() < prismAABB.second.y()) 482 extent.second.setY(prismAABB.second.y( << 496 extent.second.setY(prismAABB.second.y()); 483 if (extent.second.z() < prismAABB.second 497 if (extent.second.z() < prismAABB.second.z()) 484 extent.second.setZ(prismAABB.second.z( 498 extent.second.setZ(prismAABB.second.z()); 485 continue; << 499 continue; 486 } 500 } 487 501 488 // Check if prismAABB is outside the voxel 502 // Check if prismAABB is outside the voxel limits 489 if (prismAABB.first.x() > limits.GetMaxXE 503 if (prismAABB.first.x() > limits.GetMaxXExtent()) continue; 490 if (prismAABB.first.y() > limits.GetMaxYE 504 if (prismAABB.first.y() > limits.GetMaxYExtent()) continue; 491 if (prismAABB.first.z() > limits.GetMaxZE 505 if (prismAABB.first.z() > limits.GetMaxZExtent()) continue; 492 if (prismAABB.second.x() < limits.GetMinXE 506 if (prismAABB.second.x() < limits.GetMinXExtent()) continue; 493 if (prismAABB.second.y() < limits.GetMinYE 507 if (prismAABB.second.y() < limits.GetMinYExtent()) continue; 494 if (prismAABB.second.z() < limits.GetMinZE 508 if (prismAABB.second.z() < limits.GetMinZExtent()) continue; 495 509 496 // Clip edges of the prism by adjusted G4V 510 // Clip edges of the prism by adjusted G4VoxelLimits box 497 std::vector<G4Segment3D> vecEdges; << 511 std::vector<G4Segment3D> vecEdges; 498 CreateListOfEdges(baseA, baseB, vecEdges); << 512 CreateListOfEdges(*baseA, *baseB, vecEdges); 499 if (ClipEdgesByVoxel(vecEdges, limits, ext 513 if (ClipEdgesByVoxel(vecEdges, limits, extent)) continue; 500 514 501 // Some edges of the prism are completely 515 // Some edges of the prism are completely outside of the voxel 502 // limits, clip selected edges (see bits) 516 // limits, clip selected edges (see bits) of adjusted G4VoxelLimits 503 // by the prism << 517 // by the prism 504 G4int bits = 0x000; 518 G4int bits = 0x000; 505 if (limits.GetMinXExtent() < prismAABB.fir 519 if (limits.GetMinXExtent() < prismAABB.first.x()) 506 bits |= 0x988; // 1001 1000 1000 520 bits |= 0x988; // 1001 1000 1000 507 if (limits.GetMaxXExtent() > prismAABB.sec 521 if (limits.GetMaxXExtent() > prismAABB.second.x()) 508 bits |= 0x622; // 0110 0010 0010 522 bits |= 0x622; // 0110 0010 0010 509 523 510 if (limits.GetMinYExtent() < prismAABB.fir 524 if (limits.GetMinYExtent() < prismAABB.first.y()) 511 bits |= 0x311; // 0011 0001 0001 525 bits |= 0x311; // 0011 0001 0001 512 if (limits.GetMaxYExtent() > prismAABB.sec 526 if (limits.GetMaxYExtent() > prismAABB.second.y()) 513 bits |= 0xC44; // 1100 0100 0100 527 bits |= 0xC44; // 1100 0100 0100 514 528 515 if (limits.GetMinZExtent() < prismAABB.fir 529 if (limits.GetMinZExtent() < prismAABB.first.z()) 516 bits |= 0x00F; // 0000 0000 1111 530 bits |= 0x00F; // 0000 0000 1111 517 if (limits.GetMaxZExtent() > prismAABB.sec 531 if (limits.GetMaxZExtent() > prismAABB.second.z()) 518 bits |= 0x0F0; // 0000 1111 0000 532 bits |= 0x0F0; // 0000 1111 0000 519 if (bits == 0xFFF) continue; 533 if (bits == 0xFFF) continue; 520 534 521 std::vector<G4Plane3D> vecPlanes; << 535 std::vector<G4Plane3D> vecPlanes; 522 CreateListOfPlanes(baseA, baseB, vecPlanes << 536 CreateListOfPlanes(*baseA, *baseB, vecPlanes); 523 ClipVoxelByPlanes(bits, limits, vecPlanes, 537 ClipVoxelByPlanes(bits, limits, vecPlanes, prismAABB, extent); 524 } // End of the main loop 538 } // End of the main loop 525 539 526 // Final adjustment of the extent << 540 // Free memory 527 // 541 // >> 542 for (G4int i=0; i<nbases; ++i) { delete bases[i]; bases[i] = 0; } >> 543 >> 544 // Final adjustment of the extent >> 545 // 528 G4double emin = 0, emax = 0; 546 G4double emin = 0, emax = 0; 529 if (pAxis == kXAxis) { emin = extent.first.x 547 if (pAxis == kXAxis) { emin = extent.first.x(); emax = extent.second.x(); } 530 if (pAxis == kYAxis) { emin = extent.first.y 548 if (pAxis == kYAxis) { emin = extent.first.y(); emax = extent.second.y(); } 531 if (pAxis == kZAxis) { emin = extent.first.z 549 if (pAxis == kZAxis) { emin = extent.first.z(); emax = extent.second.z(); } 532 550 533 if (emin > emax) return false; 551 if (emin > emax) return false; 534 emin -= delta; 552 emin -= delta; 535 emax += delta; 553 emax += delta; 536 G4double minlim = pVoxelLimits.GetMinExtent( 554 G4double minlim = pVoxelLimits.GetMinExtent(pAxis); 537 G4double maxlim = pVoxelLimits.GetMaxExtent( 555 G4double maxlim = pVoxelLimits.GetMaxExtent(pAxis); 538 pMin = (emin < minlim) ? minlim-kCarToleranc 556 pMin = (emin < minlim) ? minlim-kCarTolerance : emin; 539 pMax = (emax > maxlim) ? maxlim+kCarToleranc 557 pMax = (emax > maxlim) ? maxlim+kCarTolerance : emax; 540 return true; 558 return true; 541 } 559 } 542 560 >> 561 543 ////////////////////////////////////////////// 562 /////////////////////////////////////////////////////////////////////// 544 // 563 // 545 // Find max scale factor of the transformation 564 // Find max scale factor of the transformation 546 // 565 // 547 G4double 566 G4double 548 G4BoundingEnvelope::FindScaleFactor(const G4Tr 567 G4BoundingEnvelope::FindScaleFactor(const G4Transform3D& pTransform3D) const 549 { 568 { 550 if (pTransform3D.xx() == 1. && 569 if (pTransform3D.xx() == 1. && 551 pTransform3D.yy() == 1. && 570 pTransform3D.yy() == 1. && 552 pTransform3D.zz() == 1.) return 1.; 571 pTransform3D.zz() == 1.) return 1.; 553 572 554 G4double xx = pTransform3D.xx(); 573 G4double xx = pTransform3D.xx(); 555 G4double yx = pTransform3D.yx(); 574 G4double yx = pTransform3D.yx(); 556 G4double zx = pTransform3D.zx(); 575 G4double zx = pTransform3D.zx(); 557 G4double sxsx = xx*xx + yx*yx + zx*zx; 576 G4double sxsx = xx*xx + yx*yx + zx*zx; 558 G4double xy = pTransform3D.xy(); 577 G4double xy = pTransform3D.xy(); 559 G4double yy = pTransform3D.yy(); 578 G4double yy = pTransform3D.yy(); 560 G4double zy = pTransform3D.zy(); 579 G4double zy = pTransform3D.zy(); 561 G4double sysy = xy*xy + yy*yy + zy*zy; 580 G4double sysy = xy*xy + yy*yy + zy*zy; 562 G4double xz = pTransform3D.xz(); 581 G4double xz = pTransform3D.xz(); 563 G4double yz = pTransform3D.yz(); 582 G4double yz = pTransform3D.yz(); 564 G4double zz = pTransform3D.zz(); 583 G4double zz = pTransform3D.zz(); 565 G4double szsz = xz*xz + yz*yz + zz*zz; 584 G4double szsz = xz*xz + yz*yz + zz*zz; 566 G4double ss = std::max(std::max(sxsx,sysy),s 585 G4double ss = std::max(std::max(sxsx,sysy),szsz); 567 return (ss <= 1.) ? 1. : std::sqrt(ss); 586 return (ss <= 1.) ? 1. : std::sqrt(ss); 568 } 587 } 569 588 570 ////////////////////////////////////////////// 589 /////////////////////////////////////////////////////////////////////// 571 // 590 // 572 // Transform polygonal bases 591 // Transform polygonal bases 573 // 592 // 574 void 593 void 575 G4BoundingEnvelope:: << 594 G4BoundingEnvelope::TransformVertices(const G4Transform3D& pTransform3D, 576 TransformVertices(const G4Transform3D& pTransf << 595 std::vector<G4Polygon3D*>& pBases) const 577 std::vector<G4Point3D>& pVer << 578 std::vector<std::pair<G4int, << 579 { 596 { 580 G4ThreeVectorList baseA(4), baseB(4); 597 G4ThreeVectorList baseA(4), baseB(4); 581 std::vector<const G4ThreeVectorList*> aabb(2 598 std::vector<const G4ThreeVectorList*> aabb(2); 582 aabb[0] = &baseA; << 599 aabb[0] = &baseA; aabb[1] = &baseB; 583 aabb[1] = &baseB; << 584 if (fPolygons == nullptr) 600 if (fPolygons == nullptr) 585 { 601 { 586 baseA[0].set(fMin.x(),fMin.y(),fMin.z()); 602 baseA[0].set(fMin.x(),fMin.y(),fMin.z()); 587 baseA[1].set(fMax.x(),fMin.y(),fMin.z()); 603 baseA[1].set(fMax.x(),fMin.y(),fMin.z()); 588 baseA[2].set(fMax.x(),fMax.y(),fMin.z()); 604 baseA[2].set(fMax.x(),fMax.y(),fMin.z()); 589 baseA[3].set(fMin.x(),fMax.y(),fMin.z()); 605 baseA[3].set(fMin.x(),fMax.y(),fMin.z()); 590 baseB[0].set(fMin.x(),fMin.y(),fMax.z()); 606 baseB[0].set(fMin.x(),fMin.y(),fMax.z()); 591 baseB[1].set(fMax.x(),fMin.y(),fMax.z()); 607 baseB[1].set(fMax.x(),fMin.y(),fMax.z()); 592 baseB[2].set(fMax.x(),fMax.y(),fMax.z()); 608 baseB[2].set(fMax.x(),fMax.y(),fMax.z()); 593 baseB[3].set(fMin.x(),fMax.y(),fMax.z()); 609 baseB[3].set(fMin.x(),fMax.y(),fMax.z()); 594 } 610 } 595 auto ia = (fPolygons == nullptr) ? aabb.c << 611 std::vector<const G4ThreeVectorList*>::const_iterator ia, iaend; 596 auto iaend = (fPolygons == nullptr) ? aabb.c << 612 std::vector<G4Polygon3D*>::iterator ib = pBases.begin(); 597 << 613 ia = (fPolygons == 0) ? aabb.begin() : fPolygons->begin(); 598 // Fill vector of bases << 614 iaend = (fPolygons == 0) ? aabb.end() : fPolygons->end(); 599 // << 600 G4int index = 0; << 601 for (auto i = ia; i != iaend; ++i) << 602 { << 603 auto nv = (G4int)(*i)->size(); << 604 pBases.emplace_back(index, nv); << 605 index += nv; << 606 } << 607 615 608 // Fill vector of transformed vertices << 616 if (pTransform3D.xx()==1 && pTransform3D.yy()==1 && pTransform3D.zz()==1) 609 // << 610 if (pTransform3D.xx() == 1. && << 611 pTransform3D.yy() == 1. && << 612 pTransform3D.zz() == 1.) << 613 { 617 { 614 G4ThreeVector offset = pTransform3D.getTra 618 G4ThreeVector offset = pTransform3D.getTranslation(); 615 for (auto i = ia; i != iaend; ++i) << 619 for ( ; ia != iaend; ++ia, ++ib) 616 for (auto k = (*i)->cbegin(); k != (*i)- << 620 { 617 pVertices.emplace_back((*k) + offset); << 621 G4ThreeVectorList::const_iterator ka = (*ia)->begin(); >> 622 G4Polygon3D::iterator kb = (*ib)->begin(); >> 623 for ( ; ka != (*ia)->end(); ++ka, ++kb) { (*kb) = (*ka) + offset; } >> 624 } 618 } 625 } 619 else 626 else 620 { 627 { 621 for (auto i = ia; i != iaend; ++i) << 628 for ( ; ia != iaend; ++ia, ++ib) 622 for (auto k = (*i)->cbegin(); k != (*i)- << 629 { 623 pVertices.push_back(pTransform3D*G4Poi << 630 G4ThreeVectorList::const_iterator ka = (*ia)->begin(); >> 631 G4Polygon3D::iterator kb = (*ib)->begin(); >> 632 for ( ; ka != (*ia)->end(); ++ka, ++kb) >> 633 { >> 634 (*kb) = pTransform3D*G4Point3D(*ka); >> 635 } >> 636 } 624 } 637 } 625 } 638 } 626 639 627 ////////////////////////////////////////////// 640 /////////////////////////////////////////////////////////////////////// 628 // 641 // 629 // Find bounding box of a prism 642 // Find bounding box of a prism 630 // 643 // 631 void 644 void 632 G4BoundingEnvelope::GetPrismAABB(const G4Polyg 645 G4BoundingEnvelope::GetPrismAABB(const G4Polygon3D& pBaseA, 633 const G4Polyg 646 const G4Polygon3D& pBaseB, 634 G4Segme 647 G4Segment3D& pAABB) const 635 { 648 { 636 G4double xmin = kInfinity, ymin = kInfinit 649 G4double xmin = kInfinity, ymin = kInfinity, zmin = kInfinity; 637 G4double xmax = -kInfinity, ymax = -kInfinit 650 G4double xmax = -kInfinity, ymax = -kInfinity, zmax = -kInfinity; 638 651 639 // First base 652 // First base 640 // 653 // 641 for (const auto & it1 : pBaseA) << 654 for (auto it1 = pBaseA.cbegin(); it1 != pBaseA.cend(); ++it1) 642 { << 655 { 643 G4double x = it1.x(); << 656 G4double x = it1->x(); 644 if (x < xmin) xmin = x; 657 if (x < xmin) xmin = x; 645 if (x > xmax) xmax = x; 658 if (x > xmax) xmax = x; 646 G4double y = it1.y(); << 659 G4double y = it1->y(); 647 if (y < ymin) ymin = y; 660 if (y < ymin) ymin = y; 648 if (y > ymax) ymax = y; 661 if (y > ymax) ymax = y; 649 G4double z = it1.z(); << 662 G4double z = it1->z(); 650 if (z < zmin) zmin = z; 663 if (z < zmin) zmin = z; 651 if (z > zmax) zmax = z; 664 if (z > zmax) zmax = z; 652 } 665 } 653 666 654 // Second base 667 // Second base 655 // 668 // 656 for (const auto & it2 : pBaseB) << 669 for (auto it2 = pBaseB.cbegin(); it2 != pBaseB.cend(); ++it2) 657 { << 670 { 658 G4double x = it2.x(); << 671 G4double x = it2->x(); 659 if (x < xmin) xmin = x; 672 if (x < xmin) xmin = x; 660 if (x > xmax) xmax = x; 673 if (x > xmax) xmax = x; 661 G4double y = it2.y(); << 674 G4double y = it2->y(); 662 if (y < ymin) ymin = y; 675 if (y < ymin) ymin = y; 663 if (y > ymax) ymax = y; 676 if (y > ymax) ymax = y; 664 G4double z = it2.z(); << 677 G4double z = it2->z(); 665 if (z < zmin) zmin = z; 678 if (z < zmin) zmin = z; 666 if (z > zmax) zmax = z; 679 if (z > zmax) zmax = z; 667 } 680 } 668 681 669 // Set bounding box 682 // Set bounding box 670 // 683 // 671 pAABB.first = G4Point3D(xmin,ymin,zmin); 684 pAABB.first = G4Point3D(xmin,ymin,zmin); 672 pAABB.second = G4Point3D(xmax,ymax,zmax); 685 pAABB.second = G4Point3D(xmax,ymax,zmax); 673 } 686 } 674 687 675 ////////////////////////////////////////////// 688 /////////////////////////////////////////////////////////////////////// 676 // 689 // 677 // Create list of edges of a prism 690 // Create list of edges of a prism 678 // 691 // 679 void 692 void 680 G4BoundingEnvelope::CreateListOfEdges(const G4 693 G4BoundingEnvelope::CreateListOfEdges(const G4Polygon3D& baseA, 681 const G4 694 const G4Polygon3D& baseB, 682 std::vec 695 std::vector<G4Segment3D>& pEdges) const 683 { 696 { 684 std::size_t na = baseA.size(); << 697 G4int na = baseA.size(); 685 std::size_t nb = baseB.size(); << 698 G4int nb = baseB.size(); 686 pEdges.clear(); 699 pEdges.clear(); 687 if (na == nb) 700 if (na == nb) 688 { 701 { 689 pEdges.resize(3*na); 702 pEdges.resize(3*na); 690 std::size_t k = na - 1; << 703 G4int k = na - 1; 691 for (std::size_t i=0; i<na; ++i) << 704 for (G4int i=0; i<na; ++i) 692 { 705 { 693 pEdges.emplace_back(baseA[i],baseB[i]); << 706 pEdges.push_back(G4Segment3D(baseA[i],baseB[i])); 694 pEdges.emplace_back(baseA[i],baseA[k]); << 707 pEdges.push_back(G4Segment3D(baseA[i],baseA[k])); 695 pEdges.emplace_back(baseB[i],baseB[k]); << 708 pEdges.push_back(G4Segment3D(baseB[i],baseB[k])); 696 k = i; 709 k = i; 697 } 710 } 698 } 711 } 699 else if (nb == 1) 712 else if (nb == 1) 700 { 713 { 701 pEdges.resize(2*na); 714 pEdges.resize(2*na); 702 std::size_t k = na - 1; << 715 G4int k = na - 1; 703 for (std::size_t i=0; i<na; ++i) << 716 for (G4int i=0; i<na; ++i) 704 { 717 { 705 pEdges.emplace_back(baseA[i],baseA[k]); << 718 pEdges.push_back(G4Segment3D(baseA[i],baseA[k])); 706 pEdges.emplace_back(baseA[i],baseB[0]); << 719 pEdges.push_back(G4Segment3D(baseA[i],baseB[0])); 707 k = i; 720 k = i; 708 } 721 } 709 } 722 } 710 else if (na == 1) 723 else if (na == 1) 711 { 724 { 712 pEdges.resize(2*nb); 725 pEdges.resize(2*nb); 713 std::size_t k = nb - 1; << 726 G4int k = nb - 1; 714 for (std::size_t i=0; i<nb; ++i) << 727 for (G4int i=0; i<nb; ++i) 715 { 728 { 716 pEdges.emplace_back(baseB[i],baseB[k]); << 729 pEdges.push_back(G4Segment3D(baseB[i],baseB[k])); 717 pEdges.emplace_back(baseB[i],baseA[0]); << 730 pEdges.push_back(G4Segment3D(baseB[i],baseA[0])); 718 k = i; 731 k = i; 719 } 732 } 720 } 733 } 721 } 734 } 722 735 723 ////////////////////////////////////////////// 736 /////////////////////////////////////////////////////////////////////// 724 // 737 // 725 // Create list of planes bounding a prism 738 // Create list of planes bounding a prism 726 // 739 // 727 void 740 void 728 G4BoundingEnvelope::CreateListOfPlanes(const G 741 G4BoundingEnvelope::CreateListOfPlanes(const G4Polygon3D& baseA, 729 const G 742 const G4Polygon3D& baseB, 730 std::ve 743 std::vector<G4Plane3D>& pPlanes) const 731 { 744 { 732 // Find centers of the bases and internal po 745 // Find centers of the bases and internal point of the prism 733 // 746 // 734 std::size_t na = baseA.size(); << 747 G4int na = baseA.size(); 735 std::size_t nb = baseB.size(); << 748 G4int nb = baseB.size(); 736 G4Point3D pa(0.,0.,0.), pb(0.,0.,0.), p0; 749 G4Point3D pa(0.,0.,0.), pb(0.,0.,0.), p0; 737 G4Normal3D norm; 750 G4Normal3D norm; 738 for (std::size_t i=0; i<na; ++i) pa += baseA << 751 for (G4int i=0; i<na; ++i) pa += baseA[i]; 739 for (std::size_t i=0; i<nb; ++i) pb += baseB << 752 for (G4int i=0; i<nb; ++i) pb += baseB[i]; 740 pa /= na; pb /= nb; p0 = (pa+pb)/2.; 753 pa /= na; pb /= nb; p0 = (pa+pb)/2.; 741 754 742 // Create list of planes 755 // Create list of planes 743 // 756 // 744 pPlanes.clear(); 757 pPlanes.clear(); 745 if (na == nb) // bases with equal number of 758 if (na == nb) // bases with equal number of vertices 746 { 759 { 747 std::size_t k = na - 1; << 760 G4int k = na - 1; 748 for (std::size_t i=0; i<na; ++i) << 761 for (G4int i=0; i<na; ++i) 749 { 762 { 750 norm = (baseB[k]-baseA[i]).cross(baseA[k 763 norm = (baseB[k]-baseA[i]).cross(baseA[k]-baseB[i]); 751 if (norm.mag2() > kCarTolerance) 764 if (norm.mag2() > kCarTolerance) 752 { 765 { 753 pPlanes.emplace_back(norm,baseA[i]); << 766 pPlanes.push_back(G4Plane3D(norm,baseA[i])); 754 } 767 } 755 k = i; 768 k = i; 756 } 769 } 757 norm = (baseA[2]-baseA[0]).cross(baseA[1]- 770 norm = (baseA[2]-baseA[0]).cross(baseA[1]-pa); 758 if (norm.mag2() > kCarTolerance) 771 if (norm.mag2() > kCarTolerance) 759 { 772 { 760 pPlanes.emplace_back(norm,pa); << 773 pPlanes.push_back(G4Plane3D(norm,pa)); 761 } 774 } 762 norm = (baseB[2]-baseB[0]).cross(baseB[1]- 775 norm = (baseB[2]-baseB[0]).cross(baseB[1]-pb); 763 if (norm.mag2() > kCarTolerance) 776 if (norm.mag2() > kCarTolerance) 764 { 777 { 765 pPlanes.emplace_back(norm,pb); << 778 pPlanes.push_back(G4Plane3D(norm,pb)); 766 } 779 } 767 } 780 } 768 else if (nb == 1) // baseB has one vertex 781 else if (nb == 1) // baseB has one vertex 769 { 782 { 770 std::size_t k = na - 1; << 783 G4int k = na - 1; 771 for (std::size_t i=0; i<na; ++i) << 784 for (G4int i=0; i<na; ++i) 772 { 785 { 773 norm = (baseA[i]-baseB[0]).cross(baseA[k 786 norm = (baseA[i]-baseB[0]).cross(baseA[k]-baseB[0]); 774 if (norm.mag2() > kCarTolerance) 787 if (norm.mag2() > kCarTolerance) 775 { 788 { 776 pPlanes.emplace_back(norm,baseB[0]); << 789 pPlanes.push_back(G4Plane3D(norm,baseB[0])); 777 } 790 } 778 k = i; 791 k = i; 779 } 792 } 780 norm = (baseA[2]-baseA[0]).cross(baseA[1]- 793 norm = (baseA[2]-baseA[0]).cross(baseA[1]-pa); 781 if (norm.mag2() > kCarTolerance) 794 if (norm.mag2() > kCarTolerance) 782 { 795 { 783 pPlanes.emplace_back(norm,pa); << 796 pPlanes.push_back(G4Plane3D(norm,pa)); 784 } 797 } 785 } 798 } 786 else if (na == 1) // baseA has one vertex 799 else if (na == 1) // baseA has one vertex 787 { 800 { 788 std::size_t k = nb - 1; << 801 G4int k = nb - 1; 789 for (std::size_t i=0; i<nb; ++i) << 802 for (G4int i=0; i<nb; ++i) 790 { 803 { 791 norm = (baseB[i]-baseA[0]).cross(baseB[k 804 norm = (baseB[i]-baseA[0]).cross(baseB[k]-baseA[0]); 792 if (norm.mag2() > kCarTolerance) 805 if (norm.mag2() > kCarTolerance) 793 { 806 { 794 pPlanes.emplace_back(norm,baseA[0]); << 807 pPlanes.push_back(G4Plane3D(norm,baseA[0])); 795 } 808 } 796 k = i; 809 k = i; 797 } 810 } 798 norm = (baseB[2]-baseB[0]).cross(baseB[1]- 811 norm = (baseB[2]-baseB[0]).cross(baseB[1]-pb); 799 if (norm.mag2() > kCarTolerance) 812 if (norm.mag2() > kCarTolerance) 800 { 813 { 801 pPlanes.emplace_back(norm,pb); << 814 pPlanes.push_back(G4Plane3D(norm,pb)); 802 } 815 } 803 } 816 } 804 817 805 // Ensure that normals of the planes point t 818 // Ensure that normals of the planes point to outside 806 // 819 // 807 std::size_t nplanes = pPlanes.size(); << 820 G4int nplanes = pPlanes.size(); 808 for (std::size_t i=0; i<nplanes; ++i) << 821 for (G4int i=0; i<nplanes; ++i) 809 { 822 { 810 pPlanes[i].normalize(); 823 pPlanes[i].normalize(); 811 if (pPlanes[i].distance(p0) > 0) 824 if (pPlanes[i].distance(p0) > 0) 812 { 825 { 813 pPlanes[i] = G4Plane3D(-pPlanes[i].a(),- 826 pPlanes[i] = G4Plane3D(-pPlanes[i].a(),-pPlanes[i].b(), 814 -pPlanes[i].c(),- 827 -pPlanes[i].c(),-pPlanes[i].d()); 815 } << 828 } 816 } 829 } 817 } 830 } 818 831 819 ////////////////////////////////////////////// 832 /////////////////////////////////////////////////////////////////////// 820 // 833 // 821 // Clip edges of a prism by G4VoxelLimits box. 834 // Clip edges of a prism by G4VoxelLimits box. Return true if all edges 822 // are inside or intersect the voxel, in this 835 // are inside or intersect the voxel, in this case further calculations 823 // are not needed << 836 // are not needed 824 // 837 // 825 G4bool 838 G4bool 826 G4BoundingEnvelope::ClipEdgesByVoxel(const std 839 G4BoundingEnvelope::ClipEdgesByVoxel(const std::vector<G4Segment3D>& pEdges, 827 const G4V 840 const G4VoxelLimits& pBox, 828 G4S 841 G4Segment3D& pExtent) const 829 { 842 { 830 G4bool done = true; 843 G4bool done = true; 831 G4Point3D emin = pExtent.first; 844 G4Point3D emin = pExtent.first; 832 G4Point3D emax = pExtent.second; 845 G4Point3D emax = pExtent.second; 833 846 834 std::size_t nedges = pEdges.size(); << 847 G4int nedges = pEdges.size(); 835 for (std::size_t k=0; k<nedges; ++k) << 848 for (G4int k=0; k<nedges; ++k) 836 { 849 { 837 G4Point3D p1 = pEdges[k].first; 850 G4Point3D p1 = pEdges[k].first; 838 G4Point3D p2 = pEdges[k].second; 851 G4Point3D p2 = pEdges[k].second; 839 if (std::abs(p1.x()-p2.x())+ 852 if (std::abs(p1.x()-p2.x())+ 840 std::abs(p1.y()-p2.y())+ 853 std::abs(p1.y()-p2.y())+ 841 std::abs(p1.z()-p2.z()) < kCarToleranc 854 std::abs(p1.z()-p2.z()) < kCarTolerance) continue; 842 G4double d1, d2; 855 G4double d1, d2; 843 // Clip current edge by X min 856 // Clip current edge by X min 844 d1 = pBox.GetMinXExtent() - p1.x(); << 857 d1 = pBox.GetMinXExtent() - p1.x(); 845 d2 = pBox.GetMinXExtent() - p2.x(); << 858 d2 = pBox.GetMinXExtent() - p2.x(); 846 if (d1 > 0.0) 859 if (d1 > 0.0) 847 { 860 { 848 if (d2 > 0.0) { done = false; continue; 861 if (d2 > 0.0) { done = false; continue; } // go to next edge 849 p1 = (p2*d1-p1*d2)/(d1-d2); 862 p1 = (p2*d1-p1*d2)/(d1-d2); // move p1 850 } 863 } 851 else 864 else 852 { 865 { 853 if (d2 > 0.0) { p2 = (p1*d2-p2*d1)/(d2-d 866 if (d2 > 0.0) { p2 = (p1*d2-p2*d1)/(d2-d1); } // move p2 854 } 867 } 855 868 856 // Clip current edge by X max 869 // Clip current edge by X max 857 d1 = p1.x() - pBox.GetMaxXExtent(); << 870 d1 = p1.x() - pBox.GetMaxXExtent(); 858 d2 = p2.x() - pBox.GetMaxXExtent(); << 871 d2 = p2.x() - pBox.GetMaxXExtent(); 859 if (d1 > 0.) 872 if (d1 > 0.) 860 { 873 { 861 if (d2 > 0.) { done = false; continue; } 874 if (d2 > 0.) { done = false; continue; } // go to next edge 862 p1 = (p2*d1-p1*d2)/(d1-d2); 875 p1 = (p2*d1-p1*d2)/(d1-d2); 863 } 876 } 864 else 877 else 865 { 878 { 866 if (d2 > 0.) { p2 = (p1*d2-p2*d1)/(d2-d1 879 if (d2 > 0.) { p2 = (p1*d2-p2*d1)/(d2-d1); } 867 } 880 } 868 881 869 // Clip current edge by Y min 882 // Clip current edge by Y min 870 d1 = pBox.GetMinYExtent() - p1.y(); << 883 d1 = pBox.GetMinYExtent() - p1.y(); 871 d2 = pBox.GetMinYExtent() - p2.y(); << 884 d2 = pBox.GetMinYExtent() - p2.y(); 872 if (d1 > 0.) 885 if (d1 > 0.) 873 { 886 { 874 if (d2 > 0.) { done = false; continue; } 887 if (d2 > 0.) { done = false; continue; } // go to next edge 875 p1 = (p2*d1-p1*d2)/(d1-d2); 888 p1 = (p2*d1-p1*d2)/(d1-d2); 876 } 889 } 877 else 890 else 878 { 891 { 879 if (d2 > 0.) { p2 = (p1*d2-p2*d1)/(d2-d1 892 if (d2 > 0.) { p2 = (p1*d2-p2*d1)/(d2-d1); } 880 } 893 } 881 894 882 // Clip current edge by Y max 895 // Clip current edge by Y max 883 d1 = p1.y() - pBox.GetMaxYExtent(); << 896 d1 = p1.y() - pBox.GetMaxYExtent(); 884 d2 = p2.y() - pBox.GetMaxYExtent(); << 897 d2 = p2.y() - pBox.GetMaxYExtent(); 885 if (d1 > 0.) 898 if (d1 > 0.) 886 { 899 { 887 if (d2 > 0.) { done = false; continue; } 900 if (d2 > 0.) { done = false; continue; } // go to next edge 888 p1 = (p2*d1-p1*d2)/(d1-d2); 901 p1 = (p2*d1-p1*d2)/(d1-d2); 889 } 902 } 890 else 903 else 891 { 904 { 892 if (d2 > 0.) { p2 = (p1*d2-p2*d1)/(d2-d1 905 if (d2 > 0.) { p2 = (p1*d2-p2*d1)/(d2-d1); } 893 } 906 } 894 907 895 // Clip current edge by Z min 908 // Clip current edge by Z min 896 d1 = pBox.GetMinZExtent() - p1.z(); << 909 d1 = pBox.GetMinZExtent() - p1.z(); 897 d2 = pBox.GetMinZExtent() - p2.z(); << 910 d2 = pBox.GetMinZExtent() - p2.z(); 898 if (d1 > 0.) 911 if (d1 > 0.) 899 { 912 { 900 if (d2 > 0.) { done = false; continue; } 913 if (d2 > 0.) { done = false; continue; } // go to next edge 901 p1 = (p2*d1-p1*d2)/(d1-d2); 914 p1 = (p2*d1-p1*d2)/(d1-d2); 902 } 915 } 903 else 916 else 904 { 917 { 905 if (d2 > 0.) { p2 = (p1*d2-p2*d1)/(d2-d1 918 if (d2 > 0.) { p2 = (p1*d2-p2*d1)/(d2-d1); } 906 } 919 } 907 920 908 // Clip current edge by Z max 921 // Clip current edge by Z max 909 d1 = p1.z() - pBox.GetMaxZExtent(); << 922 d1 = p1.z() - pBox.GetMaxZExtent(); 910 d2 = p2.z() - pBox.GetMaxZExtent(); << 923 d2 = p2.z() - pBox.GetMaxZExtent(); 911 if (d1 > 0.) 924 if (d1 > 0.) 912 { 925 { 913 if (d2 > 0.) { done = false; continue; } 926 if (d2 > 0.) { done = false; continue; } // go to next edge 914 p1 = (p2*d1-p1*d2)/(d1-d2); 927 p1 = (p2*d1-p1*d2)/(d1-d2); 915 } 928 } 916 else 929 else 917 { 930 { 918 if (d2 > 0.) { p2 = (p1*d2-p2*d1)/(d2-d1 931 if (d2 > 0.) { p2 = (p1*d2-p2*d1)/(d2-d1); } 919 } 932 } 920 933 921 // Adjust current extent 934 // Adjust current extent 922 emin.setX(std::min(std::min(p1.x(),p2.x()) << 935 emin.setX(std::min(std::min(p1.x(),p2.x()),emin.x())); 923 emin.setY(std::min(std::min(p1.y(),p2.y()) << 936 emin.setY(std::min(std::min(p1.y(),p2.y()),emin.y())); 924 emin.setZ(std::min(std::min(p1.z(),p2.z()) << 937 emin.setZ(std::min(std::min(p1.z(),p2.z()),emin.z())); 925 << 938 926 emax.setX(std::max(std::max(p1.x(),p2.x()) << 939 emax.setX(std::max(std::max(p1.x(),p2.x()),emax.x())); 927 emax.setY(std::max(std::max(p1.y(),p2.y()) << 940 emax.setY(std::max(std::max(p1.y(),p2.y()),emax.y())); 928 emax.setZ(std::max(std::max(p1.z(),p2.z()) << 941 emax.setZ(std::max(std::max(p1.z(),p2.z()),emax.z())); 929 } 942 } 930 943 931 // Return true if all edges (at least partia 944 // Return true if all edges (at least partially) are inside 932 // the voxel limits, otherwise return false 945 // the voxel limits, otherwise return false 933 pExtent.first = emin; 946 pExtent.first = emin; 934 pExtent.second = emax; 947 pExtent.second = emax; 935 948 936 return done; 949 return done; 937 } 950 } 938 951 939 ////////////////////////////////////////////// 952 /////////////////////////////////////////////////////////////////////// 940 // 953 // 941 // Clip G4VoxelLimits by set of planes boundin 954 // Clip G4VoxelLimits by set of planes bounding a prism 942 // 955 // 943 void 956 void 944 G4BoundingEnvelope::ClipVoxelByPlanes(G4int pB 957 G4BoundingEnvelope::ClipVoxelByPlanes(G4int pBits, 945 const G4 958 const G4VoxelLimits& pBox, 946 const st 959 const std::vector<G4Plane3D>& pPlanes, 947 const G4 960 const G4Segment3D& pAABB, 948 G4 961 G4Segment3D& pExtent) const 949 { 962 { 950 G4Point3D emin = pExtent.first; 963 G4Point3D emin = pExtent.first; 951 G4Point3D emax = pExtent.second; 964 G4Point3D emax = pExtent.second; 952 965 953 // Create edges of the voxel limits box redu 966 // Create edges of the voxel limits box reducing them where 954 // appropriate to avoid calculations with bi 967 // appropriate to avoid calculations with big numbers (kInfinity) 955 // 968 // 956 G4double xmin = std::max(pBox.GetMinXExtent( 969 G4double xmin = std::max(pBox.GetMinXExtent(),pAABB.first.x() -1.); 957 G4double xmax = std::min(pBox.GetMaxXExtent( 970 G4double xmax = std::min(pBox.GetMaxXExtent(),pAABB.second.x()+1.); 958 971 959 G4double ymin = std::max(pBox.GetMinYExtent( 972 G4double ymin = std::max(pBox.GetMinYExtent(),pAABB.first.y() -1.); 960 G4double ymax = std::min(pBox.GetMaxYExtent( 973 G4double ymax = std::min(pBox.GetMaxYExtent(),pAABB.second.y()+1.); 961 974 962 G4double zmin = std::max(pBox.GetMinZExtent( 975 G4double zmin = std::max(pBox.GetMinZExtent(),pAABB.first.z() -1.); 963 G4double zmax = std::min(pBox.GetMaxZExtent( 976 G4double zmax = std::min(pBox.GetMaxZExtent(),pAABB.second.z()+1.); 964 977 965 std::vector<G4Segment3D> edges(12); 978 std::vector<G4Segment3D> edges(12); 966 G4int i = 0, bits = pBits; 979 G4int i = 0, bits = pBits; 967 if ((bits & 0x001) == 0) << 980 if (!(bits & 0x001)) 968 { << 981 { 969 edges[i ].first.set( xmin,ymin,zmin); << 982 edges[i ].first.set( xmin,ymin,zmin); 970 edges[i++].second.set(xmax,ymin,zmin); 983 edges[i++].second.set(xmax,ymin,zmin); 971 } 984 } 972 if ((bits & 0x002) == 0) << 985 if (!(bits & 0x002)) 973 { 986 { 974 edges[i ].first.set( xmax,ymin,zmin); 987 edges[i ].first.set( xmax,ymin,zmin); 975 edges[i++].second.set(xmax,ymax,zmin); 988 edges[i++].second.set(xmax,ymax,zmin); 976 } 989 } 977 if ((bits & 0x004) == 0) << 990 if (!(bits & 0x004)) 978 { 991 { 979 edges[i ].first.set( xmax,ymax,zmin); 992 edges[i ].first.set( xmax,ymax,zmin); 980 edges[i++].second.set(xmin,ymax,zmin); 993 edges[i++].second.set(xmin,ymax,zmin); 981 } 994 } 982 if ((bits & 0x008) == 0) << 995 if (!(bits & 0x008)) 983 { 996 { 984 edges[i ].first.set( xmin,ymax,zmin); 997 edges[i ].first.set( xmin,ymax,zmin); 985 edges[i++].second.set(xmin,ymin,zmin); << 998 edges[i++].second.set(xmin,ymin,zmin); 986 } 999 } 987 1000 988 if ((bits & 0x010) == 0) << 1001 if (!(bits & 0x010)) 989 { << 1002 { 990 edges[i ].first.set( xmin,ymin,zmax); << 1003 edges[i ].first.set( xmin,ymin,zmax); 991 edges[i++].second.set(xmax,ymin,zmax); 1004 edges[i++].second.set(xmax,ymin,zmax); 992 } 1005 } 993 if ((bits & 0x020) == 0) << 1006 if (!(bits & 0x020)) 994 { 1007 { 995 edges[i ].first.set( xmax,ymin,zmax); 1008 edges[i ].first.set( xmax,ymin,zmax); 996 edges[i++].second.set(xmax,ymax,zmax); 1009 edges[i++].second.set(xmax,ymax,zmax); 997 } 1010 } 998 if ((bits & 0x040) == 0) << 1011 if (!(bits & 0x040)) 999 { 1012 { 1000 edges[i ].first.set( xmax,ymax,zmax); 1013 edges[i ].first.set( xmax,ymax,zmax); 1001 edges[i++].second.set(xmin,ymax,zmax); 1014 edges[i++].second.set(xmin,ymax,zmax); 1002 } 1015 } 1003 if ((bits & 0x080) == 0) << 1016 if (!(bits & 0x080)) 1004 { 1017 { 1005 edges[i ].first.set( xmin,ymax,zmax); 1018 edges[i ].first.set( xmin,ymax,zmax); 1006 edges[i++].second.set(xmin,ymin,zmax); << 1019 edges[i++].second.set(xmin,ymin,zmax); 1007 } 1020 } 1008 1021 1009 if ((bits & 0x100) == 0) << 1022 if (!(bits & 0x100)) 1010 { << 1023 { 1011 edges[i ].first.set( xmin,ymin,zmin); << 1024 edges[i ].first.set( xmin,ymin,zmin); 1012 edges[i++].second.set(xmin,ymin,zmax); << 1025 edges[i++].second.set(xmin,ymin,zmax); 1013 } 1026 } 1014 if ((bits & 0x200) == 0) << 1027 if (!(bits & 0x200)) 1015 { 1028 { 1016 edges[i ].first.set( xmax,ymin,zmin); 1029 edges[i ].first.set( xmax,ymin,zmin); 1017 edges[i++].second.set(xmax,ymin,zmax); 1030 edges[i++].second.set(xmax,ymin,zmax); 1018 } 1031 } 1019 if ((bits & 0x400) == 0) << 1032 if (!(bits & 0x400)) 1020 { 1033 { 1021 edges[i ].first.set( xmax,ymax,zmin); 1034 edges[i ].first.set( xmax,ymax,zmin); 1022 edges[i++].second.set(xmax,ymax,zmax); 1035 edges[i++].second.set(xmax,ymax,zmax); 1023 } 1036 } 1024 if ((bits & 0x800) == 0) << 1037 if (!(bits & 0x800)) 1025 { 1038 { 1026 edges[i ].first.set( xmin,ymax,zmin); 1039 edges[i ].first.set( xmin,ymax,zmin); 1027 edges[i++].second.set(xmin,ymax,zmax); 1040 edges[i++].second.set(xmin,ymax,zmax); 1028 } 1041 } 1029 edges.resize(i); 1042 edges.resize(i); 1030 1043 1031 // Clip the edges by the planes 1044 // Clip the edges by the planes 1032 // 1045 // 1033 for (const auto & edge : edges) << 1046 for (auto iedge = edges.cbegin(); iedge != edges.cend(); ++iedge) 1034 { 1047 { 1035 G4bool exist = true; 1048 G4bool exist = true; 1036 G4Point3D p1 = edge.first; << 1049 G4Point3D p1 = iedge->first; 1037 G4Point3D p2 = edge.second; << 1050 G4Point3D p2 = iedge->second; 1038 for (const auto & plane : pPlanes) << 1051 for (auto iplane = pPlanes.cbegin(); iplane != pPlanes.cend(); ++iplane) 1039 { 1052 { 1040 // Clip current edge 1053 // Clip current edge 1041 G4double d1 = plane.distance(p1); << 1054 G4double d1 = iplane->distance(p1); 1042 G4double d2 = plane.distance(p2); << 1055 G4double d2 = iplane->distance(p2); 1043 if (d1 > 0.0) 1056 if (d1 > 0.0) 1044 { 1057 { 1045 if (d2 > 0.0) { exist = false; break; 1058 if (d2 > 0.0) { exist = false; break; } // go to next edge 1046 p1 = (p2*d1-p1*d2)/(d1-d2); 1059 p1 = (p2*d1-p1*d2)/(d1-d2); // move p1 1047 } 1060 } 1048 else 1061 else 1049 { 1062 { 1050 if (d2 > 0.0) { p2 = (p1*d2-p2*d1)/(d 1063 if (d2 > 0.0) { p2 = (p1*d2-p2*d1)/(d2-d1); } // move p2 1051 } 1064 } 1052 } 1065 } 1053 // Adjust the extent 1066 // Adjust the extent 1054 if (exist) 1067 if (exist) 1055 { 1068 { 1056 emin.setX(std::min(std::min(p1.x(),p2.x << 1069 emin.setX(std::min(std::min(p1.x(),p2.x()),emin.x())); 1057 emin.setY(std::min(std::min(p1.y(),p2.y << 1070 emin.setY(std::min(std::min(p1.y(),p2.y()),emin.y())); 1058 emin.setZ(std::min(std::min(p1.z(),p2.z << 1071 emin.setZ(std::min(std::min(p1.z(),p2.z()),emin.z())); 1059 << 1072 1060 emax.setX(std::max(std::max(p1.x(),p2.x << 1073 emax.setX(std::max(std::max(p1.x(),p2.x()),emax.x())); 1061 emax.setY(std::max(std::max(p1.y(),p2.y << 1074 emax.setY(std::max(std::max(p1.y(),p2.y()),emax.y())); 1062 emax.setZ(std::max(std::max(p1.z(),p2.z << 1075 emax.setZ(std::max(std::max(p1.z(),p2.z()),emax.z())); 1063 } 1076 } 1064 } 1077 } 1065 1078 1066 // Copy the extent back 1079 // Copy the extent back 1067 // 1080 // 1068 pExtent.first = emin; 1081 pExtent.first = emin; 1069 pExtent.second = emax; 1082 pExtent.second = emax; 1070 } 1083 } 1071 1084