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