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