| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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()->GetSurfa
40
41 //////////////////////////////////////////////
42 //
43 // Constructor from an axis aligned bounding b
44 //
45 G4BoundingEnvelope::G4BoundingEnvelope(const G
46 const G
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 G4T
60 : fPolygons(&polygons)
61 {
62 // Check correctness of polygons
63 //
64 CheckBoundingPolygons();
65
66 // Set bounding box
67 //
68 G4double xmin = kInfinity, ymin = kInfinit
69 G4double xmax = -kInfinity, ymax = -kInfinit
70 for (const auto & polygon : *fPolygons)
71 {
72 for (auto ipoint = polygon->cbegin(); ipoi
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 seque
96 //
97 G4BoundingEnvelope::
98 G4BoundingEnvelope( const G4ThreeVector& pMin,
99 const G4ThreeVector& pMax,
100 const std::vector<const G4
101 : fMin(pMin), fMax(pMax), fPolygons(&polygon
102 {
103 // Check correctness of bounding box and pol
104 //
105 CheckBoundingBox();
106 CheckBoundingPolygons();
107 }
108
109 //////////////////////////////////////////////
110 //
111 // Check correctness of the axis aligned bound
112 //
113 void G4BoundingEnvelope::CheckBoundingBox()
114 {
115 if (fMin.x() >= fMax.x() || fMin.y() >= fMax
116 {
117 std::ostringstream message;
118 message << "Badly defined bounding box (mi
119 << "\npMin = " << fMin
120 << "\npMax = " << fMax;
121 G4Exception("G4BoundingEnvelope::CheckBoun
122 "GeomMgt0001", JustWarning, me
123 }
124 }
125
126 //////////////////////////////////////////////
127 //
128 // Check correctness of the sequence of boundi
129 // Firsf and last polygons may consist of a si
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 th
138 << "\nShould be at least two!";
139 G4Exception("G4BoundingEnvelope::CheckBoun
140 "GeomMgt0001", FatalException,
141 return;
142 }
143
144 std::size_t nsize = std::max((*fPolygons)[0
145 if (nsize < 3)
146 {
147 std::ostringstream message;
148 message << "Badly constructed polygons!"
149 << "\nNumber of polygons: " << nba
150 << "\nPolygon #0 size: " << (*fPol
151 << "\nPolygon #1 size: " << (*fPol
152 << "\n...";
153 G4Exception("G4BoundingEnvelope::CheckBoun
154 "GeomMgt0001", FatalException,
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: " << nba
167 << "\nPolygon #" << k << " size: "
168 << "\nexpected size: " << nsize;
169 G4Exception("G4BoundingEnvelope::SetBoundi
170 "GeomMgt0001", FatalException,
171 return;
172 }
173 }
174
175 //////////////////////////////////////////////
176 //
177 // Quick comparison: bounding box vs voxel, it
178 // calculations are not needed
179 //
180 G4bool
181 G4BoundingEnvelope::
182 BoundingBoxVsVoxelLimits(const EAxis pAxis,
183 const G4VoxelLimits&
184 const G4Transform3D&
185 G4double& pMin, G4dou
186 {
187 pMin = kInfinity;
188 pMax = -kInfinity;
189 G4double xminlim = pVoxelLimits.GetMinXExten
190 G4double xmaxlim = pVoxelLimits.GetMaxXExten
191 G4double yminlim = pVoxelLimits.GetMinYExten
192 G4double ymaxlim = pVoxelLimits.GetMaxYExten
193 G4double zminlim = pVoxelLimits.GetMinZExten
194 G4double zmaxlim = pVoxelLimits.GetMaxZExten
195
196 // Special case of pure translation
197 //
198 if (pTransform3D.xx()==1 && pTransform3D.yy(
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 t
208 if (xmax+kCarTolerance < xminlim) return t
209 if (ymin-kCarTolerance > ymaxlim) return t
210 if (ymax+kCarTolerance < yminlim) return t
211 if (zmin-kCarTolerance > zmaxlim) return t
212 if (zmax+kCarTolerance < zminlim) return t
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)
221 pMax = (xmax+kCarTolerance > xmaxlim)
222 }
223 else if (pAxis == kYAxis)
224 {
225 pMin = (ymin-kCarTolerance < yminlim)
226 pMax = (ymax+kCarTolerance > ymaxlim)
227 }
228 else if (pAxis == kZAxis)
229 {
230 pMin = (zmin-kCarTolerance < zminlim)
231 pMax = (zmax+kCarTolerance > zmaxlim)
232 }
233 pMin -= kCarTolerance;
234 pMax += kCarTolerance;
235 return true;
236 }
237 }
238
239 // Find max scale factor of the transformati
240 // equal to kCarTolerance multiplied by the
241 //
242 G4double scale = FindScaleFactor(pTransform3
243 G4double delta = kCarTolerance*scale;
244
245 // Set the sphere surrounding the bounding b
246 //
247 G4Point3D center = pTransform3D*G4Point3D(0.
248 G4double radius = 0.5*scale*(fMax-fMin).mag
249
250 // Check if the sphere surrounding the bound
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
265 //
266 G4bool
267 G4BoundingEnvelope::CalculateExtent(const EAxi
268 const G4Vo
269 const G4Tr
270 G4double&
271 {
272 pMin = kInfinity;
273 pMax = -kInfinity;
274 G4double xminlim = pVoxelLimits.GetMinXExten
275 G4double xmaxlim = pVoxelLimits.GetMaxXExten
276 G4double yminlim = pVoxelLimits.GetMinYExten
277 G4double ymaxlim = pVoxelLimits.GetMaxYExten
278 G4double zminlim = pVoxelLimits.GetMinZExten
279 G4double zmaxlim = pVoxelLimits.GetMaxZExten
280
281 // Special case of pure translation
282 //
283 if (pTransform3D.xx()==1 && pTransform3D.yy(
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 f
293 if (xmax+kCarTolerance < xminlim) return f
294 if (ymin-kCarTolerance > ymaxlim) return f
295 if (ymax+kCarTolerance < yminlim) return f
296 if (zmin-kCarTolerance > zmaxlim) return f
297 if (zmax+kCarTolerance < zminlim) return f
298
299 if (fPolygons == nullptr)
300 {
301 if (pAxis == kXAxis)
302 {
303 pMin = (xmin-kCarTolerance < xminlim)
304 pMax = (xmax+kCarTolerance > xmaxlim)
305 }
306 else if (pAxis == kYAxis)
307 {
308 pMin = (ymin-kCarTolerance < yminlim)
309 pMax = (ymax+kCarTolerance > ymaxlim)
310 }
311 else if (pAxis == kZAxis)
312 {
313 pMin = (zmin-kCarTolerance < zminlim)
314 pMax = (zmax+kCarTolerance > zmaxlim)
315 }
316 pMin -= kCarTolerance;
317 pMax += kCarTolerance;
318 return true;
319 }
320 }
321
322 // Find max scale factor of the transformati
323 // equal to kCarTolerance multiplied by the
324 //
325 G4double scale = FindScaleFactor(pTransform3
326 G4double delta = kCarTolerance*scale;
327
328 // Set the sphere surrounding the bounding b
329 //
330 G4Point3D center = pTransform3D*G4Point3D(0.
331 G4double radius = 0.5*scale*(fMax-fMin).mag
332
333 // Check if the sphere surrounding the bound
334 // the voxel limits, if so then transform on
335 //
336 if (center.x()-radius >= xminlim && center.x
337 center.y()-radius >= yminlim && center.y
338 center.z()-radius >= zminlim && center.z
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 = -kInfini
367 if (fPolygons == nullptr)
368 {
369 G4double coor;
370 coor = cx*fMin.x() + cy*fMin.y() + cz*fM
371 if (coor < emin) emin = coor;
372 if (coor > emax) emax = coor;
373 coor = cx*fMax.x() + cy*fMin.y() + cz*fM
374 if (coor < emin) emin = coor;
375 if (coor > emax) emax = coor;
376 coor = cx*fMax.x() + cy*fMax.y() + cz*fM
377 if (coor < emin) emin = coor;
378 if (coor > emax) emax = coor;
379 coor = cx*fMin.x() + cy*fMax.y() + cz*fM
380 if (coor < emin) emin = coor;
381 if (coor > emax) emax = coor;
382 coor = cx*fMin.x() + cy*fMin.y() + cz*fM
383 if (coor < emin) emin = coor;
384 if (coor > emax) emax = coor;
385 coor = cx*fMax.x() + cy*fMin.y() + cz*fM
386 if (coor < emin) emin = coor;
387 if (coor > emax) emax = coor;
388 coor = cx*fMax.x() + cy*fMax.y() + cz*fM
389 if (coor < emin) emin = coor;
390 if (coor > emax) emax = coor;
391 coor = cx*fMin.x() + cy*fMax.y() + cz*fM
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(); ip
400 {
401 G4double coor = ipoint->x()*cx + ipo
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 bound
413 // the voxel limits
414 //
415 if (center.x()-radius > xmaxlim) return fals
416 if (center.y()-radius > ymaxlim) return fals
417 if (center.z()-radius > zmaxlim) return fals
418 if (center.x()+radius < xminlim) return fals
419 if (center.y()+radius < yminlim) return fals
420 if (center.z()+radius < zminlim) return fals
421
422 // Transform polygons
423 //
424 std::vector<G4Point3D> vertices;
425 std::vector<std::pair<G4int, G4int>> bases;
426 TransformVertices(pTransform3D, vertices, ba
427 std::size_t nbases = bases.size();
428
429 // Create adjusted G4VoxelLimits box. New li
430 // delta, kCarTolerance multiplied by max sc
431 // the transformation
432 //
433 EAxis axes[] = { kXAxis, kYAxis, kZAxis };
434 G4VoxelLimits limits; // default is unlimite
435 for (const auto & iAxis : axes)
436 {
437 if (pVoxelLimits.IsLimited(iAxis))
438 {
439 G4double emin = pVoxelLimits.GetMinExten
440 G4double emax = pVoxelLimits.GetMaxExten
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, kInfin
450 extent.second = G4Point3D(-kInfinity,-kInfin
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; +
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
466 if (prismAABB.first.x() >= limits.GetMinXE
467 prismAABB.first.y() >= limits.GetMinYE
468 prismAABB.first.z() >= limits.GetMinZE
469 prismAABB.second.x()<= limits.GetMaxXE
470 prismAABB.second.y()<= limits.GetMaxYE
471 prismAABB.second.z()<= limits.GetMaxZE
472 {
473 if (extent.first.x() > prismAABB.first.
474 extent.first.setX( prismAABB.first.x()
475 if (extent.first.y() > prismAABB.first.
476 extent.first.setY( prismAABB.first.y()
477 if (extent.first.z() > prismAABB.first.
478 extent.first.setZ( prismAABB.first.z()
479 if (extent.second.x() < prismAABB.second
480 extent.second.setX(prismAABB.second.x(
481 if (extent.second.y() < prismAABB.second
482 extent.second.setY(prismAABB.second.y(
483 if (extent.second.z() < prismAABB.second
484 extent.second.setZ(prismAABB.second.z(
485 continue;
486 }
487
488 // Check if prismAABB is outside the voxel
489 if (prismAABB.first.x() > limits.GetMaxXE
490 if (prismAABB.first.y() > limits.GetMaxYE
491 if (prismAABB.first.z() > limits.GetMaxZE
492 if (prismAABB.second.x() < limits.GetMinXE
493 if (prismAABB.second.y() < limits.GetMinYE
494 if (prismAABB.second.z() < limits.GetMinZE
495
496 // Clip edges of the prism by adjusted G4V
497 std::vector<G4Segment3D> vecEdges;
498 CreateListOfEdges(baseA, baseB, vecEdges);
499 if (ClipEdgesByVoxel(vecEdges, limits, ext
500
501 // Some edges of the prism are completely
502 // limits, clip selected edges (see bits)
503 // by the prism
504 G4int bits = 0x000;
505 if (limits.GetMinXExtent() < prismAABB.fir
506 bits |= 0x988; // 1001 1000 1000
507 if (limits.GetMaxXExtent() > prismAABB.sec
508 bits |= 0x622; // 0110 0010 0010
509
510 if (limits.GetMinYExtent() < prismAABB.fir
511 bits |= 0x311; // 0011 0001 0001
512 if (limits.GetMaxYExtent() > prismAABB.sec
513 bits |= 0xC44; // 1100 0100 0100
514
515 if (limits.GetMinZExtent() < prismAABB.fir
516 bits |= 0x00F; // 0000 0000 1111
517 if (limits.GetMaxZExtent() > prismAABB.sec
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,
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
530 if (pAxis == kYAxis) { emin = extent.first.y
531 if (pAxis == kZAxis) { emin = extent.first.z
532
533 if (emin > emax) return false;
534 emin -= delta;
535 emax += delta;
536 G4double minlim = pVoxelLimits.GetMinExtent(
537 G4double maxlim = pVoxelLimits.GetMaxExtent(
538 pMin = (emin < minlim) ? minlim-kCarToleranc
539 pMax = (emax > maxlim) ? maxlim+kCarToleranc
540 return true;
541 }
542
543 //////////////////////////////////////////////
544 //
545 // Find max scale factor of the transformation
546 //
547 G4double
548 G4BoundingEnvelope::FindScaleFactor(const G4Tr
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),s
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& pTransf
577 std::vector<G4Point3D>& pVer
578 std::vector<std::pair<G4int,
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.c
596 auto iaend = (fPolygons == nullptr) ? aabb.c
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.getTra
615 for (auto i = ia; i != iaend; ++i)
616 for (auto k = (*i)->cbegin(); k != (*i)-
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)-
623 pVertices.push_back(pTransform3D*G4Poi
624 }
625 }
626
627 //////////////////////////////////////////////
628 //
629 // Find bounding box of a prism
630 //
631 void
632 G4BoundingEnvelope::GetPrismAABB(const G4Polyg
633 const G4Polyg
634 G4Segme
635 {
636 G4double xmin = kInfinity, ymin = kInfinit
637 G4double xmax = -kInfinity, ymax = -kInfinit
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 G4
681 const G4
682 std::vec
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 G
729 const G
730 std::ve
731 {
732 // Find centers of the bases and internal po
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
739 for (std::size_t i=0; i<nb; ++i) pb += baseB
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
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
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]-
758 if (norm.mag2() > kCarTolerance)
759 {
760 pPlanes.emplace_back(norm,pa);
761 }
762 norm = (baseB[2]-baseB[0]).cross(baseB[1]-
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
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]-
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
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]-
799 if (norm.mag2() > kCarTolerance)
800 {
801 pPlanes.emplace_back(norm,pb);
802 }
803 }
804
805 // Ensure that normals of the planes point t
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(),-
814 -pPlanes[i].c(),-
815 }
816 }
817 }
818
819 //////////////////////////////////////////////
820 //
821 // Clip edges of a prism by G4VoxelLimits box.
822 // are inside or intersect the voxel, in this
823 // are not needed
824 //
825 G4bool
826 G4BoundingEnvelope::ClipEdgesByVoxel(const std
827 const G4V
828 G4S
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()) < kCarToleranc
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;
849 p1 = (p2*d1-p1*d2)/(d1-d2);
850 }
851 else
852 {
853 if (d2 > 0.0) { p2 = (p1*d2-p2*d1)/(d2-d
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; }
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; }
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; }
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; }
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; }
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())
923 emin.setY(std::min(std::min(p1.y(),p2.y())
924 emin.setZ(std::min(std::min(p1.z(),p2.z())
925
926 emax.setX(std::max(std::max(p1.x(),p2.x())
927 emax.setY(std::max(std::max(p1.y(),p2.y())
928 emax.setZ(std::max(std::max(p1.z(),p2.z())
929 }
930
931 // Return true if all edges (at least partia
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 boundin
942 //
943 void
944 G4BoundingEnvelope::ClipVoxelByPlanes(G4int pB
945 const G4
946 const st
947 const G4
948 G4
949 {
950 G4Point3D emin = pExtent.first;
951 G4Point3D emax = pExtent.second;
952
953 // Create edges of the voxel limits box redu
954 // appropriate to avoid calculations with bi
955 //
956 G4double xmin = std::max(pBox.GetMinXExtent(
957 G4double xmax = std::min(pBox.GetMaxXExtent(
958
959 G4double ymin = std::max(pBox.GetMinYExtent(
960 G4double ymax = std::min(pBox.GetMaxYExtent(
961
962 G4double zmin = std::max(pBox.GetMinZExtent(
963 G4double zmax = std::min(pBox.GetMaxZExtent(
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;
1046 p1 = (p2*d1-p1*d2)/(d1-d2);
1047 }
1048 else
1049 {
1050 if (d2 > 0.0) { p2 = (p1*d2-p2*d1)/(d
1051 }
1052 }
1053 // Adjust the extent
1054 if (exist)
1055 {
1056 emin.setX(std::min(std::min(p1.x(),p2.x
1057 emin.setY(std::min(std::min(p1.y(),p2.y
1058 emin.setZ(std::min(std::min(p1.z(),p2.z
1059
1060 emax.setX(std::max(std::max(p1.x(),p2.x
1061 emax.setY(std::max(std::max(p1.y(),p2.y
1062 emax.setZ(std::max(std::max(p1.z(),p2.z
1063 }
1064 }
1065
1066 // Copy the extent back
1067 //
1068 pExtent.first = emin;
1069 pExtent.second = emax;
1070 }
1071