Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/management/src/G4BoundingEnvelope.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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