| Geant4 Cross Reference |
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 // 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::fDefaultVoxelsCount = -1;
52
53 //______________________________________________________________________________
54 G4Voxelizer::G4Voxelizer()
55 : fBoundingBox("VoxBBox", 1, 1, 1)
56 {
57 fCountOfVoxels = fNPerSlice = fTotalCandidates = 0;
58
59 fTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
60
61 SetMaxVoxels(fDefaultVoxelsCount);
62
63 G4SolidStore::GetInstance()->DeRegister(&fBoundingBox);
64 }
65
66 //______________________________________________________________________________
67 G4Voxelizer::~G4Voxelizer() = default;
68
69 //______________________________________________________________________________
70 void G4Voxelizer::BuildEmpty()
71 {
72 // by reserving the size of candidates, we would avoid reallocation of
73 // the vector which could cause fragmentation
74 //
75 std::vector<G4int> xyz(3), max(3), candidates(fTotalCandidates);
76 const std::vector<G4int> empty(0);
77
78 for (auto i = 0; i <= 2; ++i) max[i] = (G4int)fBoundaries[i].size();
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[0])
90 {
91 if (GetCandidatesVoxelArray(xyz, candidates) != 0)
92 {
93 G4int index = GetVoxelsIndex(xyz);
94 fEmpty.SetBitNumber(index, false);
95
96 // rather than assigning directly with:
97 // "fCandidates[index] = candidates;", in an effort to ensure that
98 // capacity would be just exact, we rather use following 3 lines
99 //
100 std::vector<G4int> &c = (fCandidates[index] = empty);
101 c.reserve(candidates.size());
102 c.assign(candidates.begin(), candidates.end());
103 }
104 }
105 }
106 }
107 #ifdef G4SPECSDEBUG
108 G4cout << "Non-empty voxels count: " << fCandidates.size() << G4endl;
109 #endif
110 }
111
112 //______________________________________________________________________________
113 void G4Voxelizer::BuildVoxelLimits(std::vector<G4VSolid*>& solids,
114 std::vector<G4Transform3D>& transforms)
115 {
116 // "BuildVoxelLimits"'s aim is to store the coordinates of the origin as
117 // well as the half lengths related to the bounding box of each node.
118 // These quantities are stored in the array "fBoxes" (6 different values per
119 // node
120 //
121 if (std::size_t numNodes = solids.size()) // Number of nodes in "multiUnion"
122 {
123 fBoxes.resize(numNodes); // Array which will store the half lengths
124 fNPerSlice = G4int(1 + (fBoxes.size() - 1) / (8 * sizeof(unsigned int)));
125
126 // related to a particular node, but also
127 // the coordinates of its origin
128
129 G4ThreeVector toleranceVector(fTolerance,fTolerance,fTolerance);
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.GetRadialTolerance() / 2.0;
143 orbToleranceVector.set(tolerance,tolerance,tolerance);
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<G4VFacet*>& facets)
162 {
163 // "BuildVoxelLimits"'s aim is to store the coordinates of the origin as well
164 // as the half lengths related to the bounding box of each node.
165 // These quantities are stored in the array "fBoxes" (6 different values per
166 // node.
167
168 if (std::size_t numNodes = facets.size()) // Number of nodes
169 {
170 fBoxes.resize(numNodes); // Array which will store the half lengths
171 fNPerSlice = G4int(1+(fBoxes.size()-1)/(8*sizeof(unsigned int)));
172
173 G4ThreeVector toleranceVector(10*fTolerance, 10*fTolerance, 10*fTolerance);
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,1);
180 G4ThreeVector extent;
181 max.set (facet.Extent(x), facet.Extent(y), facet.Extent(z));
182 min.set (-facet.Extent(-x), -facet.Extent(-y), -facet.Extent(-z));
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, dZ, pX, pY and pZ for each node
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::fixed) <<
203 " -> Node " << i+1 << ":\n" <<
204 "\t * [x,y,z] = " << fBoxes[i].hlen <<
205 "\t * [x,y,z] = " << fBoxes[i].pos << "\n";
206 }
207 G4cout.precision(oldprec);
208 }
209
210 //______________________________________________________________________________
211 void G4Voxelizer::CreateSortedBoundary(std::vector<G4double>& boundary,
212 G4int axis)
213 {
214 // "CreateBoundaries"'s aim is to determine the slices induced by the
215 // bounding fBoxes, along each axis. The created boundaries are stored
216 // in the array "boundariesRaw"
217
218 std::size_t numNodes = fBoxes.size(); // Number of nodes in structure
219
220 // Determination of the boundaries along x, y and z axis
221 //
222 for(std::size_t i = 0 ; i < numNodes; ++i)
223 {
224 // For each node, the boundaries are created by using the array "fBoxes"
225 // built in method "BuildVoxelLimits"
226 //
227 G4double p = fBoxes[i].pos[axis], d = fBoxes[i].hlen[axis];
228
229 // x boundaries
230 //
231 #ifdef G4SPECSDEBUG
232 G4cout << "Boundary " << p - d << " - " << p + d << G4endl;
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 along each axis (increasing order)
244 // and also does not take into account redundant boundaries, i.e. if two
245 // boundaries are separated by a distance strictly inferior to "tolerance".
246 // The sorted boundaries are respectively stored in:
247 // * boundaries[0..2]
248
249 // In addition, the number of elements contained in the three latter arrays
250 // are precise thanks to variables: boundariesCountX, boundariesCountY and
251 // boundariesCountZ.
252
253 if (std::size_t numNodes = fBoxes.size())
254 {
255 const G4double tolerance = fTolerance / 100.0;
256 // Minimal distance to discriminate two boundaries.
257
258 std::vector<G4double> sortedBoundary(2*numNodes);
259
260 for (auto j = 0; j <= 2; ++j)
261 {
262 CreateSortedBoundary(sortedBoundary, j);
263 std::vector<G4double> &boundary = fBoundaries[j];
264 boundary.clear();
265
266 for(std::size_t i = 0 ; i < 2*numNodes; ++i)
267 {
268 G4double newBoundary = sortedBoundary[i];
269 #ifdef G4SPECSDEBUG
270 if (j == 0) G4cout << "Examining " << newBoundary << "..." << G4endl;
271 #endif
272 auto size = (G4int)boundary.size();
273 if((size == 0) || std::abs(boundary[size-1] - newBoundary) > tolerance)
274 {
275 {
276 #ifdef G4SPECSDEBUG
277 if (j == 0) G4cout << "Adding boundary " << newBoundary << "..."
278 << G4endl;
279 #endif
280 boundary.push_back(newBoundary);
281 continue;
282 }
283 }
284 // If two successive boundaries are too close from each other,
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 be 2x bigger then 50.000.
293 // therefore only from 100.000 reduced
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 == size - 1)
300 {
301 // this condition of merging boundaries was wrong,
302 // it did not count with right part, which can be
303 // completely ommited and not included in final consideration.
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:" << G4endl << " | ";
322 DisplayBoundaries(fBoundaries[i]);
323 }
324 }
325
326 //______________________________________________________________________________
327 void G4Voxelizer::DisplayBoundaries(std::vector<G4double> &boundaries)
328 {
329 // Prints the positions of the boundaries of the slices on the three axes
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::fixed) << boundaries[i];
336 if(i != count-1) G4cout << "-> ";
337 }
338 G4cout << "|" << G4endl << "Number of boundaries: " << count << G4endl;
339 G4cout.precision(oldprec);
340 }
341
342 //______________________________________________________________________________
343 void G4Voxelizer::BuildBitmasks(std::vector<G4double> boundaries[],
344 G4SurfBits bitmasks[], G4bool countsOnly)
345 {
346 // "BuildListNodes" stores in the bitmasks solids present in each slice
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 = boundaries[k];
355 G4int voxelsCount = (G4int)boundary.size() - 1;
356 G4SurfBits& bitmask = bitmasks[k];
357
358 if (!countsOnly)
359 {
360 bitmask.Clear();
361 #ifdef G4SPECSDEBUG
362 G4cout << "Allocating bitmask..." << G4endl;
363 #endif
364 bitmask.SetBitNumber(voxelsCount*bitsPerSlice-1, false);
365 // it is here so we can set the maximum number of bits. this line
366 // will rellocate the memory and set all to zero
367 }
368 std::vector<G4int>& candidatesCount = fCandidatesCounts[k];
369 candidatesCount.resize(voxelsCount);
370
371 for(G4int i = 0 ; i < voxelsCount; ++i) { candidatesCount[i] = 0; }
372
373 // Loop on the nodes, number of slices per axis
374 //
375 for(std::size_t j = 0 ; j < numNodes; ++j)
376 {
377 // Determination of the minimum and maximum position along x
378 // of the bounding boxe of each node
379 //
380 G4double p = fBoxes[j].pos[k], d = fBoxes[j].hlen[k];
381
382 G4double min = p - d; // - localTolerance;
383 G4double max = p + d; // + localTolerance;
384
385 G4int i = BinarySearch(boundary, min);
386 if (i < 0) { i = 0; }
387
388 do // Loop checking, 13.08.2015, G.Cosmo
389 {
390 if (!countsOnly)
391 {
392 bitmask.SetBitNumber(i*bitsPerSlice+(G4int)j);
393 }
394 candidatesCount[i]++;
395 ++i;
396 }
397 while (max > boundary[i] && i < voxelsCount);
398 }
399 }
400 #ifdef G4SPECSDEBUG
401 G4cout << "Build list nodes completed." << G4endl;
402 #endif
403 }
404
405 //______________________________________________________________________________
406 G4String G4Voxelizer::GetCandidatesAsString(const G4SurfBits& bits) const
407 {
408 // Decodes the candidates in mask as G4String.
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 slices previously elaborated.
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:" << G4endl;
432 auto count = (G4int)fBoundaries[j].size();
433 for(G4int i=0; i < count-1; ++i)
434 {
435 G4cout << " Slice #" << i+1 << ": [" << fBoundaries[j][i]
436 << " ; " << fBoundaries[j][i+1] << "] -> ";
437 bits.set(size,(const char *)fBitmasks[j].fAllBits+i
438 *fNPerSlice*sizeof(G4int));
439 G4String result = GetCandidatesAsString(bits);
440 G4cout << "[ " << result.c_str() << "] " << G4endl;
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(G4ThreeVector& amin,
459 G4ThreeVector& amax,
460 G4double tolerance)
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 + tolerance * 0.5;
467 fBoundingBoxCenter[i] = min + fBoundingBoxSize[i];
468 }
469 fBoundingBox.SetXHalfLength(fBoundingBoxSize.x());
470 fBoundingBox.SetYHalfLength(fBoundingBoxSize.y());
471 fBoundingBox.SetZHalfLength(fBoundingBoxSize.z());
472 }
473
474 // algorithm -
475
476 // in order to get balanced voxels, merge should always unite those regions,
477 // where the number of voxels is least the number.
478 // We will keep sorted list (std::set) with all voxels. there will be
479 // comparator function between two voxels, which will tell if voxel is less
480 // by looking at his right neighbor.
481 // First, we will add all the voxels into the tree.
482 // We will be pick the first item in the tree, merging it, adding the right
483 // merged voxel into the a list for future reduction (fBitmasks will be
484 // rebuilded later, therefore they need not to be updated).
485 // The merged voxel need to be added to the tree again, so it's position
486 // would be updated.
487
488 //______________________________________________________________________________
489 void G4Voxelizer::SetReductionRatio(G4int maxVoxels,
490 G4ThreeVector& reductionRatio)
491 {
492 G4double maxTotal = (G4double) fCandidatesCounts[0].size()
493 * fCandidatesCounts[1].size() * fCandidatesCounts[2].size();
494
495 if (maxVoxels > 0 && maxVoxels < maxTotal)
496 {
497 G4double ratio = (G4double) maxVoxels / maxTotal;
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::vector<G4double> boundaries[],
506 G4ThreeVector reductionRatio)
507 {
508 for (auto k = 0; k <= 2; ++k)
509 {
510 std::vector<G4int> &candidatesCount = fCandidatesCounts[k];
511 auto max = (G4int)candidatesCount.size();
512 std::vector<G4VoxelInfo> voxels(max);
513 G4VoxelComparator comp(voxels);
514 std::set<G4int, G4VoxelComparator> voxelSet(comp);
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) { voxelSet.insert(j); }
527 // we go to size-1 to make sure we will not merge the last element
528
529 G4double reduction = reductionRatio[k];
530 if (reduction != 0)
531 {
532 G4int count = 0, currentCount;
533 while ((currentCount = (G4int)voxelSet.size()) > 2)
534 {
535 G4double currentRatio = 1 - (G4double) count / max;
536 if ((currentRatio <= reduction) && (currentCount <= 1000))
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.next];
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) != 1)
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.next;
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 = boundaries[k];
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(reducedBoundary);
595 }
596 /*
597 G4int count = 0;
598 while (true) // Loop checking, 13.08.2015, G.Cosmo
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) count / max;
607 if (currentRatio <= reduction && currentCount <= 1000)
608 break;
609 const G4int pos = *voxelSet.begin();
610 mergings.push_back(pos);
611
612 G4VoxelInfo &voxel = voxels[pos];
613 G4VoxelInfo &nextVoxel = voxels[voxel.next];
614
615 voxelSet.erase(pos);
616 if (voxel.next != max - 1) { voxelSet.erase(voxel.next); }
617 if (voxel.previous != -1) { voxelSet.erase(voxel.previous); }
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.next;
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 = boundaries[k];
639 std::vector<G4double> reducedBoundary(boundary.size() - mergings.size());
640 G4int skip = mergings[0] + 1, cur = 0, i = 0;
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()) { skip = mergings[i] + 1; }
651 }
652 }
653 boundaries[k] = reducedBoundary;
654 }
655 */
656 }
657 }
658
659 //______________________________________________________________________________
660 void G4Voxelizer::BuildReduceVoxels2(std::vector<G4double> boundaries[],
661 G4ThreeVector reductionRatio)
662 {
663 for (auto k = 0; k <= 2; ++k)
664 {
665 std::vector<G4int> &candidatesCount = fCandidatesCounts[k];
666 auto max = (G4int)candidatesCount.size();
667 G4int total = 0;
668 for (G4int i = 0; i < max; ++i) total += candidatesCount[i];
669
670 G4double reduction = reductionRatio[k];
671 if (reduction == 0)
672 break;
673
674 G4int destination = (G4int) (reduction * max) + 1;
675 if (destination > 1000) destination = 1000;
676 if (destination < 2) destination = 2;
677 G4double average = ((G4double)total / max) / reduction;
678
679 std::vector<G4int> mergings;
680
681 std::vector<G4double> &boundary = boundaries[k];
682 std::vector<G4double> reducedBoundary(destination);
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[max];
698 boundaries[k] = std::move(reducedBoundary);
699 }
700 }
701
702 //______________________________________________________________________________
703 void G4Voxelizer::Voxelize(std::vector<G4VSolid*>& solids,
704 std::vector<G4Transform3D>& transforms)
705 {
706 BuildVoxelLimits(solids, transforms);
707 BuildBoundaries();
708 BuildBitmasks(fBoundaries, fBitmasks);
709 BuildBoundingBox();
710 BuildEmpty(); // this does not work well for multi-union,
711 // actually only makes performance slower,
712 // these are only pre-calculated but not used by multi-union
713
714 for (auto & fCandidatesCount : fCandidatesCounts)
715 {
716 fCandidatesCount.resize(0);
717 }
718 }
719
720 //______________________________________________________________________________
721 void G4Voxelizer::CreateMiniVoxels(std::vector<G4double> boundaries[],
722 G4SurfBits bitmasks[])
723 {
724 std::vector<G4int> voxel(3), maxVoxels(3);
725 for (auto i = 0; i <= 2; ++i) maxVoxels[i] = (G4int)boundaries[i].size();
726
727 G4ThreeVector point;
728 for (voxel[2] = 0; voxel[2] < maxVoxels[2] - 1; ++voxel[2])
729 {
730 for (voxel[1] = 0; voxel[1] < maxVoxels[1] - 1; ++voxel[1])
731 {
732 for (voxel[0] = 0; voxel[0] < maxVoxels[0] - 1; ++voxel[0])
733 {
734 std::vector<G4int> candidates;
735 if (GetCandidatesVoxelArray(voxel, bitmasks, candidates, nullptr) != 0)
736 {
737 // find a box for corresponding non-empty voxel
738 G4VoxelBox box;
739 for (auto i = 0; i <= 2; ++i)
740 {
741 G4int index = voxel[i];
742 const std::vector<G4double> &boundary = boundaries[i];
743 G4double hlen = 0.5 * (boundary[index+1] - boundary[index]);
744 box.hlen[i] = hlen;
745 box.pos[i] = boundary[index] + hlen;
746 }
747 fVoxelBoxes.push_back(box);
748 std::vector<G4int>(candidates).swap(candidates);
749 fVoxelBoxesCandidates.push_back(std::move(candidates));
750 }
751 }
752 }
753 }
754 }
755
756 //______________________________________________________________________________
757 void G4Voxelizer::Voxelize(std::vector<G4VFacet*>& facets)
758 {
759 G4int maxVoxels = fMaxVoxels;
760 G4ThreeVector reductionRatio = fReductionRatio;
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) ++size;
768 }
769 }
770
771 if ((size >= 10 || maxVoxels > 0) && maxVoxels != 0 && maxVoxels != 1)
772 {
773 #ifdef G4SPECSDEBUG
774 G4cout << "Building voxel limits..." << G4endl;
775 #endif
776
777 BuildVoxelLimits(facets);
778
779 #ifdef G4SPECSDEBUG
780 G4cout << "Building boundaries..." << G4endl;
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 == G4ThreeVector())
792 {
793 maxVoxels = fTotalCandidates;
794 if (fTotalCandidates > 1000000) maxVoxels = 1000000;
795 }
796
797 SetReductionRatio(maxVoxels, reductionRatio);
798
799 fCountOfVoxels = CountVoxels(fBoundaries);
800
801 #ifdef G4SPECSDEBUG
802 G4cout << "Total number of voxels: " << fCountOfVoxels << G4endl;
803 #endif
804
805 BuildReduceVoxels2(fBoundaries, reductionRatio);
806
807 fCountOfVoxels = CountVoxels(fBoundaries);
808
809 #ifdef G4SPECSDEBUG
810 G4cout << "Total number of voxels after reduction: "
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) { miniBoundaries[i] = fBoundaries[i]; }
829
830 G4int voxelsCountMini = (fCountOfVoxels >= 1000)
831 ? 100 : G4int(fCountOfVoxels / 10);
832
833 SetReductionRatio(voxelsCountMini, reductionRatioMini);
834
835 #ifdef G4SPECSDEBUG
836 G4cout << "Building reduced voxels..." << G4endl;
837 #endif
838
839 BuildReduceVoxels(miniBoundaries, reductionRatioMini);
840
841 #ifdef G4SPECSDEBUG
842 G4int total = CountVoxels(miniBoundaries);
843 G4cout << "Total number of mini voxels: " << total << G4endl;
844 #endif
845
846 #ifdef G4SPECSDEBUG
847 G4cout << "Building mini bitmasks..." << G4endl;
848 #endif
849
850 BuildBitmasks(miniBoundaries, bitmasksMini);
851
852 #ifdef G4SPECSDEBUG
853 G4cout << "Creating Mini Voxels..." << G4endl;
854 #endif
855
856 CreateMiniVoxels(miniBoundaries, bitmasksMini);
857
858 #ifdef G4SPECSDEBUG
859 G4cout << "Building bounding box..." << G4endl;
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 during runtime..." << G4endl;
872 #endif
873 // deallocate fields unnecessary during runtime
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::vector<G4int>& voxels)
887 {
888 // "GetCandidates" should compute which solids are possibly contained in
889 // the voxel defined by the three slices characterized by the passed indexes.
890
891 G4cout << " Candidates in voxel [" << voxels[0] << " ; " << voxels[1]
892 << " ; " << voxels[2] << "]: ";
893 std::vector<G4int> candidates;
894 G4int count = GetCandidatesVoxelArray(voxels, candidates);
895 G4cout << "[ ";
896 for (G4int i = 0; i < count; ++i) G4cout << candidates[i];
897 G4cout << "] " << G4endl;
898 }
899
900 //______________________________________________________________________________
901 void G4Voxelizer::FindComponentsFastest(unsigned int mask,
902 std::vector<G4int>& list, G4int i)
903 {
904 for (G4int byte = 0; byte < (G4int) (sizeof(unsigned int)); ++byte)
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 int)*i+ byte) + bit); }
912 if ((maskByte >>= 1) == 0) break;
913 }
914 }
915 mask >>= 8;
916 }
917 }
918
919 //______________________________________________________________________________
920 void G4Voxelizer::TransformLimits(G4ThreeVector& min, G4ThreeVector& max,
921 const G4Transform3D& transformation) const
922 {
923 // The goal of this method is to convert the quantities min and max
924 // (representing the bounding box of a given solid in its local frame)
925 // to the main frame, using "transformation"
926
927 G4ThreeVector vertices[8] = // Detemination of the vertices thanks to
928 { // the extension of each solid:
929 G4ThreeVector(min.x(), min.y(), min.z()), // 1st vertice:
930 G4ThreeVector(min.x(), max.y(), min.z()), // 2nd vertice:
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(G4ThreeVector);
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(transformation, vertices[i]);
949
950 // If need be, replacement of the min & max values:
951 if (current.x() > max.x()) max.setX(current.x());
952 if (current.x() < min.x()) min.setX(current.x());
953
954 if (current.y() > max.y()) max.setY(current.y());
955 if (current.y() < min.y()) min.setY(current.y());
956
957 if (current.z() > max.z()) max.setZ(current.z());
958 if (current.z() < min.z()) min.setZ(current.z());
959 }
960 }
961
962 //______________________________________________________________________________
963 G4int G4Voxelizer::GetCandidatesVoxelArray(const G4ThreeVector &point,
964 std::vector<G4int> &list, G4SurfBits *crossed) const
965 {
966 // Method returning the candidates corresponding to the passed point
967
968 list.clear();
969
970 for (auto i = 0; i <= 2; ++i)
971 {
972 if(point[i] < fBoundaries[i].front() || point[i] >= fBoundaries[i].back())
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], point.x());
990 if ((mask = ((unsigned int*) fBitmasks[0].fAllBits)[slice]) == 0u)
991 return 0;
992 }
993 if (fBoundaries[1].size() > 2)
994 {
995 slice = BinarySearch(fBoundaries[1], point.y());
996 if ((mask &= ((unsigned int*) fBitmasks[1].fAllBits)[slice]) == 0u)
997 return 0;
998 }
999 if (fBoundaries[2].size() > 2)
1000 {
1001 slice = BinarySearch(fBoundaries[2], point.z());
1002 if ((mask &= ((unsigned int*) fBitmasks[2].fAllBits)[slice]) == 0u)
1003 return 0;
1004 }
1005 if ((crossed != nullptr) && ((mask &= ~((unsigned int*)crossed->fAllBits)[0]) == 0u))
1006 return 0;
1007
1008 FindComponentsFastest(mask, list, 0);
1009 }
1010 else
1011 {
1012 unsigned int* masks[3], mask; // masks for X,Y,Z axis
1013 for (auto i = 0; i <= 2; ++i)
1014 {
1015 G4int slice = BinarySearch(fBoundaries[i], point[i]);
1016 masks[i] = ((unsigned int*) fBitmasks[i].fAllBits)
1017 + slice * fNPerSlice;
1018 }
1019 unsigned int* maskCrossed = crossed != nullptr
1020 ? (unsigned int*)crossed->fAllBits : nullptr;
1021
1022 for (G4int i = 0 ; i < fNPerSlice; ++i)
1023 {
1024 // Logic "and" of the masks along the 3 axes x, y, z:
1025 // removing "if (!" and ") continue" => slightly slower
1026 //
1027 if ((mask = masks[0][i]) == 0u) continue;
1028 if ((mask &= masks[1][i]) == 0u) continue;
1029 if ((mask &= masks[2][i]) == 0u) continue;
1030 if ((maskCrossed != nullptr) && ((mask &= ~maskCrossed[i]) == 0u)) continue;
1031
1032 FindComponentsFastest(mask, list, i);
1033 }
1034 }
1035 /*
1036 if (fNPerSlice == 1)
1037 {
1038 unsigned int mask;
1039 G4int slice = BinarySearch(fBoundaries[0], point.x());
1040 if (!(mask = ((unsigned int *) fBitmasks[0].fAllBits)[slice]
1041 )) return 0;
1042 slice = BinarySearch(fBoundaries[1], point.y());
1043 if (!(mask &= ((unsigned int *) fBitmasks[1].fAllBits)[slice]
1044 )) return 0;
1045 slice = BinarySearch(fBoundaries[2], point.z());
1046 if (!(mask &= ((unsigned int *) fBitmasks[2].fAllBits)[slice]
1047 )) return 0;
1048 if (crossed && (!(mask &= ~((unsigned int *)crossed->fAllBits)[0])))
1049 return 0;
1050
1051 FindComponentsFastest(mask, list, 0);
1052 }
1053 else
1054 {
1055 unsigned int *masks[3], mask; // masks for X,Y,Z axis
1056 for (auto i = 0; i <= 2; ++i)
1057 {
1058 G4int slice = BinarySearch(fBoundaries[i], point[i]);
1059 masks[i] = ((unsigned int *) fBitmasks[i].fAllBits) + slice*fNPerSlice;
1060 }
1061 unsigned int *maskCrossed =
1062 crossed ? (unsigned int *)crossed->fAllBits : 0;
1063
1064 for (G4int i = 0 ; i < fNPerSlice; ++i)
1065 {
1066 // Logic "and" of the masks along the 3 axes x, y, z:
1067 // removing "if (!" and ") continue" => slightly slower
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 &= ~maskCrossed[i])) continue;
1073
1074 FindComponentsFastest(mask, list, i);
1075 }
1076 }
1077 */
1078 }
1079 return (G4int)list.size();
1080 }
1081
1082 //______________________________________________________________________________
1083 G4int
1084 G4Voxelizer::GetCandidatesVoxelArray(const std::vector<G4int>& voxels,
1085 const G4SurfBits bitmasks[],
1086 std::vector<G4int>& list,
1087 G4SurfBits* crossed) const
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[0].fAllBits)[voxels[0]]) == 0u)
1102 return 0;
1103 if ((mask &= ((unsigned int *) bitmasks[1].fAllBits)[voxels[1]]) == 0u)
1104 return 0;
1105 if ((mask &= ((unsigned int *) bitmasks[2].fAllBits)[voxels[2]]) == 0u)
1106 return 0;
1107 if ((crossed != nullptr) && ((mask &= ~((unsigned int *)crossed->fAllBits)[0]) == 0u))
1108 return 0;
1109
1110 FindComponentsFastest(mask, list, 0);
1111 }
1112 else
1113 {
1114 unsigned int *masks[3], mask; // masks for X,Y,Z axis
1115 for (auto i = 0; i <= 2; ++i)
1116 {
1117 masks[i] = ((unsigned int *) bitmasks[i].fAllBits)
1118 + voxels[i]*fNPerSlice;
1119 }
1120 unsigned int *maskCrossed = crossed != nullptr
1121 ? (unsigned int *)crossed->fAllBits : nullptr;
1122
1123 for (G4int i = 0 ; i < fNPerSlice; ++i)
1124 {
1125 // Logic "and" of the masks along the 3 axes x, y, z:
1126 // removing "if (!" and ") continue" => slightly slower
1127 //
1128 if ((mask = masks[0][i]) == 0u) continue;
1129 if ((mask &= masks[1][i]) == 0u) continue;
1130 if ((mask &= masks[2][i]) == 0u) continue;
1131 if ((maskCrossed != nullptr) && ((mask &= ~maskCrossed[i]) == 0u)) continue;
1132
1133 FindComponentsFastest(mask, list, i);
1134 }
1135 }
1136 }
1137 return (G4int)list.size();
1138 }
1139
1140 //______________________________________________________________________________
1141 G4int
1142 G4Voxelizer::GetCandidatesVoxelArray(const std::vector<G4int>& voxels,
1143 std::vector<G4int>& list, G4SurfBits* crossed) const
1144 {
1145 // Method returning the candidates corresponding to the passed point
1146
1147 return GetCandidatesVoxelArray(voxels, fBitmasks, list, crossed);
1148 }
1149
1150 //______________________________________________________________________________
1151 G4bool G4Voxelizer::Contains(const G4ThreeVector& point) const
1152 {
1153 for (auto i = 0; i < 3; ++i)
1154 {
1155 if (point[i] < fBoundaries[i].front() || point[i] > fBoundaries[i].back())
1156 return false;
1157 }
1158 return true;
1159 }
1160
1161 //______________________________________________________________________________
1162 G4double
1163 G4Voxelizer::DistanceToFirst(const G4ThreeVector& point,
1164 const G4ThreeVector& direction) const
1165 {
1166 G4ThreeVector pointShifted = point - fBoundingBoxCenter;
1167 G4double shift = fBoundingBox.DistanceToIn(pointShifted, direction);
1168 return shift;
1169 }
1170
1171 //______________________________________________________________________________
1172 G4double
1173 G4Voxelizer::DistanceToBoundingBox(const G4ThreeVector& point) const
1174 {
1175 G4ThreeVector pointShifted = point - fBoundingBoxCenter;
1176 G4double shift = MinDistanceToBox(pointShifted, fBoundingBoxSize);
1177 return shift;
1178 }
1179
1180 //______________________________________________________________________________
1181 G4double
1182 G4Voxelizer::MinDistanceToBox (const G4ThreeVector& aPoint,
1183 const G4ThreeVector& f)
1184 {
1185 // Estimates the isotropic safety from a point outside the current solid to
1186 // any of its surfaces. The algorithm may be accurate or should provide a
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 inside
1196
1197 G4double safsq = 0.0;
1198 G4int count = 0;
1199 if ( safx > 0 ) { safsq += safx*safx; ++count; }
1200 if ( safy > 0 ) { safsq += safy*safy; ++count; }
1201 if ( safz > 0 ) { safsq += safz*safz; ++count; }
1202 if (count == 1) return safe;
1203 return std::sqrt(safsq);
1204 }
1205
1206 //______________________________________________________________________________
1207 G4double
1208 G4Voxelizer::DistanceToNext(const G4ThreeVector& point,
1209 const G4ThreeVector& direction,
1210 std::vector<G4int>& curVoxel) const
1211 {
1212 G4double shift = kInfinity;
1213
1214 G4int cur = 0; // the smallest index, which would be than increased
1215 for (G4int i = 0; i <= 2; ++i)
1216 {
1217 // Looking for the next voxels on the considered direction X,Y,Z axis
1218 //
1219 const std::vector<G4double>& boundary = fBoundaries[i];
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 corresponding
1243 // to the closest voxel boundary on the ray
1244
1245 if (direction[cur] > 0)
1246 {
1247 if (++curVoxel[cur] >= (G4int) fBoundaries[cur].size() - 1)
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 considered direction X,Y,Z axis
1261 //
1262 const std::vector<G4double> &boundary = fBoundaries[i];
1263 G4int cur = curVoxel[i];
1264 if(direction[i] >= 1e-10)
1265 {
1266 if (boundary[++cur] - point[i] < fTolerance) // make sure shift would
1267 if (++cur >= (G4int) boundary.size()) // be non-zero
1268 continue;
1269 }
1270 else
1271 {
1272 if(direction[i] <= -1e-10)
1273 {
1274 if (point[i] - boundary[cur] < fTolerance) // make sure shift would
1275 if (--cur < 0) // be non-zero
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 G4ThreeVector& point,
1294 const G4ThreeVector& direction,
1295 std::vector<G4int>& curVoxel) const
1296 {
1297 for (auto i = 0; i <= 2; ++i)
1298 {
1299 G4int index = curVoxel[i];
1300 const std::vector<G4double> &boundary = fBoundaries[i];
1301
1302 if (direction[i] > 0)
1303 {
1304 if (point[i] >= boundary[++index])
1305 if (++curVoxel[i] >= (G4int) boundary.size() - 1)
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, point[i]);
1316 if (curVoxel[i] != indexOK)
1317 curVoxel[i] = indexOK; // put breakpoint here
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 G4ThreeVector& ratioOfReduction)
1332 {
1333 fMaxVoxels = -1;
1334 fReductionRatio = ratioOfReduction;
1335 }
1336
1337 //______________________________________________________________________________
1338 void G4Voxelizer::SetDefaultVoxelsCount(G4int count)
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(G4VoxelBox);
1354 size += sizeof(G4double) * (fBoundaries[0].capacity()
1355 + fBoundaries[1].capacity() + fBoundaries[2].capacity());
1356 size += sizeof(G4int) * (fCandidatesCounts[0].capacity()
1357 + fCandidatesCounts[1].capacity() + fCandidatesCounts[2].capacity());
1358 size += fBitmasks[0].GetNbytes() + fBitmasks[1].GetNbytes()
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>) + fCandidates[i].capacity() * sizeof(G4int);
1365 }
1366
1367 return size;
1368 }
1369