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