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