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