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