Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4SurfaceVoxelizer.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4SurfaceVoxelizer.cc 106567 2017-10-13 09:45:14Z gcosmo $
27 //
28 // --------------------------------------------------------------------
29 // GEANT 4 class header file
30 //
31 // G4SurfaceVoxelizer implementation
32 //
33 // History:
34 // 19.10.12 Marek Gayer, created
35 // --------------------------------------------------------------------
36 
37 #include <iostream>
38 #include <iomanip>
39 #include <sstream>
40 #include <algorithm>
41 #include <set>
42 
43 #include "G4VSolid.hh"
44 
45 #include "G4Orb.hh"
46 #include "G4SurfaceVoxelizer.hh"
47 #include "G4SolidStore.hh"
48 #include "Randomize.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4GeometryTolerance.hh"
51 #include "G4CSGSolid.hh"
52 #include "G4Types.hh"
53 #include "geomdefs.hh"
54 
55 using namespace std;
56 
57 G4ThreadLocal G4int G4SurfaceVoxelizer::fDefaultVoxelsCount = -1;
58 
59 //______________________________________________________________________________
61  : fBoundingBox("TessBBox", 1, 1, 1)
62 {
63  fCountOfVoxels = fNPerSlice = fTotalCandidates = 0;
64 
66 
67  SetMaxVoxels(fDefaultVoxelsCount);
68 
69  G4SolidStore::GetInstance()->DeRegister(&fBoundingBox);
70 }
71 
72 //______________________________________________________________________________
74 {
75 }
76 
77 //______________________________________________________________________________
78 void G4SurfaceVoxelizer::BuildEmpty()
79 {
80  // by reserving the size of candidates, we would avoid reallocation of
81  // the vector which could cause fragmentation
82  //
83  vector<G4int> xyz(3), max(3), candidates(fTotalCandidates);
84  const vector<G4int> empty(0);
85 
86  for (G4int i = 0; i <= 2; i++) max[i] = fBoundaries[i].size();
87  unsigned int size = max[0] * max[1] * max[2];
88 
89  fEmpty.Clear();
90  fEmpty.ResetBitNumber(size-1);
91  fEmpty.ResetAllBits(true);
92 
93  for (xyz[2] = 0; xyz[2] < max[2]; ++xyz[2])
94  {
95  for (xyz[1] = 0; xyz[1] < max[1]; ++xyz[1])
96  {
97  for (xyz[0] = 0; xyz[0] < max[0]; ++xyz[0])
98  {
99  G4int index = GetVoxelsIndex(xyz);
100  if (GetCandidatesVoxelArray(xyz, candidates))
101  {
102  fEmpty.SetBitNumber(index, false);
103 
104  // rather than assigning directly with:
105  // "fCandidates[index] = candidates;", in an effort to ensure that
106  // capacity would be just exact, we rather use following 3 lines
107  //
108  vector<G4int> &c = (fCandidates[index] = empty);
109  c.reserve(candidates.size());
110  c.assign(candidates.begin(), candidates.end());
111  }
112  }
113  }
114  }
115 #ifdef G4SPECSDEBUG
116  G4cout << "Non-empty voxels count: " << fCandidates.size() << G4endl;
117 #endif
118 }
119 
120 //______________________________________________________________________________
121 void G4SurfaceVoxelizer::BuildVoxelLimits(std::vector<G4VFacet *> &facets)
122 {
123  // "BuildVoxelLimits"'s aim is to store the coordinates of the origin as well
124  // as the half lengths related to the bounding box of each node.
125  // These quantities are stored in the array "fBoxes" (6 different values per
126  // node.
127 
128  if (G4int numNodes = facets.size()) // Number of nodes
129  {
130  fBoxes.resize(numNodes); // Array which will store the half lengths
131  fNPerSlice = 1+(fBoxes.size()-1)/(8*sizeof(unsigned int));
132 
133  G4ThreeVector toleranceVector;
134  toleranceVector.set(10 * fTolerance,10 * fTolerance,10 * fTolerance);
135 
136  for (G4int i = 0; i < numNodes; ++i)
137  {
138  G4VFacet &facet = *facets[i];
140  G4ThreeVector x(1,0,0), y(0,1,0), z(0,0,1);
141  G4ThreeVector extent;
142  max.set (facet.Extent(x), facet.Extent(y), facet.Extent(z));
143  min.set (-facet.Extent(-x), -facet.Extent(-y), -facet.Extent(-z));
144  min -= toleranceVector; max += toleranceVector;
145  G4ThreeVector hlen = (max - min) / 2;
146  fBoxes[i].hlen = hlen;
147  fBoxes[i].pos = min + hlen;
148  }
149  fTotalCandidates = fBoxes.size();
150  }
151 }
152 
153 //______________________________________________________________________________
155 {
156  // "DisplayVoxelLimits" displays the dX, dY, dZ, pX, pY and pZ for each node
157 
158  G4int numNodes = fBoxes.size();
159  G4int oldprec = G4cout.precision(16);
160  for(G4int i = 0; i < numNodes; ++i)
161  {
162  G4cout << setw(10) << setiosflags(ios::fixed) <<
163  " -> Node " << i+1 << ":\n" <<
164  "\t * [x,y,z] = " << fBoxes[i].hlen <<
165  "\t * [x,y,z] = " << fBoxes[i].pos << "\n";
166  }
167  G4cout.precision(oldprec);
168 }
169 
170 //______________________________________________________________________________
171 void G4SurfaceVoxelizer::CreateSortedBoundary(std::vector<G4double> &boundary,
172  G4int axis)
173 {
174  // "CreateBoundaries"'s aim is to determine the slices induced by the
175  // bounding fBoxes, along each axis. The created boundaries are stored
176  // in the array "boundariesRaw"
177 
178  G4int numNodes = fBoxes.size(); // Number of nodes in structure
179 
180  // Determination of the boundaries along x, y and z axis
181  //
182  for(G4int i = 0 ; i < numNodes; ++i)
183  {
184  // For each node, the boundaries are created by using the array "fBoxes"
185  // built in method "BuildVoxelLimits"
186  //
187  G4double p = fBoxes[i].pos[axis], d = fBoxes[i].hlen[axis];
188 
189  // x boundaries
190  //
191 #ifdef G4SPECSDEBUG
192  G4cout << "Boundary " << p - d << " - " << p + d << G4endl;
193 #endif
194  boundary[2*i] = p - d;
195  boundary[2*i+1] = p + d;
196  }
197  sort(boundary.begin(), boundary.end());
198 }
199 
200 //______________________________________________________________________________
201 void G4SurfaceVoxelizer::BuildBoundaries()
202 {
203  // "SortBoundaries" orders the boundaries along each axis (increasing order)
204  // and also does not take into account redundant boundaries, i.e. if two
205  // boundaries are separated by a distance strictly inferior to "tolerance".
206  // The sorted boundaries are respectively stored in:
207  // * boundaries[0..2]
208 
209  // In addition, the number of elements contained in the three latter arrays
210  // are precise thanks to variables: boundariesCountX, boundariesCountY and
211  // boundariesCountZ.
212 
213  if (G4int numNodes = fBoxes.size())
214  {
215  const G4double tolerance = fTolerance / 100.0;
216  // Minimal distance to discriminate two boundaries.
217 
218  vector<G4double> sortedBoundary(2*numNodes);
219 
220  G4int considered;
221 
222  for (G4int j = 0; j <= 2; ++j)
223  {
224  CreateSortedBoundary(sortedBoundary, j);
225  vector<G4double> &boundary = fBoundaries[j];
226  boundary.clear();
227 
228  considered = 0;
229 
230  for(G4int i = 0 ; i < 2*numNodes; ++i)
231  {
232  G4double newBoundary = sortedBoundary[i];
233 #ifdef G4SPECSDEBUG
234  if (j == 0) G4cout << "Examining " << newBoundary << "..." << G4endl;
235 #endif
236  G4int size = boundary.size();
237  if(!size || abs(boundary[size-1] - newBoundary) > tolerance)
238  {
239  considered++;
240  {
241 #ifdef G4SPECSDEBUG
242  if (j == 0) G4cout << "Adding boundary " << newBoundary << "..."
243  << G4endl;
244 #endif
245  boundary.push_back(newBoundary);
246  continue;
247  }
248  }
249  // If two successive boundaries are too close from each other,
250  // only the first one is considered
251  }
252 
253  G4int n = boundary.size();
254  G4int max = 100000;
255  if (n > max/2)
256  {
257  G4int skip = n / (max /2); // n has to be 2x bigger then 50.000.
258  // therefore only from 100.000 reduced
259  vector<G4double> reduced;
260  for (G4int i = 0; i < n; i++)
261  {
262  // 50 ok for 2k, 1000, 2000
263  G4int size = boundary.size();
264  if (i % skip == 0 || i == 0 || i == size - 1)
265  {
266  // this condition of merging boundaries was wrong,
267  // it did not count with right part, which can be
268  // completely ommited and not included in final consideration.
269  // Now should be OK
270  //
271  reduced.push_back(boundary[i]);
272  }
273  }
274  boundary = reduced;
275  }
276  }
277  }
278 }
279 
280 //______________________________________________________________________________
282 {
283  char axis[3] = {'X', 'Y', 'Z'};
284  for (G4int i = 0; i <= 2; ++i)
285  {
286  G4cout << " * " << axis[i] << " axis:" << G4endl << " | ";
287  DisplayBoundaries(fBoundaries[i]);
288  }
289 }
290 
291 //______________________________________________________________________________
292 void G4SurfaceVoxelizer::DisplayBoundaries(std::vector<G4double> &boundaries)
293 {
294  // Prints the positions of the boundaries of the slices on the three axes
295 
296  G4int count = boundaries.size();
297  G4int oldprec = G4cout.precision(16);
298  for(G4int i = 0; i < count; ++i)
299  {
300  G4cout << setw(10) << setiosflags(ios::fixed) << boundaries[i];
301  if(i != count-1) G4cout << "-> ";
302  }
303  G4cout << "|" << G4endl << "Number of boundaries: " << count << G4endl;
304  G4cout.precision(oldprec);
305 }
306 
307 //______________________________________________________________________________
308 void G4SurfaceVoxelizer::BuildBitmasks(std::vector<G4double> boundaries[],
309  G4SurfBits bitmasks[], G4bool countsOnly)
310 {
311  // "BuildListNodes" stores in the bitmasks solids present in each slice
312  // along an axis.
313 
314  G4int numNodes = fBoxes.size();
315  G4int bitsPerSlice = GetBitsPerSlice();
316 
317  for (G4int k = 0; k < 3; k++)
318  {
319  G4int total = 0;
320  vector<G4double> &boundary = boundaries[k];
321  G4int voxelsCount = boundary.size() - 1;
322  G4SurfBits &bitmask = bitmasks[k];
323 
324  if (!countsOnly)
325  {
326  bitmask.Clear();
327 #ifdef G4SPECSDEBUG
328  G4cout << "Allocating bitmask..." << G4endl;
329 #endif
330  bitmask.SetBitNumber(voxelsCount*bitsPerSlice-1, false);
331  // it is here so we can set the maximum number of bits. this line
332  // will rellocate the memory and set all to zero
333  }
334  vector<G4int> &candidatesCount = fCandidatesCounts[k];
335  candidatesCount.resize(voxelsCount);
336 
337  for(G4int i = 0 ; i < voxelsCount; ++i) { candidatesCount[i] = 0; }
338 
339  // Loop on the nodes, number of slices per axis
340  //
341  for(G4int j = 0 ; j < numNodes; j++)
342  {
343  // Determination of the minimum and maximum position along x
344  // of the bounding boxe of each node
345  //
346  G4double p = fBoxes[j].pos[k], d = fBoxes[j].hlen[k];
347 
348  G4double min = p - d; // - localTolerance;
349  G4double max = p + d; // + localTolerance;
350 
351  G4int i = BinarySearch(boundary, min);
352  if (i < 0) { i = 0; }
353 
354  do // Loop checking, 13.08.2015, G.Cosmo
355  {
356  if (!countsOnly)
357  {
358  G4int index = i*bitsPerSlice+j;
359  bitmask.SetBitNumber(index);
360  }
361  candidatesCount[i]++;
362  total++;
363  i++;
364  }
365  while (max > boundary[i] && i < voxelsCount);
366  }
367  }
368 #ifdef G4SPECSDEBUG
369  G4cout << "Build list nodes completed." << G4endl;
370 #endif
371 }
372 
373 //______________________________________________________________________________
374 G4String G4SurfaceVoxelizer::GetCandidatesAsString(const G4SurfBits &bits) const
375 {
376  // Decodes the candidates in mask as G4String.
377 
378  stringstream ss;
379  G4int numNodes = fBoxes.size();
380 
381  for(G4int i=0; i<numNodes; ++i)
382  {
383  if (bits.TestBitNumber(i)) { ss << i+1 << " "; }
384  }
385  return ss.str();
386 }
387 
388 //______________________________________________________________________________
390 {
391  // Prints which solids are present in the slices previously elaborated.
392 
393  char axis[3] = {'X', 'Y', 'Z'};
394  G4int size=8*sizeof(G4int)*fNPerSlice;
395  G4SurfBits bits(size);
396 
397  for (G4int j = 0; j <= 2; ++j)
398  {
399  G4cout << " * " << axis[j] << " axis:" << G4endl;
400  G4int count = fBoundaries[j].size();
401  for(G4int i=0; i < count-1; ++i)
402  {
403  G4cout << " Slice #" << i+1 << ": [" << fBoundaries[j][i]
404  << " ; " << fBoundaries[j][i+1] << "] -> ";
405  bits.set(size,(const char *)fBitmasks[j].fAllBits+i
406  *fNPerSlice*sizeof(G4int));
407  G4String result = GetCandidatesAsString(bits);
408  G4cout << "[ " << result.c_str() << "] " << G4endl;
409  }
410  }
411 }
412 
413 //______________________________________________________________________________
414 void G4SurfaceVoxelizer::BuildBoundingBox()
415 {
416  G4ThreeVector toleranceVector;
417  toleranceVector.set(fTolerance/100,fTolerance/100,fTolerance/100);
418  for (G4int i = 0; i <= 2; ++i)
419  {
420  G4double min = fBoundaries[i].front();
421  G4double max = fBoundaries[i].back();
422  fBoundingBoxSize[i] = (max-min)/2;
423  fBoundingBoxCenter[i] = min + fBoundingBoxSize[i];
424  }
425  fBoundingBox.SetXHalfLength(fBoundingBoxSize.x());
426  fBoundingBox.SetYHalfLength(fBoundingBoxSize.y());
427  fBoundingBox.SetZHalfLength(fBoundingBoxSize.z());
428 }
429 
430 // algorithm -
431 
432 // in order to get balanced voxels, merge should always unite those regions,
433 // where the number of voxels is least the number.
434 // We will keep sorted list (std::set) with all voxels. there will be
435 // comparator function between two voxels, which will tell if voxel is less
436 // by looking at his right neighbor.
437 // First, we will add all the voxels into the tree.
438 // We will be pick the first item in the tree, merging it, adding the right
439 // merged voxel into the a list for future reduction (fBitmasks will be
440 // rebuilded later, therefore they need not to be updated).
441 // The merged voxel need to be added to the tree again, so it's position
442 // would be updated.
443 
444 //______________________________________________________________________________
445 void G4SurfaceVoxelizer::SetReductionRatio(G4int maxVoxels,
446  G4ThreeVector &reductionRatio)
447 {
448  G4double maxTotal = (G4double) fCandidatesCounts[0].size()
449  * fCandidatesCounts[1].size() * fCandidatesCounts[2].size();
450 
451  if (maxVoxels > 0 && maxVoxels < maxTotal)
452  {
453  G4double ratio = (G4double) maxVoxels / maxTotal;
454  ratio = pow(ratio, 1./3.);
455  if (ratio > 1) { ratio = 1; }
456  reductionRatio.set(ratio,ratio,ratio);
457  }
458 }
459 
460 //______________________________________________________________________________
461 void G4SurfaceVoxelizer::BuildReduceVoxels(std::vector<G4double> boundaries[],
462  G4ThreeVector reductionRatio)
463 {
464  for (G4int k = 0; k <= 2; ++k)
465  {
466  vector<G4int> &candidatesCount = fCandidatesCounts[k];
467  G4int max = candidatesCount.size();
468  vector<G4VoxelInfo> voxels(max);
469  G4VoxelComparator comp(voxels);
470  set<G4int, G4VoxelComparator> voxelSet(comp);
471  vector<G4int> mergings;
472 
473  for (G4int j = 0; j < max; ++j)
474  {
475  G4VoxelInfo &voxel = voxels[j];
476  voxel.count = candidatesCount[j];
477  voxel.previous = j - 1;
478  voxel.next = j + 1;
479  voxels[j] = voxel;
480  }
481 
482  for (G4int j = 0; j < max - 1; ++j) { voxelSet.insert(j); }
483  // we go to size-1 to make sure we will not merge the last element
484 
485  G4int count = 0;
486  while (true) // Loop checking, 13.08.2015, G.Cosmo
487  {
488  G4double reduction = reductionRatio[k];
489  if (reduction == 0)
490  break;
491  G4int currentCount = voxelSet.size();
492  if (currentCount <= 2)
493  break;
494  G4double currentRatio = 1 - (G4double) count / max;
495  if (currentRatio <= reduction && currentCount <= 1000)
496  break;
497  const G4int pos = *voxelSet.begin();
498  mergings.push_back(pos);
499 
500  G4VoxelInfo &voxel = voxels[pos];
501  G4VoxelInfo &nextVoxel = voxels[voxel.next];
502 
503  voxelSet.erase(pos);
504  if (voxel.next != max - 1) { voxelSet.erase(voxel.next); }
505  if (voxel.previous != -1) { voxelSet.erase(voxel.previous); }
506 
507  nextVoxel.count += voxel.count;
508  voxel.count = 0;
509  nextVoxel.previous = voxel.previous;
510 
511  if (voxel.next != max - 1)
512  voxelSet.insert(voxel.next);
513 
514  if (voxel.previous != -1)
515  {
516  voxels[voxel.previous].next = voxel.next;
517  voxelSet.insert(voxel.previous);
518  }
519  count++;
520  }
521 
522  if (mergings.size())
523  {
524  sort(mergings.begin(), mergings.end());
525 
526  vector<G4double> &boundary = boundaries[k];
527  vector<G4double> reducedBoundary(boundary.size() - mergings.size());
528  G4int skip = mergings[0] + 1, cur = 0, i = 0;
529  max = boundary.size();
530  for (G4int j = 0; j < max; ++j)
531  {
532  if (j != skip)
533  {
534  reducedBoundary[cur++] = boundary[j];
535  }
536  else
537  {
538  if (++i < (G4int)mergings.size()) { skip = mergings[i] + 1; }
539  }
540  }
541  boundaries[k] = reducedBoundary;
542  }
543  }
544 }
545 
546 //______________________________________________________________________________
547 void G4SurfaceVoxelizer::BuildReduceVoxels2(std::vector<G4double> boundaries[],
548  G4ThreeVector reductionRatio)
549 {
550  for (G4int k = 0; k <= 2; ++k)
551  {
552  vector<G4int> &candidatesCount = fCandidatesCounts[k];
553  G4int max = candidatesCount.size();
554  G4int total = 0;
555  for (G4int i = 0; i < max; i++) total += candidatesCount[i];
556 
557  G4double reduction = reductionRatio[k];
558  if (reduction == 0)
559  break;
560 
561  G4int destination = (G4int) (reduction * max) + 1;
562  if (destination > 1000) destination = 1000;
563  if (destination < 2) destination = 2;
564  G4double average = ((G4double)total / max) / reduction;
565 
566  vector<G4int> mergings;
567 
568  vector<G4double> &boundary = boundaries[k];
569  vector<G4double> reducedBoundary(destination);
570 
571  G4int sum = 0, cur = 0;
572  for (G4int i = 0; i < max; i++)
573  {
574  sum += candidatesCount[i];
575  if (sum > average * (cur + 1) || i == 0)
576  {
577  G4double val = boundary[i];
578  reducedBoundary[cur] = val;
579  cur++;
580  if (cur == destination)
581  break;
582  }
583  }
584  reducedBoundary[destination-1] = boundary[max];
585  boundaries[k] = reducedBoundary;
586  }
587 }
588 
589 //______________________________________________________________________________
590 void G4SurfaceVoxelizer::CreateMiniVoxels(std::vector<G4double> boundaries[],
591  G4SurfBits bitmasks[])
592 {
593  vector<G4int> voxel(3), maxVoxels(3);
594  for (G4int i = 0; i <= 2; ++i) maxVoxels[i] = boundaries[i].size();
595 
596  G4ThreeVector point;
597  for (voxel[2] = 0; voxel[2] < maxVoxels[2] - 1; ++voxel[2])
598  {
599  for (voxel[1] = 0; voxel[1] < maxVoxels[1] - 1; ++voxel[1])
600  {
601  for (voxel[0] = 0; voxel[0] < maxVoxels[0] - 1; ++voxel[0])
602  {
603  vector<G4int> candidates;
604  if (GetCandidatesVoxelArray(voxel, bitmasks, candidates, 0))
605  {
606  // find a box for corresponding non-empty voxel
607  G4VoxelBox box;
608  for (G4int i = 0; i <= 2; ++i)
609  {
610  G4int index = voxel[i];
611  const vector<G4double> &boundary = boundaries[i];
612  G4double hlen = 0.5 * (boundary[index+1] - boundary[index]);
613  box.hlen[i] = hlen;
614  box.pos[i] = boundary[index] + hlen;
615  }
616  fVoxelBoxes.push_back(box);
617  vector<G4int>(candidates).swap(candidates);
618  fVoxelBoxesCandidates.push_back(candidates);
619  }
620  }
621  }
622  }
623 }
624 
625 //______________________________________________________________________________
626 void G4SurfaceVoxelizer::Voxelize(std::vector<G4VFacet *> &facets)
627 {
628  G4int maxVoxels = fMaxVoxels;
629  G4ThreeVector reductionRatio = fReductionRatio;
630 
631  G4int size = facets.size();
632  if (size < 10)
633  {
634  for (G4int i = 0; i < (G4int) facets.size(); i++)
635  {
636  if (facets[i]->GetNumberOfVertices() > 3) size++;
637  }
638  }
639 
640  if ((size >= 10 || maxVoxels > 0) && maxVoxels != 0 && maxVoxels != 1)
641  {
642 #ifdef G4SPECSDEBUG
643  G4cout << "Building voxel limits..." << G4endl;
644 #endif
645 
646  BuildVoxelLimits(facets);
647 
648 #ifdef G4SPECSDEBUG
649  G4cout << "Building boundaries..." << G4endl;
650 #endif
651 
652  BuildBoundaries();
653 
654 #ifdef G4SPECSDEBUG
655  G4cout << "Building bitmasks..." << G4endl;
656 #endif
657 
658  BuildBitmasks(fBoundaries, 0, true);
659 
660  if (maxVoxels < 0 && reductionRatio == G4ThreeVector())
661  {
662  maxVoxels = fTotalCandidates;
663  if (fTotalCandidates > 1000000) maxVoxels = 1000000;
664  }
665 
666  SetReductionRatio(maxVoxels, reductionRatio);
667 
668  fCountOfVoxels = CountVoxels(fBoundaries);
669 
670 #ifdef G4SPECSDEBUG
671  G4cout << "Total number of voxels: " << fCountOfVoxels << G4endl;
672 #endif
673 
674  BuildReduceVoxels2(fBoundaries, reductionRatio);
675 
676  fCountOfVoxels = CountVoxels(fBoundaries);
677 
678 #ifdef G4SPECSDEBUG
679  G4cout << "Total number of voxels after reduction: "
680  << fCountOfVoxels << G4endl;
681 #endif
682 
683 #ifdef G4SPECSDEBUG
684  G4cout << "Building bitmasks..." << G4endl;
685 #endif
686 
687  BuildBitmasks(fBoundaries, fBitmasks);
688 
689  G4ThreeVector reductionRatioMini;
690 
691  G4SurfBits bitmasksMini[3];
692 
693  // section for building mini voxels
694 
695  vector<G4double> miniBoundaries[3];
696 
697  for (G4int i = 0; i <= 2; ++i) { miniBoundaries[i] = fBoundaries[i]; }
698 
699  G4int voxelsCountMini = (fCountOfVoxels >= 1000)
700  ? 100 : fCountOfVoxels / 10;
701 
702  // if (voxelCountMini < 8) voxelCountMini = 8;
703  // voxelsCountMini = 1;
704 
705  SetReductionRatio(voxelsCountMini, reductionRatioMini);
706 
707 #ifdef G4SPECSDEBUG
708  G4cout << "Building reduced voxels..." << G4endl;
709 #endif
710 
711  BuildReduceVoxels(miniBoundaries, reductionRatioMini);
712 
713 #ifdef G4SPECSDEBUG
714  G4int total = CountVoxels(miniBoundaries);
715  G4cout << "Total number of mini voxels: " << total << G4endl;
716 #endif
717 
718 #ifdef G4SPECSDEBUG
719  G4cout << "Building mini bitmasks..." << G4endl;
720 #endif
721 
722  BuildBitmasks(miniBoundaries, bitmasksMini);
723 
724 #ifdef G4SPECSDEBUG
725  G4cout << "Creating Mini Voxels..." << G4endl;
726 #endif
727 
728  CreateMiniVoxels(miniBoundaries, bitmasksMini);
729 
730 #ifdef G4SPECSDEBUG
731  G4cout << "Building bounding box..." << G4endl;
732 #endif
733 
734  BuildBoundingBox();
735 
736 #ifdef G4SPECSDEBUG
737  G4cout << "Building empty..." << G4endl;
738 #endif
739 
740  BuildEmpty();
741 
742 #ifdef G4SPECSDEBUG
743  G4cout << "Deallocating unnecessary fields during runtime..." << G4endl;
744 #endif
745  // deallocate fields unnecessary during runtime
746  //
747  fBoxes.resize(0);
748  for (G4int i = 0; i < 3; i++)
749  {
750  fCandidatesCounts[i].resize(0);
751  fBitmasks[i].Clear();
752  }
753  }
754 }
755 
756 
757 //______________________________________________________________________________
758 void G4SurfaceVoxelizer::GetCandidatesVoxel(std::vector<G4int> &voxels)
759 {
760  // "GetCandidates" should compute which solids are possibly contained in
761  // the voxel defined by the three slices characterized by the passed indexes.
762 
763  G4cout << " Candidates in voxel [" << voxels[0] << " ; " << voxels[1]
764  << " ; " << voxels[2] << "]: ";
765  vector<G4int> candidates;
766  G4int count = GetCandidatesVoxelArray(voxels, candidates);
767  G4cout << "[ ";
768  for (G4int i = 0; i < count; ++i) G4cout << candidates[i];
769  G4cout << "] " << G4endl;
770 }
771 
772 //______________________________________________________________________________
773 void G4SurfaceVoxelizer::FindComponentsFastest(unsigned int mask,
774  std::vector<G4int> &list, G4int i)
775 {
776  for (G4int byte = 0; byte < (G4int) (sizeof(unsigned int)); byte++)
777  {
778  if (G4int maskByte = mask & 0xFF)
779  {
780  for (G4int bit = 0; bit < 8; bit++)
781  {
782  if (maskByte & 1)
783  { list.push_back(8*(sizeof(unsigned int)*i+ byte) + bit); }
784  if (!(maskByte >>= 1)) break;
785  }
786  }
787  mask >>= 8;
788  }
789 }
790 
791 //______________________________________________________________________________
793  std::vector<G4int> &list, G4SurfBits *crossed) const
794 {
795  // Method returning the candidates corresponding to the passed point
796 
797  list.clear();
798 
799  for (G4int i = 0; i <= 2; ++i)
800  {
801  if(point[i] < fBoundaries[i].front() || point[i] >= fBoundaries[i].back())
802  return 0;
803  }
804 
805  if (fTotalCandidates == 1)
806  {
807  list.push_back(0);
808  return 1;
809  }
810  else
811  {
812  if (fNPerSlice == 1)
813  {
814  unsigned int mask;
815  G4int slice = BinarySearch(fBoundaries[0], point.x());
816  if (!(mask = ((unsigned int *) fBitmasks[0].fAllBits)[slice]
817  )) return 0;
818  slice = BinarySearch(fBoundaries[1], point.y());
819  if (!(mask &= ((unsigned int *) fBitmasks[1].fAllBits)[slice]
820  )) return 0;
821  slice = BinarySearch(fBoundaries[2], point.z());
822  if (!(mask &= ((unsigned int *) fBitmasks[2].fAllBits)[slice]
823  )) return 0;
824  if (crossed && (!(mask &= ~((unsigned int *)crossed->fAllBits)[0])))
825  return 0;
826 
827  FindComponentsFastest(mask, list, 0);
828  }
829  else
830  {
831  unsigned int *masks[3], mask; // masks for X,Y,Z axis
832  for (G4int i = 0; i <= 2; ++i)
833  {
834  G4int slice = BinarySearch(fBoundaries[i], point[i]);
835  masks[i] = ((unsigned int *) fBitmasks[i].fAllBits) + slice*fNPerSlice;
836  }
837  unsigned int *maskCrossed =
838  crossed ? (unsigned int *)crossed->fAllBits : 0;
839 
840  for (G4int i = 0 ; i < fNPerSlice; ++i)
841  {
842  // Logic "and" of the masks along the 3 axes x, y, z:
843  // removing "if (!" and ") continue" => slightly slower
844  //
845  if (!(mask = masks[0][i])) continue;
846  if (!(mask &= masks[1][i])) continue;
847  if (!(mask &= masks[2][i])) continue;
848  if (maskCrossed && !(mask &= ~maskCrossed[i])) continue;
849 
850  FindComponentsFastest(mask, list, i);
851  }
852  }
853  }
854  return list.size();
855 }
856 
857 //______________________________________________________________________________
858 G4int
859 G4SurfaceVoxelizer::GetCandidatesVoxelArray(const std::vector<G4int> &voxels,
860  const G4SurfBits bitmasks[],
861  std::vector<G4int> &list,
862  G4SurfBits *crossed) const
863 {
864  list.clear();
865 
866  if (fTotalCandidates == 1)
867  {
868  list.push_back(0);
869  return 1;
870  }
871  else
872  {
873  if (fNPerSlice == 1)
874  {
875  unsigned int mask;
876  if (!(mask = ((unsigned int *) bitmasks[0].fAllBits)[voxels[0]]
877  )) return 0;
878  if (!(mask &= ((unsigned int *) bitmasks[1].fAllBits)[voxels[1]]
879  )) return 0;
880  if (!(mask &= ((unsigned int *) bitmasks[2].fAllBits)[voxels[2]]
881  )) return 0;
882  if (crossed && (!(mask &= ~((unsigned int *)crossed->fAllBits)[0])))
883  return 0;
884 
885  FindComponentsFastest(mask, list, 0);
886  }
887  else
888  {
889  unsigned int *masks[3], mask; // masks for X,Y,Z axis
890  for (G4int i = 0; i <= 2; ++i)
891  {
892  masks[i] = ((unsigned int *) bitmasks[i].fAllBits)
893  + voxels[i]*fNPerSlice;
894  }
895  unsigned int *maskCrossed =
896  crossed ? (unsigned int *)crossed->fAllBits : 0;
897 
898  for (G4int i = 0 ; i < fNPerSlice; ++i)
899  {
900  // Logic "and" of the masks along the 3 axes x, y, z:
901  // removing "if (!" and ") continue" => slightly slower
902  //
903  if (!(mask = masks[0][i])) continue;
904  if (!(mask &= masks[1][i])) continue;
905  if (!(mask &= masks[2][i])) continue;
906  if (maskCrossed && !(mask &= ~maskCrossed[i])) continue;
907 
908  FindComponentsFastest(mask, list, i);
909  }
910  }
911  }
912  return list.size();
913 }
914 
915 //______________________________________________________________________________
916 G4int
917 G4SurfaceVoxelizer::GetCandidatesVoxelArray(const std::vector<G4int> &voxels,
918  std::vector<G4int> &list, G4SurfBits *crossed) const
919 {
920  // Method returning the candidates corresponding to the passed point
921 
922  return GetCandidatesVoxelArray(voxels, fBitmasks, list, crossed);
923 }
924 
925 //______________________________________________________________________________
927 {
928  for (G4int i = 0; i < 3; ++i)
929  {
930  if (point[i] < fBoundaries[i].front() || point[i] > fBoundaries[i].back())
931  return false;
932  }
933  return true;
934 }
935 
936 //______________________________________________________________________________
937 G4double
939  const G4ThreeVector &direction) const
940 {
941  G4ThreeVector pointShifted = point - fBoundingBoxCenter;
942  G4double shift = fBoundingBox.DistanceToIn(pointShifted, direction);
943  return shift;
944 }
945 
946 //______________________________________________________________________________
947 G4double
949 {
950  G4ThreeVector pointShifted = point - fBoundingBoxCenter;
951  G4double shift = MinDistanceToBox(pointShifted, fBoundingBoxSize);
952  return shift;
953 }
954 
955 //______________________________________________________________________________
956 G4double
958  const G4ThreeVector &f)
959 {
960  // Estimates the isotropic safety from a point outside the current solid to
961  // any of its surfaces. The algorithm may be accurate or should provide a
962  // fast underestimate.
963 
964  G4double safe, safx, safy, safz;
965  safe = safx = -f.x() + std::abs(aPoint.x());
966  safy = -f.y() + std::abs(aPoint.y());
967  if ( safy > safe ) safe = safy;
968  safz = -f.z() + std::abs(aPoint.z());
969  if ( safz > safe ) safe = safz;
970  if (safe < 0.0) return 0.0; // point is inside
971  // if (!aAccurate) return safe;
972  G4double safsq = 0.0;
973  G4int count = 0;
974  if ( safx > 0 ) { safsq += safx*safx; count++; }
975  if ( safy > 0 ) { safsq += safy*safy; count++; }
976  if ( safz > 0 ) { safsq += safz*safz; count++; }
977  if (count == 1) return safe;
978  return std::sqrt(safsq);
979 }
980 
981 //______________________________________________________________________________
982 G4double
984  const G4ThreeVector &direction,
985  const std::vector<G4int> &curVoxel) const
986 {
987  G4double shift = kInfinity;
988 
989  for (G4int i = 0; i <= 2; ++i)
990  {
991  // Looking for the next voxels on the considered direction X,Y,Z axis
992  //
993  const vector<G4double> &boundary = fBoundaries[i];
994  G4int cur = curVoxel[i];
995  if(direction[i] >= 1e-10)
996  {
997  if (boundary[++cur] - point[i] < fTolerance) // make sure shift would
998  if (++cur >= (G4int) boundary.size()) // be non-zero
999  continue;
1000  }
1001  else
1002  {
1003  if(direction[i] <= -1e-10)
1004  {
1005  if (point[i] - boundary[cur] < fTolerance) // make sure shift would
1006  if (--cur < 0) // be non-zero
1007  continue;
1008  }
1009  else
1010  continue;
1011  }
1012  G4double dif = boundary[cur] - point[i];
1013  G4double distance = dif / direction[i];
1014 
1015  if (shift > distance)
1016  shift = distance;
1017  }
1018 
1019  return shift;
1020 }
1021 
1022 //______________________________________________________________________________
1023 G4bool
1025  const G4ThreeVector &direction,
1026  std::vector<G4int> &curVoxel) const
1027 {
1028  for (G4int i = 0; i <= 2; ++i)
1029  {
1030  G4int index = curVoxel[i];
1031  const vector<G4double> &boundary = fBoundaries[i];
1032 
1033  if (direction[i] > 0)
1034  {
1035  if (point[i] >= boundary[++index])
1036  if (++curVoxel[i] >= (G4int) boundary.size() - 1)
1037  return false;
1038  }
1039  else
1040  {
1041  if (point[i] < boundary[index])
1042  if (--curVoxel[i] < 0)
1043  return false;
1044  }
1045 #ifdef G4SPECSDEBUG
1046  G4int indexOK = BinarySearch(boundary, point[i]);
1047  if (curVoxel[i] != indexOK)
1048  curVoxel[i] = indexOK; // put breakpoint here
1049 #endif
1050  }
1051  return true;
1052 }
1053 
1054 //______________________________________________________________________________
1056 {
1057  fMaxVoxels = max;
1058  fReductionRatio.set(0,0,0);
1059 }
1060 
1061 //______________________________________________________________________________
1063 {
1064  fMaxVoxels = -1;
1065  fReductionRatio = ratioOfReduction;
1066 }
1067 
1068 //______________________________________________________________________________
1070 {
1071  G4int res = fDefaultVoxelsCount;
1072  fDefaultVoxelsCount = count;
1073  return res;
1074 }
1075 
1076 //______________________________________________________________________________
1078 {
1079  return fDefaultVoxelsCount;
1080 }
1081 
1082 //______________________________________________________________________________
1084 {
1085  G4int size = fEmpty.GetNbytes();
1086  size += fBoxes.capacity() * sizeof(G4VoxelBox);
1087  size += sizeof(G4double) * (fBoundaries[0].capacity()
1088  + fBoundaries[1].capacity() + fBoundaries[2].capacity());
1089  size += sizeof(G4int) * (fCandidatesCounts[0].capacity()
1090  + fCandidatesCounts[1].capacity() + fCandidatesCounts[2].capacity());
1091  size += fBitmasks[0].GetNbytes() + fBitmasks[1].GetNbytes()
1092  + fBitmasks[2].GetNbytes();
1093 
1094  G4int csize = fCandidates.size();
1095  for (G4int i = 0; i < csize; i++)
1096  {
1097  size += sizeof(vector<G4int>) + fCandidates[i].capacity() * sizeof(G4int);
1098  }
1099 
1100  return size;
1101 }
1102 
G4double G4ParticleHPJENDLHEData::G4double result
void set(double x, double y, double z)
G4bool Contains(const G4ThreeVector &point) const
void ResetBitNumber(unsigned int bitnumber)
Definition: G4SurfBits.hh:161
void Clear()
Definition: G4SurfBits.cc:92
void DisplayListNodes() const
void SetZHalfLength(G4double dz)
Definition: G4Box.cc:174
static const G4double kInfinity
Definition: geomdefs.hh:42
G4int GetCandidatesVoxelArray(const G4ThreeVector &point, std::vector< G4int > &list, G4SurfBits *crossed=0) const
void Voxelize(std::vector< G4VFacet * > &facets)
CLHEP::Hep3Vector G4ThreeVector
static void DeRegister(G4VSolid *pSolid)
double x() const
static ush mask[]
Definition: csz_inflate.cc:317
const char * p
Definition: xmltok.h:285
G4double DistanceToFirst(const G4ThreeVector &point, const G4ThreeVector &direction) const
G4double GetSurfaceTolerance() const
G4double DistanceToBoundingBox(const G4ThreeVector &point) const
G4int GetVoxelsIndex(G4int x, G4int y, G4int z) const
tuple x
Definition: test.py:50
static G4double MinDistanceToBox(const G4ThreeVector &aPoint, const G4ThreeVector &f)
#define G4ThreadLocal
Definition: tls.hh:89
G4bool UpdateCurrentVoxel(const G4ThreeVector &point, const G4ThreeVector &direction, std::vector< G4int > &curVoxel) const
int G4int
Definition: G4Types.hh:78
G4int GetBitsPerSlice() const
double z() const
void SetBitNumber(unsigned int bitnumber, G4bool value=true)
Definition: G4SurfBits.hh:123
unsigned char * fAllBits
Definition: G4SurfBits.hh:113
static G4int BinarySearch(const std::vector< T > &vec, T value)
G4bool TestBitNumber(unsigned int bitnumber) const
Definition: G4SurfBits.hh:148
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
G4double DistanceToNext(const G4ThreeVector &point, const G4ThreeVector &direction, const std::vector< G4int > &curVoxel) const
static G4SolidStore * GetInstance()
void SetMaxVoxels(G4int max)
const G4int n
void ResetAllBits(G4bool value=false)
Definition: G4SurfBits.cc:154
void set(unsigned int nbits, const char *array)
Definition: G4SurfBits.cc:174
static G4int SetDefaultVoxelsCount(G4int count)
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
Definition: G4Box.cc:436
G4double total(Particle const *const p1, Particle const *const p2)
virtual G4double Extent(const G4ThreeVector)=0
G4ThreeVector pos
T max(const T t1, const T t2)
brief Return the largest of the two arguments
double y() const
void SetYHalfLength(G4double dy)
Definition: G4Box.cc:154
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void SetXHalfLength(G4double dx)
Definition: G4Box.cc:134
tuple z
Definition: test.py:28
void GetCandidatesVoxel(std::vector< G4int > &voxels)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
G4ThreeVector hlen
long long CountVoxels(std::vector< G4double > boundaries[]) const
static const G4double pos
for(G4int i1=0;i1< theStableOnes.GetNumberOfIsotopes(static_cast< G4int >(anE->GetZ()));i1++)
static G4GeometryTolerance * GetInstance()
static G4int GetDefaultVoxelsCount()
unsigned int GetNbytes() const
Definition: G4SurfBits.hh:102