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