| 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 // G4Voxelizer implementation
27 //
28 // 19.10.12 Marek Gayer, created
29 // -------------------------------------------
30
31 #include <iostream>
32 #include <iomanip>
33 #include <sstream>
34 #include <algorithm>
35 #include <set>
36
37 #include "G4VSolid.hh"
38
39 #include "G4CSGSolid.hh"
40 #include "G4GeometryTolerance.hh"
41 #include "G4Orb.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SolidStore.hh"
44 #include "G4Types.hh"
45 #include "G4Voxelizer.hh"
46 #include "Randomize.hh"
47 #include "geomdefs.hh"
48
49 using namespace std;
50
51 G4ThreadLocal G4int G4Voxelizer::fDefaultVoxel
52
53 //____________________________________________
54 G4Voxelizer::G4Voxelizer()
55 : fBoundingBox("VoxBBox", 1, 1, 1)
56 {
57 fCountOfVoxels = fNPerSlice = fTotalCandidat
58
59 fTolerance = G4GeometryTolerance::GetInstanc
60
61 SetMaxVoxels(fDefaultVoxelsCount);
62
63 G4SolidStore::GetInstance()->DeRegister(&fBo
64 }
65
66 //____________________________________________
67 G4Voxelizer::~G4Voxelizer() = default;
68
69 //____________________________________________
70 void G4Voxelizer::BuildEmpty()
71 {
72 // by reserving the size of candidates, we w
73 // the vector which could cause fragmentatio
74 //
75 std::vector<G4int> xyz(3), max(3), candidate
76 const std::vector<G4int> empty(0);
77
78 for (auto i = 0; i <= 2; ++i) max[i] = (G4in
79 unsigned int size = max[0] * max[1] * max[2]
80
81 fEmpty.Clear();
82 fEmpty.ResetBitNumber(size-1);
83 fEmpty.ResetAllBits(true);
84
85 for (xyz[2] = 0; xyz[2] < max[2]; ++xyz[2])
86 {
87 for (xyz[1] = 0; xyz[1] < max[1]; ++xyz[1]
88 {
89 for (xyz[0] = 0; xyz[0] < max[0]; ++xyz[
90 {
91 if (GetCandidatesVoxelArray(xyz, candi
92 {
93 G4int index = GetVoxelsIndex(xyz);
94 fEmpty.SetBitNumber(index, false);
95
96 // rather than assigning directly wi
97 // "fCandidates[index] = candidates;
98 // capacity would be just exact, we
99 //
100 std::vector<G4int> &c = (fCandidates
101 c.reserve(candidates.size());
102 c.assign(candidates.begin(), candida
103 }
104 }
105 }
106 }
107 #ifdef G4SPECSDEBUG
108 G4cout << "Non-empty voxels count: " << fCan
109 #endif
110 }
111
112 //____________________________________________
113 void G4Voxelizer::BuildVoxelLimits(std::vector
114 std::vector
115 {
116 // "BuildVoxelLimits"'s aim is to store the
117 // well as the half lengths related to the b
118 // These quantities are stored in the array
119 // node
120 //
121 if (std::size_t numNodes = solids.size()) //
122 {
123 fBoxes.resize(numNodes); // Array which wi
124 fNPerSlice = G4int(1 + (fBoxes.size() - 1)
125
126 // related to a particular node, but also
127 // the coordinates of its origin
128
129 G4ThreeVector toleranceVector(fTolerance,f
130
131 for (std::size_t i = 0; i < numNodes; ++i)
132 {
133 G4VSolid& solid = *solids[i];
134 G4Transform3D transform = transforms[i];
135 G4ThreeVector min, max;
136
137 solid.BoundingLimits(min, max);
138 if (solid.GetEntityType() == "G4Orb")
139 {
140 G4Orb& orb = *(G4Orb*) &solid;
141 G4ThreeVector orbToleranceVector;
142 G4double tolerance = orb.GetRadialTole
143 orbToleranceVector.set(tolerance,toler
144 min -= orbToleranceVector;
145 max += orbToleranceVector;
146 }
147 else
148 {
149 min -= toleranceVector;
150 max += toleranceVector;
151 }
152 TransformLimits(min, max, transform);
153 fBoxes[i].hlen = (max - min) / 2.;
154 fBoxes[i].pos = (max + min) / 2.;
155 }
156 fTotalCandidates = (G4int)fBoxes.size();
157 }
158 }
159
160 //____________________________________________
161 void G4Voxelizer::BuildVoxelLimits(std::vector
162 {
163 // "BuildVoxelLimits"'s aim is to store the
164 // as the half lengths related to the boundi
165 // These quantities are stored in the array
166 // node.
167
168 if (std::size_t numNodes = facets.size()) //
169 {
170 fBoxes.resize(numNodes); // Array which wi
171 fNPerSlice = G4int(1+(fBoxes.size()-1)/(8*
172
173 G4ThreeVector toleranceVector(10*fToleranc
174
175 for (std::size_t i = 0; i < numNodes; ++i)
176 {
177 G4VFacet &facet = *facets[i];
178 G4ThreeVector min, max;
179 G4ThreeVector x(1,0,0), y(0,1,0), z(0,0,
180 G4ThreeVector extent;
181 max.set (facet.Extent(x), facet.Extent(y
182 min.set (-facet.Extent(-x), -facet.Exten
183 min -= toleranceVector;
184 max += toleranceVector;
185 G4ThreeVector hlen = (max - min) / 2;
186 fBoxes[i].hlen = hlen;
187 fBoxes[i].pos = min + hlen;
188 }
189 fTotalCandidates = (G4int)fBoxes.size();
190 }
191 }
192
193 //____________________________________________
194 void G4Voxelizer::DisplayVoxelLimits() const
195 {
196 // "DisplayVoxelLimits" displays the dX, dY,
197
198 std::size_t numNodes = fBoxes.size();
199 G4long oldprec = G4cout.precision(16);
200 for(std::size_t i = 0; i < numNodes; ++i)
201 {
202 G4cout << setw(10) << setiosflags(ios::fix
203 " -> Node " << i+1 << ":\n" <<
204 "\t * [x,y,z] = " << fBoxes[i].hlen <<
205 "\t * [x,y,z] = " << fBoxes[i].pos << "\
206 }
207 G4cout.precision(oldprec);
208 }
209
210 //____________________________________________
211 void G4Voxelizer::CreateSortedBoundary(std::ve
212 G4int a
213 {
214 // "CreateBoundaries"'s aim is to determine
215 // bounding fBoxes, along each axis. The cre
216 // in the array "boundariesRaw"
217
218 std::size_t numNodes = fBoxes.size(); // Num
219
220 // Determination of the boundaries along x,
221 //
222 for(std::size_t i = 0 ; i < numNodes; ++i)
223 {
224 // For each node, the boundaries are creat
225 // built in method "BuildVoxelLimits"
226 //
227 G4double p = fBoxes[i].pos[axis], d = fBox
228
229 // x boundaries
230 //
231 #ifdef G4SPECSDEBUG
232 G4cout << "Boundary " << p - d << " - " <<
233 #endif
234 boundary[2*i] = p - d;
235 boundary[2*i+1] = p + d;
236 }
237 std::sort(boundary.begin(), boundary.end());
238 }
239
240 //____________________________________________
241 void G4Voxelizer::BuildBoundaries()
242 {
243 // "SortBoundaries" orders the boundaries al
244 // and also does not take into account redun
245 // boundaries are separated by a distance st
246 // The sorted boundaries are respectively st
247 // * boundaries[0..2]
248
249 // In addition, the number of elements conta
250 // are precise thanks to variables: boundari
251 // boundariesCountZ.
252
253 if (std::size_t numNodes = fBoxes.size())
254 {
255 const G4double tolerance = fTolerance / 10
256 // Minimal distance to discriminate two
257
258 std::vector<G4double> sortedBoundary(2*num
259
260 for (auto j = 0; j <= 2; ++j)
261 {
262 CreateSortedBoundary(sortedBoundary, j);
263 std::vector<G4double> &boundary = fBound
264 boundary.clear();
265
266 for(std::size_t i = 0 ; i < 2*numNodes;
267 {
268 G4double newBoundary = sortedBoundary[
269 #ifdef G4SPECSDEBUG
270 if (j == 0) G4cout << "Examining " <<
271 #endif
272 auto size = (G4int)boundary.size();
273 if((size == 0) || std::abs(boundary[si
274 {
275 {
276 #ifdef G4SPECSDEBUG
277 if (j == 0) G4cout << "Adding boun
278 << G4endl;
279 #endif
280 boundary.push_back(newBoundary);
281 continue;
282 }
283 }
284 // If two successive boundaries are to
285 // only the first one is considered
286 }
287
288 auto n = (G4int)boundary.size();
289 G4int max = 100000;
290 if (n > max/2)
291 {
292 G4int skip = n / (max /2); // n has to
293 // therefor
294 std::vector<G4double> reduced;
295 for (G4int i = 0; i < n; ++i)
296 {
297 // 50 ok for 2k, 1000, 2000
298 auto size = (G4int)boundary.size();
299 if (i % skip == 0 || i == 0 || i ==
300 {
301 // this condition of merging bound
302 // it did not count with right par
303 // completely ommited and not incl
304 // Now should be OK
305 //
306 reduced.push_back(boundary[i]);
307 }
308 }
309 boundary = std::move(reduced);
310 }
311 }
312 }
313 }
314
315 //____________________________________________
316 void G4Voxelizer::DisplayBoundaries()
317 {
318 char axis[3] = {'X', 'Y', 'Z'};
319 for (auto i = 0; i <= 2; ++i)
320 {
321 G4cout << " * " << axis[i] << " axis:" <<
322 DisplayBoundaries(fBoundaries[i]);
323 }
324 }
325
326 //____________________________________________
327 void G4Voxelizer::DisplayBoundaries(std::vecto
328 {
329 // Prints the positions of the boundaries of
330
331 std::size_t count = boundaries.size();
332 G4long oldprec = G4cout.precision(16);
333 for(std::size_t i = 0; i < count; ++i)
334 {
335 G4cout << setw(10) << setiosflags(ios::fix
336 if(i != count-1) G4cout << "-> ";
337 }
338 G4cout << "|" << G4endl << "Number of bounda
339 G4cout.precision(oldprec);
340 }
341
342 //____________________________________________
343 void G4Voxelizer::BuildBitmasks(std::vector<G4
344 G4SurfBits bit
345 {
346 // "BuildListNodes" stores in the bitmasks s
347 // along an axis.
348
349 std::size_t numNodes = fBoxes.size();
350 G4int bitsPerSlice = GetBitsPerSlice();
351
352 for (auto k = 0; k < 3; ++k)
353 {
354 std::vector<G4double>& boundary = boundari
355 G4int voxelsCount = (G4int)boundary.size()
356 G4SurfBits& bitmask = bitmasks[k];
357
358 if (!countsOnly)
359 {
360 bitmask.Clear();
361 #ifdef G4SPECSDEBUG
362 G4cout << "Allocating bitmask..." << G4e
363 #endif
364 bitmask.SetBitNumber(voxelsCount*bitsPer
365 // it is here so we can set the maximu
366 // will rellocate the memory and set a
367 }
368 std::vector<G4int>& candidatesCount = fCan
369 candidatesCount.resize(voxelsCount);
370
371 for(G4int i = 0 ; i < voxelsCount; ++i) {
372
373 // Loop on the nodes, number of slices per
374 //
375 for(std::size_t j = 0 ; j < numNodes; ++j)
376 {
377 // Determination of the minimum and maxi
378 // of the bounding boxe of each node
379 //
380 G4double p = fBoxes[j].pos[k], d = fBoxe
381
382 G4double min = p - d; // - localToleranc
383 G4double max = p + d; // + localToleranc
384
385 G4int i = BinarySearch(boundary, min);
386 if (i < 0) { i = 0; }
387
388 do // Loop checking, 13.08.2015, G.Co
389 {
390 if (!countsOnly)
391 {
392 bitmask.SetBitNumber(i*bitsPerSlice+
393 }
394 candidatesCount[i]++;
395 ++i;
396 }
397 while (max > boundary[i] && i < voxelsCo
398 }
399 }
400 #ifdef G4SPECSDEBUG
401 G4cout << "Build list nodes completed." << G
402 #endif
403 }
404
405 //____________________________________________
406 G4String G4Voxelizer::GetCandidatesAsString(co
407 {
408 // Decodes the candidates in mask as G4Strin
409
410 stringstream ss;
411 auto numNodes = (G4int)fBoxes.size();
412
413 for(auto i=0; i<numNodes; ++i)
414 {
415 if (bits.TestBitNumber(i)) { ss << i+1 <<
416 }
417 return ss.str();
418 }
419
420 //____________________________________________
421 void G4Voxelizer::DisplayListNodes() const
422 {
423 // Prints which solids are present in the sl
424
425 char axis[3] = {'X', 'Y', 'Z'};
426 G4int size=8*sizeof(G4int)*fNPerSlice;
427 G4SurfBits bits(size);
428
429 for (auto j = 0; j <= 2; ++j)
430 {
431 G4cout << " * " << axis[j] << " axis:" <<
432 auto count = (G4int)fBoundaries[j].size();
433 for(G4int i=0; i < count-1; ++i)
434 {
435 G4cout << " Slice #" << i+1 << ": ["
436 << " ; " << fBoundaries[j][i+1] <
437 bits.set(size,(const char *)fBitmasks[j]
438 *fNPerSlice*s
439 G4String result = GetCandidatesAsString(
440 G4cout << "[ " << result.c_str() << "]
441 }
442 }
443 }
444
445 //____________________________________________
446 void G4Voxelizer::BuildBoundingBox()
447 {
448 G4ThreeVector min(fBoundaries[0].front(),
449 fBoundaries[1].front(),
450 fBoundaries[2].front());
451 G4ThreeVector max(fBoundaries[0].back(),
452 fBoundaries[1].back(),
453 fBoundaries[2].back());
454 BuildBoundingBox(min, max);
455 }
456
457 //____________________________________________
458 void G4Voxelizer::BuildBoundingBox(G4ThreeVect
459 G4ThreeVect
460 G4double to
461 {
462 for (auto i = 0; i <= 2; ++i)
463 {
464 G4double min = amin[i];
465 G4double max = amax[i];
466 fBoundingBoxSize[i] = (max - min) / 2 + to
467 fBoundingBoxCenter[i] = min + fBoundingBox
468 }
469 fBoundingBox.SetXHalfLength(fBoundingBoxSize
470 fBoundingBox.SetYHalfLength(fBoundingBoxSize
471 fBoundingBox.SetZHalfLength(fBoundingBoxSize
472 }
473
474 // algorithm -
475
476 // in order to get balanced voxels, merge shou
477 // where the number of voxels is least the num
478 // We will keep sorted list (std::set) with al
479 // comparator function between two voxels, whi
480 // by looking at his right neighbor.
481 // First, we will add all the voxels into the
482 // We will be pick the first item in the tree,
483 // merged voxel into the a list for future red
484 // rebuilded later, therefore they need not to
485 // The merged voxel need to be added to the tr
486 // would be updated.
487
488 //____________________________________________
489 void G4Voxelizer::SetReductionRatio(G4int maxV
490 G4ThreeVec
491 {
492 G4double maxTotal = (G4double) fCandidatesCo
493 * fCandidatesCounts[1].siz
494
495 if (maxVoxels > 0 && maxVoxels < maxTotal)
496 {
497 G4double ratio = (G4double) maxVoxels / ma
498 ratio = std::pow(ratio, 1./3.);
499 if (ratio > 1) { ratio = 1; }
500 reductionRatio.set(ratio,ratio,ratio);
501 }
502 }
503
504 //____________________________________________
505 void G4Voxelizer::BuildReduceVoxels(std::vecto
506 G4ThreeVec
507 {
508 for (auto k = 0; k <= 2; ++k)
509 {
510 std::vector<G4int> &candidatesCount = fCan
511 auto max = (G4int)candidatesCount.size();
512 std::vector<G4VoxelInfo> voxels(max);
513 G4VoxelComparator comp(voxels);
514 std::set<G4int, G4VoxelComparator> voxelSe
515 std::vector<G4int> mergings;
516
517 for (G4int j = 0; j < max; ++j)
518 {
519 G4VoxelInfo &voxel = voxels[j];
520 voxel.count = candidatesCount[j];
521 voxel.previous = j - 1;
522 voxel.next = j + 1;
523 voxels[j] = voxel;
524 }
525
526 for (G4int j = 0; j < max - 1; ++j) { vox
527 // we go to size-1 to make sure we will
528
529 G4double reduction = reductionRatio[k];
530 if (reduction != 0)
531 {
532 G4int count = 0, currentCount;
533 while ((currentCount = (G4int)voxelSet.s
534 {
535 G4double currentRatio = 1 - (G4double)
536 if ((currentRatio <= reduction) && (cu
537 break;
538 const G4int pos = *voxelSet.begin();
539 mergings.push_back(pos + 1);
540
541 G4VoxelInfo& voxel = voxels[pos];
542 G4VoxelInfo& nextVoxel = voxels[voxel.
543
544 if (voxelSet.erase(pos) != 1)
545 {
546 ;// k = k;
547 }
548 if (voxel.next != max - 1)
549 if (voxelSet.erase(voxel.next) != 1)
550 {
551 ;// k = k;
552 }
553 if (voxel.previous != -1)
554 if (voxelSet.erase(voxel.previous) !
555 {
556 ;// k = k;
557 }
558 nextVoxel.count += voxel.count;
559 voxel.count = 0;
560 nextVoxel.previous = voxel.previous;
561
562 if (voxel.next != max - 1)
563 voxelSet.insert(voxel.next);
564
565 if (voxel.previous != -1)
566 {
567 voxels[voxel.previous].next = voxel.
568 voxelSet.insert(voxel.previous);
569 }
570 ++count;
571 } // Loop checking, 13.08.2015, G.Cosmo
572 }
573
574 if (!mergings.empty())
575 {
576 std::sort(mergings.begin(), mergings.end
577
578 const std::vector<G4double>& boundary =
579 auto mergingsSize = (G4int)mergings.size
580 vector<G4double> reducedBoundary;
581 G4int skip = mergings[0], i = 0;
582 max = (G4int)boundary.size();
583 for (G4int j = 0; j < max; ++j)
584 {
585 if (j != skip)
586 {
587 reducedBoundary.push_back(boundary[j
588 }
589 else if (++i < mergingsSize)
590 {
591 skip = mergings[i];
592 }
593 }
594 boundaries[k] = std::move(reducedBoundar
595 }
596 /*
597 G4int count = 0;
598 while (true) // Loop checking, 13.08.20
599 {
600 G4double reduction = reductionRatio[k];
601 if (reduction == 0)
602 break;
603 G4int currentCount = voxelSet.size();
604 if (currentCount <= 2)
605 break;
606 G4double currentRatio = 1 - (G4double) c
607 if (currentRatio <= reduction && current
608 break;
609 const G4int pos = *voxelSet.begin();
610 mergings.push_back(pos);
611
612 G4VoxelInfo &voxel = voxels[pos];
613 G4VoxelInfo &nextVoxel = voxels[voxel.ne
614
615 voxelSet.erase(pos);
616 if (voxel.next != max - 1) { voxelSet.er
617 if (voxel.previous != -1) { voxelSet.er
618
619 nextVoxel.count += voxel.count;
620 voxel.count = 0;
621 nextVoxel.previous = voxel.previous;
622
623 if (voxel.next != max - 1)
624 voxelSet.insert(voxel.next);
625
626 if (voxel.previous != -1)
627 {
628 voxels[voxel.previous].next = voxel.ne
629 voxelSet.insert(voxel.previous);
630 }
631 ++count;
632 }
633
634 if (mergings.size())
635 {
636 std::sort(mergings.begin(), mergings.end
637
638 std::vector<G4double> &boundary = bounda
639 std::vector<G4double> reducedBoundary(bo
640 G4int skip = mergings[0] + 1, cur = 0, i
641 max = boundary.size();
642 for (G4int j = 0; j < max; ++j)
643 {
644 if (j != skip)
645 {
646 reducedBoundary[cur++] = boundary[j]
647 }
648 else
649 {
650 if (++i < (G4int)mergings.size()) {
651 }
652 }
653 boundaries[k] = reducedBoundary;
654 }
655 */
656 }
657 }
658
659 //____________________________________________
660 void G4Voxelizer::BuildReduceVoxels2(std::vect
661 G4ThreeVe
662 {
663 for (auto k = 0; k <= 2; ++k)
664 {
665 std::vector<G4int> &candidatesCount = fCan
666 auto max = (G4int)candidatesCount.size();
667 G4int total = 0;
668 for (G4int i = 0; i < max; ++i) total += c
669
670 G4double reduction = reductionRatio[k];
671 if (reduction == 0)
672 break;
673
674 G4int destination = (G4int) (reduction * m
675 if (destination > 1000) destination = 1000
676 if (destination < 2) destination = 2;
677 G4double average = ((G4double)total / max)
678
679 std::vector<G4int> mergings;
680
681 std::vector<G4double> &boundary = boundari
682 std::vector<G4double> reducedBoundary(dest
683
684 G4int sum = 0, cur = 0;
685 for (G4int i = 0; i < max; ++i)
686 {
687 sum += candidatesCount[i];
688 if (sum > average * (cur + 1) || i == 0)
689 {
690 G4double val = boundary[i];
691 reducedBoundary[cur] = val;
692 ++cur;
693 if (cur == destination)
694 break;
695 }
696 }
697 reducedBoundary[destination-1] = boundary[
698 boundaries[k] = std::move(reducedBoundary)
699 }
700 }
701
702 //____________________________________________
703 void G4Voxelizer::Voxelize(std::vector<G4VSoli
704 std::vector<G4Trans
705 {
706 BuildVoxelLimits(solids, transforms);
707 BuildBoundaries();
708 BuildBitmasks(fBoundaries, fBitmasks);
709 BuildBoundingBox();
710 BuildEmpty(); // this does not work well for
711 // actually only makes perform
712 // these are only pre-calculat
713
714 for (auto & fCandidatesCount : fCandidatesCo
715 {
716 fCandidatesCount.resize(0);
717 }
718 }
719
720 //____________________________________________
721 void G4Voxelizer::CreateMiniVoxels(std::vector
722 G4SurfBits
723 {
724 std::vector<G4int> voxel(3), maxVoxels(3);
725 for (auto i = 0; i <= 2; ++i) maxVoxels[i] =
726
727 G4ThreeVector point;
728 for (voxel[2] = 0; voxel[2] < maxVoxels[2] -
729 {
730 for (voxel[1] = 0; voxel[1] < maxVoxels[1]
731 {
732 for (voxel[0] = 0; voxel[0] < maxVoxels[
733 {
734 std::vector<G4int> candidates;
735 if (GetCandidatesVoxelArray(voxel, bit
736 {
737 // find a box for corresponding non-
738 G4VoxelBox box;
739 for (auto i = 0; i <= 2; ++i)
740 {
741 G4int index = voxel[i];
742 const std::vector<G4double> &bound
743 G4double hlen = 0.5 * (boundary[in
744 box.hlen[i] = hlen;
745 box.pos[i] = boundary[index] + hle
746 }
747 fVoxelBoxes.push_back(box);
748 std::vector<G4int>(candidates).swap(
749 fVoxelBoxesCandidates.push_back(std:
750 }
751 }
752 }
753 }
754 }
755
756 //____________________________________________
757 void G4Voxelizer::Voxelize(std::vector<G4VFace
758 {
759 G4int maxVoxels = fMaxVoxels;
760 G4ThreeVector reductionRatio = fReductionRat
761
762 std::size_t size = facets.size();
763 if (size < 10)
764 {
765 for (const auto & facet : facets)
766 {
767 if (facet->GetNumberOfVertices() > 3) ++
768 }
769 }
770
771 if ((size >= 10 || maxVoxels > 0) && maxVoxe
772 {
773 #ifdef G4SPECSDEBUG
774 G4cout << "Building voxel limits..." << G4
775 #endif
776
777 BuildVoxelLimits(facets);
778
779 #ifdef G4SPECSDEBUG
780 G4cout << "Building boundaries..." << G4en
781 #endif
782
783 BuildBoundaries();
784
785 #ifdef G4SPECSDEBUG
786 G4cout << "Building bitmasks..." << G4endl
787 #endif
788
789 BuildBitmasks(fBoundaries, nullptr, true);
790
791 if (maxVoxels < 0 && reductionRatio == G4T
792 {
793 maxVoxels = fTotalCandidates;
794 if (fTotalCandidates > 1000000) maxVoxel
795 }
796
797 SetReductionRatio(maxVoxels, reductionRati
798
799 fCountOfVoxels = CountVoxels(fBoundaries);
800
801 #ifdef G4SPECSDEBUG
802 G4cout << "Total number of voxels: " << fC
803 #endif
804
805 BuildReduceVoxels2(fBoundaries, reductionR
806
807 fCountOfVoxels = CountVoxels(fBoundaries);
808
809 #ifdef G4SPECSDEBUG
810 G4cout << "Total number of voxels after re
811 << fCountOfVoxels << G4endl;
812 #endif
813
814 #ifdef G4SPECSDEBUG
815 G4cout << "Building bitmasks..." << G4endl
816 #endif
817
818 BuildBitmasks(fBoundaries, fBitmasks);
819
820 G4ThreeVector reductionRatioMini;
821
822 G4SurfBits bitmasksMini[3];
823
824 // section for building mini voxels
825
826 std::vector<G4double> miniBoundaries[3];
827
828 for (auto i = 0; i <= 2; ++i) { miniBound
829
830 G4int voxelsCountMini = (fCountOfVoxels >=
831 ? 100 : G4int(fCount
832
833 SetReductionRatio(voxelsCountMini, reducti
834
835 #ifdef G4SPECSDEBUG
836 G4cout << "Building reduced voxels..." <<
837 #endif
838
839 BuildReduceVoxels(miniBoundaries, reductio
840
841 #ifdef G4SPECSDEBUG
842 G4int total = CountVoxels(miniBoundaries);
843 G4cout << "Total number of mini voxels: "
844 #endif
845
846 #ifdef G4SPECSDEBUG
847 G4cout << "Building mini bitmasks..." << G
848 #endif
849
850 BuildBitmasks(miniBoundaries, bitmasksMini
851
852 #ifdef G4SPECSDEBUG
853 G4cout << "Creating Mini Voxels..." << G4e
854 #endif
855
856 CreateMiniVoxels(miniBoundaries, bitmasksM
857
858 #ifdef G4SPECSDEBUG
859 G4cout << "Building bounding box..." << G4
860 #endif
861
862 BuildBoundingBox();
863
864 #ifdef G4SPECSDEBUG
865 G4cout << "Building empty..." << G4endl;
866 #endif
867
868 BuildEmpty();
869
870 #ifdef G4SPECSDEBUG
871 G4cout << "Deallocating unnecessary fields
872 #endif
873 // deallocate fields unnecessary during ru
874 //
875 fBoxes.resize(0);
876 for (auto i = 0; i < 3; ++i)
877 {
878 fCandidatesCounts[i].resize(0);
879 fBitmasks[i].Clear();
880 }
881 }
882 }
883
884
885 //____________________________________________
886 void G4Voxelizer::GetCandidatesVoxel(std::vect
887 {
888 // "GetCandidates" should compute which soli
889 // the voxel defined by the three slices cha
890
891 G4cout << " Candidates in voxel [" << voxe
892 << " ; " << voxels[2] << "]: ";
893 std::vector<G4int> candidates;
894 G4int count = GetCandidatesVoxelArray(voxels
895 G4cout << "[ ";
896 for (G4int i = 0; i < count; ++i) G4cout <<
897 G4cout << "] " << G4endl;
898 }
899
900 //____________________________________________
901 void G4Voxelizer::FindComponentsFastest(unsign
902 std::
903 {
904 for (G4int byte = 0; byte < (G4int) (sizeof(
905 {
906 if (G4int maskByte = mask & 0xFF)
907 {
908 for (G4int bit = 0; bit < 8; ++bit)
909 {
910 if ((maskByte & 1) != 0)
911 { list.push_back(8*(sizeof(unsigned
912 if ((maskByte >>= 1) == 0) break;
913 }
914 }
915 mask >>= 8;
916 }
917 }
918
919 //____________________________________________
920 void G4Voxelizer::TransformLimits(G4ThreeVecto
921 const G4Tran
922 {
923 // The goal of this method is to convert the
924 // (representing the bounding box of a given
925 // to the main frame, using "transformation"
926
927 G4ThreeVector vertices[8] = // Deteminat
928 { // the exten
929 G4ThreeVector(min.x(), min.y(), min.z()),
930 G4ThreeVector(min.x(), max.y(), min.z()),
931 G4ThreeVector(max.x(), max.y(), min.z()),
932 G4ThreeVector(max.x(), min.y(), min.z()),
933 G4ThreeVector(min.x(), min.y(), max.z()),
934 G4ThreeVector(min.x(), max.y(), max.z()),
935 G4ThreeVector(max.x(), max.y(), max.z()),
936 G4ThreeVector(max.x(), min.y(), max.z())
937 };
938
939 min.set(kInfinity,kInfinity,kInfinity);
940 max.set(-kInfinity,-kInfinity,-kInfinity);
941
942 // Loop on th vertices
943 G4int limit = sizeof(vertices) / sizeof(G4Th
944 for (G4int i = 0 ; i < limit; ++i)
945 {
946 // From local frame to the global one:
947 // Current positions on the three axis:
948 G4ThreeVector current = GetGlobalPoint(tra
949
950 // If need be, replacement of the min & ma
951 if (current.x() > max.x()) max.setX(curren
952 if (current.x() < min.x()) min.setX(curren
953
954 if (current.y() > max.y()) max.setY(curren
955 if (current.y() < min.y()) min.setY(curren
956
957 if (current.z() > max.z()) max.setZ(curren
958 if (current.z() < min.z()) min.setZ(curren
959 }
960 }
961
962 //____________________________________________
963 G4int G4Voxelizer::GetCandidatesVoxelArray(con
964 std::vector<G4int> &
965 {
966 // Method returning the candidates correspon
967
968 list.clear();
969
970 for (auto i = 0; i <= 2; ++i)
971 {
972 if(point[i] < fBoundaries[i].front() || po
973 return 0;
974 }
975
976 if (fTotalCandidates == 1)
977 {
978 list.push_back(0);
979 return 1;
980 }
981 else
982 {
983 if (fNPerSlice == 1)
984 {
985 unsigned int mask = 0xFFffFFff;
986 G4int slice;
987 if (fBoundaries[0].size() > 2)
988 {
989 slice = BinarySearch(fBoundaries[0], p
990 if ((mask = ((unsigned int*) fBitmasks
991 return 0;
992 }
993 if (fBoundaries[1].size() > 2)
994 {
995 slice = BinarySearch(fBoundaries[1], p
996 if ((mask &= ((unsigned int*) fBitmask
997 return 0;
998 }
999 if (fBoundaries[2].size() > 2)
1000 {
1001 slice = BinarySearch(fBoundaries[2],
1002 if ((mask &= ((unsigned int*) fBitmas
1003 return 0;
1004 }
1005 if ((crossed != nullptr) && ((mask &= ~
1006 return 0;
1007
1008 FindComponentsFastest(mask, list, 0);
1009 }
1010 else
1011 {
1012 unsigned int* masks[3], mask; // masks
1013 for (auto i = 0; i <= 2; ++i)
1014 {
1015 G4int slice = BinarySearch(fBoundarie
1016 masks[i] = ((unsigned int*) fBitmasks
1017 + slice * fNPerSlice;
1018 }
1019 unsigned int* maskCrossed = crossed !=
1020 ? (unsigned i
1021
1022 for (G4int i = 0 ; i < fNPerSlice; ++i)
1023 {
1024 // Logic "and" of the masks along the
1025 // removing "if (!" and ") continue"
1026 //
1027 if ((mask = masks[0][i]) == 0u) conti
1028 if ((mask &= masks[1][i]) == 0u) cont
1029 if ((mask &= masks[2][i]) == 0u) cont
1030 if ((maskCrossed != nullptr) && ((mas
1031
1032 FindComponentsFastest(mask, list, i);
1033 }
1034 }
1035 /*
1036 if (fNPerSlice == 1)
1037 {
1038 unsigned int mask;
1039 G4int slice = BinarySearch(fBoundaries[
1040 if (!(mask = ((unsigned int *) fBitmask
1041 )) return 0;
1042 slice = BinarySearch(fBoundaries[1], po
1043 if (!(mask &= ((unsigned int *) fBitmas
1044 )) return 0;
1045 slice = BinarySearch(fBoundaries[2], po
1046 if (!(mask &= ((unsigned int *) fBitmas
1047 )) return 0;
1048 if (crossed && (!(mask &= ~((unsigned i
1049 return 0;
1050
1051 FindComponentsFastest(mask, list, 0);
1052 }
1053 else
1054 {
1055 unsigned int *masks[3], mask; // masks
1056 for (auto i = 0; i <= 2; ++i)
1057 {
1058 G4int slice = BinarySearch(fBoundarie
1059 masks[i] = ((unsigned int *) fBitmask
1060 }
1061 unsigned int *maskCrossed =
1062 crossed ? (unsigned int *)crossed->fA
1063
1064 for (G4int i = 0 ; i < fNPerSlice; ++i)
1065 {
1066 // Logic "and" of the masks along the
1067 // removing "if (!" and ") continue"
1068 //
1069 if (!(mask = masks[0][i])) continue;
1070 if (!(mask &= masks[1][i])) continue;
1071 if (!(mask &= masks[2][i])) continue;
1072 if (maskCrossed && !(mask &= ~maskCro
1073
1074 FindComponentsFastest(mask, list, i);
1075 }
1076 }
1077 */
1078 }
1079 return (G4int)list.size();
1080 }
1081
1082 //___________________________________________
1083 G4int
1084 G4Voxelizer::GetCandidatesVoxelArray(const st
1085 const G4
1086 std::vec
1087 G4SurfBi
1088 {
1089 list.clear();
1090
1091 if (fTotalCandidates == 1)
1092 {
1093 list.push_back(0);
1094 return 1;
1095 }
1096 else
1097 {
1098 if (fNPerSlice == 1)
1099 {
1100 unsigned int mask;
1101 if ((mask = ((unsigned int *) bitmasks[
1102 return 0;
1103 if ((mask &= ((unsigned int *) bitmasks
1104 return 0;
1105 if ((mask &= ((unsigned int *) bitmasks
1106 return 0;
1107 if ((crossed != nullptr) && ((mask &= ~
1108 return 0;
1109
1110 FindComponentsFastest(mask, list, 0);
1111 }
1112 else
1113 {
1114 unsigned int *masks[3], mask; // masks
1115 for (auto i = 0; i <= 2; ++i)
1116 {
1117 masks[i] = ((unsigned int *) bitmasks
1118 + voxels[i]*fNPerSlice;
1119 }
1120 unsigned int *maskCrossed = crossed !=
1121 ? (unsigned i
1122
1123 for (G4int i = 0 ; i < fNPerSlice; ++i)
1124 {
1125 // Logic "and" of the masks along the
1126 // removing "if (!" and ") continue"
1127 //
1128 if ((mask = masks[0][i]) == 0u) conti
1129 if ((mask &= masks[1][i]) == 0u) cont
1130 if ((mask &= masks[2][i]) == 0u) cont
1131 if ((maskCrossed != nullptr) && ((mas
1132
1133 FindComponentsFastest(mask, list, i);
1134 }
1135 }
1136 }
1137 return (G4int)list.size();
1138 }
1139
1140 //___________________________________________
1141 G4int
1142 G4Voxelizer::GetCandidatesVoxelArray(const st
1143 std::vector<G4int>&
1144 {
1145 // Method returning the candidates correspo
1146
1147 return GetCandidatesVoxelArray(voxels, fBit
1148 }
1149
1150 //___________________________________________
1151 G4bool G4Voxelizer::Contains(const G4ThreeVec
1152 {
1153 for (auto i = 0; i < 3; ++i)
1154 {
1155 if (point[i] < fBoundaries[i].front() ||
1156 return false;
1157 }
1158 return true;
1159 }
1160
1161 //___________________________________________
1162 G4double
1163 G4Voxelizer::DistanceToFirst(const G4ThreeVec
1164 const G4ThreeVec
1165 {
1166 G4ThreeVector pointShifted = point - fBound
1167 G4double shift = fBoundingBox.DistanceToIn(
1168 return shift;
1169 }
1170
1171 //___________________________________________
1172 G4double
1173 G4Voxelizer::DistanceToBoundingBox(const G4Th
1174 {
1175 G4ThreeVector pointShifted = point - fBound
1176 G4double shift = MinDistanceToBox(pointShif
1177 return shift;
1178 }
1179
1180 //___________________________________________
1181 G4double
1182 G4Voxelizer::MinDistanceToBox (const G4ThreeV
1183 const G4ThreeV
1184 {
1185 // Estimates the isotropic safety from a po
1186 // any of its surfaces. The algorithm may b
1187 // fast underestimate.
1188
1189 G4double safe, safx, safy, safz;
1190 safe = safx = -f.x() + std::abs(aPoint.x())
1191 safy = -f.y() + std::abs(aPoint.y());
1192 if ( safy > safe ) safe = safy;
1193 safz = -f.z() + std::abs(aPoint.z());
1194 if ( safz > safe ) safe = safz;
1195 if (safe < 0.0) return 0.0; // point is ins
1196
1197 G4double safsq = 0.0;
1198 G4int count = 0;
1199 if ( safx > 0 ) { safsq += safx*safx; ++cou
1200 if ( safy > 0 ) { safsq += safy*safy; ++cou
1201 if ( safz > 0 ) { safsq += safz*safz; ++cou
1202 if (count == 1) return safe;
1203 return std::sqrt(safsq);
1204 }
1205
1206 //___________________________________________
1207 G4double
1208 G4Voxelizer::DistanceToNext(const G4ThreeVect
1209 const G4ThreeVect
1210 std::vector
1211 {
1212 G4double shift = kInfinity;
1213
1214 G4int cur = 0; // the smallest index, which
1215 for (G4int i = 0; i <= 2; ++i)
1216 {
1217 // Looking for the next voxels on the con
1218 //
1219 const std::vector<G4double>& boundary = f
1220 G4int index = curVoxel[i];
1221 if (direction[i] >= 1e-10)
1222 {
1223 ++index;
1224 }
1225 else
1226 {
1227 if (direction[i] > -1e-10)
1228 continue;
1229 }
1230 G4double dif = boundary[index] - point[i]
1231 G4double distance = dif / direction[i];
1232
1233 if (shift > distance)
1234 {
1235 shift = distance;
1236 cur = i;
1237 }
1238 }
1239
1240 if (shift != kInfinity)
1241 {
1242 // updating current voxel using the index
1243 // to the closest voxel boundary on the r
1244
1245 if (direction[cur] > 0)
1246 {
1247 if (++curVoxel[cur] >= (G4int) fBoundar
1248 shift = kInfinity;
1249 }
1250 else
1251 {
1252 if (--curVoxel[cur] < 0)
1253 shift = kInfinity;
1254 }
1255 }
1256
1257 /*
1258 for (auto i = 0; i <= 2; ++i)
1259 {
1260 // Looking for the next voxels on the con
1261 //
1262 const std::vector<G4double> &boundary = f
1263 G4int cur = curVoxel[i];
1264 if(direction[i] >= 1e-10)
1265 {
1266 if (boundary[++cur] - point[i] < fTol
1267 if (++cur >= (G4int) boundary.size())
1268 continue;
1269 }
1270 else
1271 {
1272 if(direction[i] <= -1e-10)
1273 {
1274 if (point[i] - boundary[cur] < fToler
1275 if (--cur < 0)
1276 continue;
1277 }
1278 else
1279 continue;
1280 }
1281 G4double dif = boundary[cur] - point[i];
1282 G4double distance = dif / direction[i];
1283
1284 if (shift > distance)
1285 shift = distance;
1286 }
1287 */
1288 return shift;
1289 }
1290
1291 //___________________________________________
1292 G4bool
1293 G4Voxelizer::UpdateCurrentVoxel(const G4Three
1294 const G4Three
1295 std::vector<G
1296 {
1297 for (auto i = 0; i <= 2; ++i)
1298 {
1299 G4int index = curVoxel[i];
1300 const std::vector<G4double> &boundary = f
1301
1302 if (direction[i] > 0)
1303 {
1304 if (point[i] >= boundary[++index])
1305 if (++curVoxel[i] >= (G4int) boundary
1306 return false;
1307 }
1308 else
1309 {
1310 if (point[i] < boundary[index])
1311 if (--curVoxel[i] < 0)
1312 return false;
1313 }
1314 #ifdef G4SPECSDEBUG
1315 G4int indexOK = BinarySearch(boundary, po
1316 if (curVoxel[i] != indexOK)
1317 curVoxel[i] = indexOK; // put breakpoin
1318 #endif
1319 }
1320 return true;
1321 }
1322
1323 //___________________________________________
1324 void G4Voxelizer::SetMaxVoxels(G4int max)
1325 {
1326 fMaxVoxels = max;
1327 fReductionRatio.set(0,0,0);
1328 }
1329
1330 //___________________________________________
1331 void G4Voxelizer::SetMaxVoxels(const G4ThreeV
1332 {
1333 fMaxVoxels = -1;
1334 fReductionRatio = ratioOfReduction;
1335 }
1336
1337 //___________________________________________
1338 void G4Voxelizer::SetDefaultVoxelsCount(G4int
1339 {
1340 fDefaultVoxelsCount = count;
1341 }
1342
1343 //___________________________________________
1344 G4int G4Voxelizer::GetDefaultVoxelsCount()
1345 {
1346 return fDefaultVoxelsCount;
1347 }
1348
1349 //___________________________________________
1350 G4int G4Voxelizer::AllocatedMemory()
1351 {
1352 G4int size = fEmpty.GetNbytes();
1353 size += fBoxes.capacity() * sizeof(G4VoxelB
1354 size += sizeof(G4double) * (fBoundaries[0].
1355 + fBoundaries[1].capacity() + fBounda
1356 size += sizeof(G4int) * (fCandidatesCounts[
1357 + fCandidatesCounts[1].capacity() + f
1358 size += fBitmasks[0].GetNbytes() + fBitmask
1359 + fBitmasks[2].GetNbytes();
1360
1361 auto csize = (G4int)fCandidates.size();
1362 for (G4int i = 0; i < csize; ++i)
1363 {
1364 size += sizeof(vector<G4int>) + fCandidat
1365 }
1366
1367 return size;
1368 }
1369