61 : fBoundingBox(
"TessBBox", 1, 1, 1)
63 fCountOfVoxels = fNPerSlice = fTotalCandidates = 0;
78 void G4SurfaceVoxelizer::BuildEmpty()
83 vector<G4int> xyz(3),
max(3), candidates(fTotalCandidates);
84 const vector<G4int> empty(0);
86 for (
G4int i = 0; i <= 2; i++)
max[i] = fBoundaries[i].size();
87 unsigned int size =
max[0] *
max[1] *
max[2];
93 for (xyz[2] = 0; xyz[2] < max[2]; ++xyz[2])
95 for (xyz[1] = 0; xyz[1] < max[1]; ++xyz[1])
97 for (xyz[0] = 0; xyz[0] < max[0]; ++xyz[0])
108 vector<G4int> &
c = (fCandidates[
index] = empty);
109 c.reserve(candidates.size());
110 c.assign(candidates.begin(), candidates.end());
116 G4cout <<
"Non-empty voxels count: " << fCandidates.size() <<
G4endl;
121 void G4SurfaceVoxelizer::BuildVoxelLimits(std::vector<G4VFacet *> &facets)
128 if (
G4int numNodes = facets.size())
130 fBoxes.resize(numNodes);
131 fNPerSlice = 1+(fBoxes.size()-1)/(8*
sizeof(
unsigned int));
134 toleranceVector.
set(10 * fTolerance,10 * fTolerance,10 * fTolerance);
136 for (
G4int i = 0; i < numNodes; ++i)
144 min -= toleranceVector; max += toleranceVector;
146 fBoxes[i].hlen = hlen;
147 fBoxes[i].pos = min + hlen;
149 fTotalCandidates = fBoxes.size();
158 G4int numNodes = fBoxes.size();
160 for(
G4int i = 0; i < numNodes; ++i)
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";
167 G4cout.precision(oldprec);
171 void G4SurfaceVoxelizer::CreateSortedBoundary(std::vector<G4double> &boundary,
178 G4int numNodes = fBoxes.size();
182 for(
G4int i = 0 ; i < numNodes; ++i)
187 G4double p = fBoxes[i].pos[axis],
d = fBoxes[i].hlen[axis];
194 boundary[2*i] = p -
d;
195 boundary[2*i+1] = p +
d;
197 sort(boundary.begin(), boundary.end());
201 void G4SurfaceVoxelizer::BuildBoundaries()
213 if (
G4int numNodes = fBoxes.size())
215 const G4double tolerance = fTolerance / 100.0;
218 vector<G4double> sortedBoundary(2*numNodes);
222 for (
G4int j = 0; j <= 2; ++j)
224 CreateSortedBoundary(sortedBoundary, j);
225 vector<G4double> &boundary = fBoundaries[j];
230 for(
G4int i = 0 ; i < 2*numNodes; ++i)
232 G4double newBoundary = sortedBoundary[i];
234 if (j == 0)
G4cout <<
"Examining " << newBoundary <<
"..." <<
G4endl;
236 G4int size = boundary.size();
237 if(!size || abs(boundary[size-1] - newBoundary) > tolerance)
242 if (j == 0)
G4cout <<
"Adding boundary " << newBoundary <<
"..."
245 boundary.push_back(newBoundary);
253 G4int n = boundary.size();
257 G4int skip = n / (max /2);
259 vector<G4double> reduced;
260 for (
G4int i = 0; i <
n; i++)
263 G4int size = boundary.size();
264 if (i % skip == 0 || i == 0 || i == size - 1)
271 reduced.push_back(boundary[i]);
283 char axis[3] = {
'X',
'Y',
'Z'};
284 for (
G4int i = 0; i <= 2; ++i)
286 G4cout <<
" * " << axis[i] <<
" axis:" << G4endl <<
" | ";
296 G4int count = boundaries.size();
298 for(
G4int i = 0; i < count; ++i)
300 G4cout << setw(10) << setiosflags(ios::fixed) << boundaries[i];
301 if(i != count-1)
G4cout <<
"-> ";
303 G4cout <<
"|" << G4endl <<
"Number of boundaries: " << count <<
G4endl;
304 G4cout.precision(oldprec);
308 void G4SurfaceVoxelizer::BuildBitmasks(std::vector<G4double> boundaries[],
314 G4int numNodes = fBoxes.size();
317 for (
G4int k = 0; k < 3; k++)
320 vector<G4double> &boundary = boundaries[k];
321 G4int voxelsCount = boundary.size() - 1;
330 bitmask.
SetBitNumber(voxelsCount*bitsPerSlice-1,
false);
334 vector<G4int> &candidatesCount = fCandidatesCounts[k];
335 candidatesCount.resize(voxelsCount);
337 for(
G4int i = 0 ; i < voxelsCount; ++i) { candidatesCount[i] = 0; }
341 for(
G4int j = 0 ; j < numNodes; j++)
346 G4double p = fBoxes[j].pos[k],
d = fBoxes[j].hlen[k];
352 if (i < 0) { i = 0; }
358 G4int index = i*bitsPerSlice+j;
361 candidatesCount[i]++;
365 while (max > boundary[i] && i < voxelsCount);
379 G4int numNodes = fBoxes.size();
381 for(
G4int i=0; i<numNodes; ++i)
393 char axis[3] = {
'X',
'Y',
'Z'};
397 for (
G4int j = 0; j <= 2; ++j)
400 G4int count = fBoundaries[j].size();
401 for(
G4int i=0; i < count-1; ++i)
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));
414 void G4SurfaceVoxelizer::BuildBoundingBox()
417 toleranceVector.
set(fTolerance/100,fTolerance/100,fTolerance/100);
418 for (
G4int i = 0; i <= 2; ++i)
420 G4double min = fBoundaries[i].front();
421 G4double max = fBoundaries[i].back();
422 fBoundingBoxSize[i] = (max-
min)/2;
423 fBoundingBoxCenter[i] = min + fBoundingBoxSize[i];
426 fBoundingBox =
G4Box(
"TessBBox", fBoundingBoxSize.
x(),
427 fBoundingBoxSize.
y(), fBoundingBoxSize.
z());
445 void G4SurfaceVoxelizer::SetReductionRatio(
G4int maxVoxels,
449 * fCandidatesCounts[1].size() * fCandidatesCounts[2].size();
451 if (maxVoxels > 0 && maxVoxels < maxTotal)
454 ratio = pow(ratio, 1./3.);
455 if (ratio > 1) { ratio = 1; }
456 reductionRatio.
set(ratio,ratio,ratio);
461 void G4SurfaceVoxelizer::BuildReduceVoxels(std::vector<G4double> boundaries[],
464 for (
G4int k = 0; k <= 2; ++k)
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;
476 voxel.
count = candidatesCount[j];
482 for (
G4int j = 0; j < max - 1; ++j) { voxelSet.insert(j); }
488 G4double reduction = reductionRatio[k];
491 G4int currentCount = voxelSet.size();
492 if (currentCount <= 2)
495 if (currentRatio <= reduction && currentCount <= 1000)
497 const G4int pos = *voxelSet.begin();
498 mergings.push_back(pos);
504 if (voxel.
next != max - 1) { voxelSet.erase(voxel.
next); }
511 if (voxel.
next != max - 1)
512 voxelSet.insert(voxel.
next);
524 sort(mergings.begin(), mergings.end());
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();
534 reducedBoundary[cur++] = boundary[j];
538 if (++i < (
G4int)mergings.size()) { skip = mergings[i] + 1; }
541 boundaries[k] = reducedBoundary;
547 void G4SurfaceVoxelizer::BuildReduceVoxels2(std::vector<G4double> boundaries[],
550 for (
G4int k = 0; k <= 2; ++k)
552 vector<G4int> &candidatesCount = fCandidatesCounts[k];
553 G4int max = candidatesCount.size();
555 for (
G4int i = 0; i <
max; i++) total += candidatesCount[i];
557 G4double reduction = reductionRatio[k];
561 G4int destination = (
G4int) (reduction * max) + 1;
562 if (destination > 1000) destination = 1000;
563 if (destination < 2) destination = 2;
566 vector<G4int> mergings;
568 vector<G4double> &boundary = boundaries[k];
569 vector<G4double> reducedBoundary(destination);
571 G4int sum = 0, cur = 0;
574 sum += candidatesCount[i];
575 if (sum > average * (cur + 1) || i == 0)
578 reducedBoundary[cur] = val;
580 if (cur == destination)
584 reducedBoundary[destination-1] = boundary[
max];
585 boundaries[k] = reducedBoundary;
590 void G4SurfaceVoxelizer::CreateMiniVoxels(std::vector<G4double> boundaries[],
593 vector<G4int> voxel(3), maxVoxels(3);
594 for (
G4int i = 0; i <= 2; ++i) maxVoxels[i] = boundaries[i].size();
597 for (voxel[2] = 0; voxel[2] < maxVoxels[2] - 1; ++voxel[2])
599 for (voxel[1] = 0; voxel[1] < maxVoxels[1] - 1; ++voxel[1])
601 for (voxel[0] = 0; voxel[0] < maxVoxels[0] - 1; ++voxel[0])
603 vector<G4int> candidates;
608 for (
G4int i = 0; i <= 2; ++i)
610 G4int index = voxel[i];
611 const vector<G4double> &boundary = boundaries[i];
614 box.
pos[i] = boundary[
index] + hlen;
616 fVoxelBoxes.push_back(box);
617 vector<G4int>(candidates).
swap(candidates);
618 fVoxelBoxesCandidates.push_back(candidates);
628 G4int maxVoxels = fMaxVoxels;
631 G4int size = facets.size();
634 for (
G4int i = 0; i < (
G4int) facets.size(); i++)
636 if (facets[i]->GetNumberOfVertices() > 3) size++;
640 if ((size >= 10 || maxVoxels > 0) && maxVoxels != 0 && maxVoxels != 1)
646 BuildVoxelLimits(facets);
658 BuildBitmasks(fBoundaries, 0,
true);
662 maxVoxels = fTotalCandidates;
663 if (fTotalCandidates > 1000000) maxVoxels = 1000000;
666 SetReductionRatio(maxVoxels, reductionRatio);
671 G4cout <<
"Total number of voxels: " << fCountOfVoxels <<
G4endl;
674 BuildReduceVoxels2(fBoundaries, reductionRatio);
679 G4cout <<
"Total number of voxels after reduction: "
680 << fCountOfVoxels <<
G4endl;
687 BuildBitmasks(fBoundaries, fBitmasks);
695 vector<G4double> miniBoundaries[3];
697 for (
G4int i = 0; i <= 2; ++i) { miniBoundaries[i] = fBoundaries[i]; }
699 G4int voxelsCountMini = (fCountOfVoxels >= 1000)
700 ? 100 : fCountOfVoxels / 10;
705 SetReductionRatio(voxelsCountMini, reductionRatioMini);
711 BuildReduceVoxels(miniBoundaries, reductionRatioMini);
715 G4cout <<
"Total number of mini voxels: " << total <<
G4endl;
722 BuildBitmasks(miniBoundaries, bitmasksMini);
728 CreateMiniVoxels(miniBoundaries, bitmasksMini);
743 G4cout <<
"Deallocating unnecessary fields during runtime..." <<
G4endl;
748 for (
G4int i = 0; i < 3; i++)
750 fCandidatesCounts[i].resize(0);
751 fBitmasks[i].
Clear();
763 G4cout <<
" Candidates in voxel [" << voxels[0] <<
" ; " << voxels[1]
764 <<
" ; " << voxels[2] <<
"]: ";
765 vector<G4int> candidates;
768 for (
G4int i = 0; i < count; ++i)
G4cout << candidates[i];
773 void G4SurfaceVoxelizer::FindComponentsFastest(
unsigned int mask,
774 std::vector<G4int> &list,
G4int i)
776 for (
G4int byte = 0; byte < (
G4int) (
sizeof(
unsigned int)); byte++)
778 if (
G4int maskByte = mask & 0xFF)
780 for (
G4int bit = 0; bit < 8; bit++)
783 { list.push_back(8*(
sizeof(
unsigned int)*i+ byte) + bit); }
784 if (!(maskByte >>= 1))
break;
793 std::vector<G4int> &list,
G4SurfBits *crossed)
const
799 for (
G4int i = 0; i <= 2; ++i)
801 if(point[i] < fBoundaries[i].front() || point[i] >= fBoundaries[i].back())
805 if (fTotalCandidates == 1)
816 if (!(mask = ((
unsigned int *) fBitmasks[0].fAllBits)[slice]
819 if (!(mask &= ((
unsigned int *) fBitmasks[1].fAllBits)[slice]
822 if (!(mask &= ((
unsigned int *) fBitmasks[2].fAllBits)[slice]
824 if (crossed && (!(mask &= ~((
unsigned int *)crossed->
fAllBits)[0])))
827 FindComponentsFastest(mask, list, 0);
831 unsigned int *masks[3], mask;
832 for (
G4int i = 0; i <= 2; ++i)
835 masks[i] = ((
unsigned int *) fBitmasks[i].fAllBits) + slice*fNPerSlice;
837 unsigned int *maskCrossed =
838 crossed ? (
unsigned int *)crossed->
fAllBits : 0;
840 for (
G4int i = 0 ; i < fNPerSlice; ++i)
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;
850 FindComponentsFastest(mask, list, i);
861 std::vector<G4int> &list,
866 if (fTotalCandidates == 1)
876 if (!(mask = ((
unsigned int *) bitmasks[0].fAllBits)[voxels[0]]
878 if (!(mask &= ((
unsigned int *) bitmasks[1].fAllBits)[voxels[1]]
880 if (!(mask &= ((
unsigned int *) bitmasks[2].fAllBits)[voxels[2]]
882 if (crossed && (!(mask &= ~((
unsigned int *)crossed->
fAllBits)[0])))
885 FindComponentsFastest(mask, list, 0);
889 unsigned int *masks[3], mask;
890 for (
G4int i = 0; i <= 2; ++i)
892 masks[i] = ((
unsigned int *) bitmasks[i].fAllBits)
893 + voxels[i]*fNPerSlice;
895 unsigned int *maskCrossed =
896 crossed ? (
unsigned int *)crossed->
fAllBits : 0;
898 for (
G4int i = 0 ; i < fNPerSlice; ++i)
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;
908 FindComponentsFastest(mask, list, i);
918 std::vector<G4int> &list,
G4SurfBits *crossed)
const
928 for (
G4int i = 0; i < 3; ++i)
930 if (point[i] < fBoundaries[i].front() || point[i] > fBoundaries[i].back())
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;
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);
985 const std::vector<G4int> &curVoxel)
const
989 for (
G4int i = 0; i <= 2; ++i)
993 const vector<G4double> &boundary = fBoundaries[i];
994 G4int cur = curVoxel[i];
995 if(direction[i] >= 1e-10)
997 if (boundary[++cur] - point[i] < fTolerance)
998 if (++cur >= (
G4int) boundary.size())
1003 if(direction[i] <= -1e-10)
1005 if (point[i] - boundary[cur] < fTolerance)
1012 G4double dif = boundary[cur] - point[i];
1013 G4double distance = dif / direction[i];
1015 if (shift > distance)
1026 std::vector<G4int> &curVoxel)
const
1028 for (
G4int i = 0; i <= 2; ++i)
1030 G4int index = curVoxel[i];
1031 const vector<G4double> &boundary = fBoundaries[i];
1033 if (direction[i] > 0)
1035 if (point[i] >= boundary[++index])
1036 if (++curVoxel[i] >= (
G4int) boundary.size() - 1)
1041 if (point[i] < boundary[index])
1042 if (--curVoxel[i] < 0)
1047 if (curVoxel[i] != indexOK)
1048 curVoxel[i] = indexOK;
1058 fReductionRatio.
set(0,0,0);
1065 fReductionRatio = ratioOfReduction;
1071 G4int res = fDefaultVoxelsCount;
1072 fDefaultVoxelsCount = count;
1079 return fDefaultVoxelsCount;
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());
1094 G4int csize = fCandidates.size();
1095 for (
G4int i = 0; i < csize; i++)
1097 size +=
sizeof(vector<G4int>) + fCandidates[i].capacity() *
sizeof(
G4int);
void set(double x, double y, double z)
G4bool Contains(const G4ThreeVector &point) const
void ResetBitNumber(unsigned int bitnumber)
void DisplayListNodes() const
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)
for(Int_t i=0;i< nentries;i++)
G4double DistanceToFirst(const G4ThreeVector &point, const G4ThreeVector &direction) const
G4double GetSurfaceTolerance() const
G4double DistanceToBoundingBox(const G4ThreeVector &point) const
G4double G4NeutronHPJENDLHEData::G4double result
G4int GetVoxelsIndex(G4int x, G4int y, G4int z) const
void DisplayVoxelLimits()
static G4double MinDistanceToBox(const G4ThreeVector &aPoint, const G4ThreeVector &f)
G4bool UpdateCurrentVoxel(const G4ThreeVector &point, const G4ThreeVector &direction, std::vector< G4int > &curVoxel) const
G4int GetBitsPerSlice() const
void SetBitNumber(unsigned int bitnumber, G4bool value=true)
static G4int BinarySearch(const std::vector< T > &vec, T value)
G4bool TestBitNumber(unsigned int bitnumber) const
G4GLOB_DLL std::ostream G4cout
G4double DistanceToNext(const G4ThreeVector &point, const G4ThreeVector &direction, const std::vector< G4int > &curVoxel) const
static G4SolidStore * GetInstance()
void swap(shared_ptr< P > &, shared_ptr< P > &)
void SetMaxVoxels(G4int max)
void ResetAllBits(G4bool value=false)
void set(unsigned int nbits, const char *array)
static G4int SetDefaultVoxelsCount(G4int count)
G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
G4double total(Particle const *const p1, Particle const *const p2)
virtual G4double Extent(const G4ThreeVector)=0
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void GetCandidatesVoxel(std::vector< G4int > &voxels)
long long CountVoxels(std::vector< G4double > boundaries[]) const
static G4GeometryTolerance * GetInstance()
static G4int GetDefaultVoxelsCount()
unsigned int GetNbytes() const