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