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