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