00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include <iostream>
00038 #include <iomanip>
00039 #include <sstream>
00040 #include <algorithm>
00041 #include <set>
00042
00043 #include "G4VSolid.hh"
00044
00045 #include "G4Orb.hh"
00046 #include "G4SurfaceVoxelizer.hh"
00047 #include "G4SolidStore.hh"
00048 #include "Randomize.hh"
00049 #include "G4PhysicalConstants.hh"
00050 #include "G4GeometryTolerance.hh"
00051 #include "G4CSGSolid.hh"
00052 #include "G4Types.hh"
00053 #include "geomdefs.hh"
00054
00055 using namespace std;
00056
00057 G4int G4SurfaceVoxelizer::fDefaultVoxelsCount = -1;
00058
00059
00060 G4SurfaceVoxelizer::G4SurfaceVoxelizer()
00061 : fBoundingBox("TessBBox", 1, 1, 1)
00062 {
00063 fCountOfVoxels = fNPerSlice = fTotalCandidates = 0;
00064
00065 fTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00066
00067 SetMaxVoxels(fDefaultVoxelsCount);
00068
00069 G4SolidStore::GetInstance()->DeRegister(&fBoundingBox);
00070 }
00071
00072
00073 G4SurfaceVoxelizer::~G4SurfaceVoxelizer()
00074 {
00075 }
00076
00077
00078 void G4SurfaceVoxelizer::BuildEmpty()
00079 {
00080
00081
00082
00083 vector<G4int> xyz(3), max(3), candidates(fTotalCandidates);
00084 const vector<G4int> empty(0);
00085
00086 for (G4int i = 0; i <= 2; i++) max[i] = fBoundaries[i].size();
00087 unsigned int size = max[0] * max[1] * max[2];
00088
00089 fEmpty.Clear();
00090 fEmpty.ResetBitNumber(size-1);
00091 fEmpty.ResetAllBits(true);
00092
00093 for (xyz[2] = 0; xyz[2] < max[2]; ++xyz[2])
00094 {
00095 for (xyz[1] = 0; xyz[1] < max[1]; ++xyz[1])
00096 {
00097 for (xyz[0] = 0; xyz[0] < max[0]; ++xyz[0])
00098 {
00099 G4int index = GetVoxelsIndex(xyz);
00100 if (GetCandidatesVoxelArray(xyz, candidates))
00101 {
00102 fEmpty.SetBitNumber(index, false);
00103
00104
00105
00106
00107
00108 vector<G4int> &c = (fCandidates[index] = empty);
00109 c.reserve(candidates.size());
00110 c.assign(candidates.begin(), candidates.end());
00111 }
00112 }
00113 }
00114 }
00115 #ifdef G4SPECSDEBUG
00116 G4cout << "Non-empty voxels count: " << fCandidates.size() << G4endl;
00117 #endif
00118 }
00119
00120
00121 void G4SurfaceVoxelizer::BuildVoxelLimits(std::vector<G4VFacet *> &facets)
00122 {
00123
00124
00125
00126
00127
00128 if (G4int numNodes = facets.size())
00129 {
00130 fBoxes.resize(numNodes);
00131 fNPerSlice = 1+(fBoxes.size()-1)/(8*sizeof(unsigned int));
00132
00133 G4ThreeVector toleranceVector;
00134 toleranceVector.set(10 * fTolerance,10 * fTolerance,10 * fTolerance);
00135
00136 for (G4int i = 0; i < numNodes; ++i)
00137 {
00138 G4VFacet &facet = *facets[i];
00139 G4ThreeVector min, max;
00140 G4ThreeVector x(1,0,0), y(0,1,0), z(0,0,1);
00141 G4ThreeVector extent;
00142 max.set (facet.Extent(x), facet.Extent(y), facet.Extent(z));
00143 min.set (-facet.Extent(-x), -facet.Extent(-y), -facet.Extent(-z));
00144 min -= toleranceVector; max += toleranceVector;
00145 G4ThreeVector hlen = (max - min) / 2;
00146 fBoxes[i].hlen = hlen;
00147 fBoxes[i].pos = min + hlen;
00148 }
00149 fTotalCandidates = fBoxes.size();
00150 }
00151 }
00152
00153
00154 void G4SurfaceVoxelizer::DisplayVoxelLimits()
00155 {
00156
00157
00158 G4int numNodes = fBoxes.size();
00159 G4int oldprec = G4cout.precision(16);
00160 for(G4int i = 0; i < numNodes; ++i)
00161 {
00162 G4cout << setw(10) << setiosflags(ios::fixed) <<
00163 " -> Node " << i+1 << ":\n" <<
00164 "\t * [x,y,z] = " << fBoxes[i].hlen <<
00165 "\t * [x,y,z] = " << fBoxes[i].pos << "\n";
00166 }
00167 G4cout.precision(oldprec);
00168 }
00169
00170
00171 void G4SurfaceVoxelizer::CreateSortedBoundary(std::vector<G4double> &boundary,
00172 G4int axis)
00173 {
00174
00175
00176
00177
00178 G4int numNodes = fBoxes.size();
00179
00180
00181
00182 for(G4int i = 0 ; i < numNodes; ++i)
00183 {
00184
00185
00186
00187 G4double p = fBoxes[i].pos[axis], d = fBoxes[i].hlen[axis];
00188
00189
00190
00191 #ifdef G4SPECSDEBUG
00192 G4cout << "Boundary " << p - d << " - " << p + d << G4endl;
00193 #endif
00194 boundary[2*i] = p - d;
00195 boundary[2*i+1] = p + d;
00196 }
00197 sort(boundary.begin(), boundary.end());
00198 }
00199
00200
00201 void G4SurfaceVoxelizer::BuildBoundaries()
00202 {
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213 if (G4int numNodes = fBoxes.size())
00214 {
00215 const G4double tolerance = fTolerance / 100.0;
00216
00217
00218 vector<G4double> sortedBoundary(2*numNodes);
00219
00220 G4int considered;
00221
00222 for (G4int j = 0; j <= 2; ++j)
00223 {
00224 CreateSortedBoundary(sortedBoundary, j);
00225 vector<G4double> &boundary = fBoundaries[j];
00226 boundary.clear();
00227
00228 considered = 0;
00229
00230 for(G4int i = 0 ; i < 2*numNodes; ++i)
00231 {
00232 G4double newBoundary = sortedBoundary[i];
00233 #ifdef G4SPECSDEBUG
00234 if (j == 0) G4cout << "Examining " << newBoundary << "..." << G4endl;
00235 #endif
00236 G4int size = boundary.size();
00237 if(!size || abs(boundary[size-1] - newBoundary) > tolerance)
00238 {
00239 considered++;
00240 {
00241 #ifdef G4SPECSDEBUG
00242 if (j == 0) G4cout << "Adding boundary " << newBoundary << "..."
00243 << G4endl;
00244 #endif
00245 boundary.push_back(newBoundary);
00246 continue;
00247 }
00248 }
00249
00250
00251 }
00252
00253 G4int n = boundary.size();
00254 G4int max = 100000;
00255 if (n > max/2)
00256 {
00257 G4int skip = n / (max /2);
00258
00259 vector<G4double> reduced;
00260 for (G4int i = 0; i < n; i++)
00261 {
00262
00263 G4int size = boundary.size();
00264 if (i % skip == 0 || i == 0 || i == size - 1)
00265 {
00266
00267
00268
00269
00270
00271 reduced.push_back(boundary[i]);
00272 }
00273 }
00274 boundary = reduced;
00275 }
00276 }
00277 }
00278 }
00279
00280
00281 void G4SurfaceVoxelizer::DisplayBoundaries()
00282 {
00283 char axis[3] = {'X', 'Y', 'Z'};
00284 for (G4int i = 0; i <= 2; ++i)
00285 {
00286 G4cout << " * " << axis[i] << " axis:" << G4endl << " | ";
00287 DisplayBoundaries(fBoundaries[i]);
00288 }
00289 }
00290
00291
00292 void G4SurfaceVoxelizer::DisplayBoundaries(std::vector<G4double> &boundaries)
00293 {
00294
00295
00296 G4int count = boundaries.size();
00297 G4int oldprec = G4cout.precision(16);
00298 for(G4int i = 0; i < count; ++i)
00299 {
00300 G4cout << setw(10) << setiosflags(ios::fixed) << boundaries[i];
00301 if(i != count-1) G4cout << "-> ";
00302 }
00303 G4cout << "|" << G4endl << "Number of boundaries: " << count << G4endl;
00304 G4cout.precision(oldprec);
00305 }
00306
00307
00308 void G4SurfaceVoxelizer::BuildBitmasks(std::vector<G4double> boundaries[],
00309 G4SurfBits bitmasks[], G4bool countsOnly)
00310 {
00311
00312
00313
00314 G4int numNodes = fBoxes.size();
00315 G4int bitsPerSlice = GetBitsPerSlice();
00316
00317 for (G4int k = 0; k < 3; k++)
00318 {
00319 G4int total = 0;
00320 vector<G4double> &boundary = boundaries[k];
00321 G4int voxelsCount = boundary.size() - 1;
00322 G4SurfBits &bitmask = bitmasks[k];
00323
00324 if (!countsOnly)
00325 {
00326 bitmask.Clear();
00327 #ifdef G4SPECSDEBUG
00328 G4cout << "Allocating bitmask..." << G4endl;
00329 #endif
00330 bitmask.SetBitNumber(voxelsCount*bitsPerSlice-1, false);
00331
00332
00333 }
00334 vector<G4int> &candidatesCount = fCandidatesCounts[k];
00335 candidatesCount.resize(voxelsCount);
00336
00337 for(G4int i = 0 ; i < voxelsCount; ++i) { candidatesCount[i] = 0; }
00338
00339
00340
00341 for(G4int j = 0 ; j < numNodes; j++)
00342 {
00343
00344
00345
00346 G4double p = fBoxes[j].pos[k], d = fBoxes[j].hlen[k];
00347
00348 G4double min = p - d;
00349 G4double max = p + d;
00350
00351 G4int i = BinarySearch(boundary, min);
00352 if (i < 0) { i = 0; }
00353
00354 do
00355 {
00356 if (!countsOnly)
00357 {
00358 G4int index = i*bitsPerSlice+j;
00359 bitmask.SetBitNumber(index);
00360 }
00361 candidatesCount[i]++;
00362 total++;
00363 i++;
00364 }
00365 while (max > boundary[i] && i < voxelsCount);
00366 }
00367 }
00368 #ifdef G4SPECSDEBUG
00369 G4cout << "Build list nodes completed." << G4endl;
00370 #endif
00371 }
00372
00373
00374 G4String G4SurfaceVoxelizer::GetCandidatesAsString(const G4SurfBits &bits)
00375 {
00376
00377
00378 stringstream ss;
00379 G4int numNodes = fBoxes.size();
00380
00381 for(G4int i=0; i<numNodes; ++i)
00382 {
00383 if (bits.TestBitNumber(i)) { ss << i+1 << " "; }
00384 }
00385 return ss.str();
00386 }
00387
00388
00389 void G4SurfaceVoxelizer::DisplayListNodes()
00390 {
00391
00392
00393 char axis[3] = {'X', 'Y', 'Z'};
00394 G4int size=8*sizeof(G4int)*fNPerSlice;
00395 G4SurfBits bits(size);
00396
00397 for (G4int j = 0; j <= 2; ++j)
00398 {
00399 G4cout << " * " << axis[j] << " axis:" << G4endl;
00400 G4int count = fBoundaries[j].size();
00401 for(G4int i=0; i < count-1; ++i)
00402 {
00403 G4cout << " Slice #" << i+1 << ": [" << fBoundaries[j][i]
00404 << " ; " << fBoundaries[j][i+1] << "] -> ";
00405 bits.set(size,(const char *)fBitmasks[j].fAllBits+i
00406 *fNPerSlice*sizeof(G4int));
00407 G4String result = GetCandidatesAsString(bits);
00408 G4cout << "[ " << result.c_str() << "] " << G4endl;
00409 }
00410 }
00411 }
00412
00413
00414 void G4SurfaceVoxelizer::BuildBoundingBox()
00415 {
00416 G4ThreeVector toleranceVector;
00417 toleranceVector.set(fTolerance/100,fTolerance/100,fTolerance/100);
00418 for (G4int i = 0; i <= 2; ++i)
00419 {
00420 G4double min = fBoundaries[i].front();
00421 G4double max = fBoundaries[i].back();
00422 fBoundingBoxSize[i] = (max-min)/2;
00423 fBoundingBoxCenter[i] = min + fBoundingBoxSize[i];
00424 }
00425
00426 fBoundingBox = G4Box("TessBBox", fBoundingBoxSize.x(),
00427 fBoundingBoxSize.y(), fBoundingBoxSize.z());
00428 }
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445 void G4SurfaceVoxelizer::SetReductionRatio(G4int maxVoxels,
00446 G4ThreeVector &reductionRatio)
00447 {
00448 G4double maxTotal = (G4double) fCandidatesCounts[0].size()
00449 * fCandidatesCounts[1].size() * fCandidatesCounts[2].size();
00450
00451 if (maxVoxels > 0 && maxVoxels < maxTotal)
00452 {
00453 G4double ratio = (G4double) maxVoxels / maxTotal;
00454 ratio = pow(ratio, 1./3.);
00455 if (ratio > 1) { ratio = 1; }
00456 reductionRatio.set(ratio,ratio,ratio);
00457 }
00458 }
00459
00460
00461 void G4SurfaceVoxelizer::BuildReduceVoxels(std::vector<G4double> boundaries[],
00462 G4ThreeVector reductionRatio)
00463 {
00464 for (G4int k = 0; k <= 2; ++k)
00465 {
00466 vector<G4int> &candidatesCount = fCandidatesCounts[k];
00467 G4int max = candidatesCount.size();
00468 vector<G4VoxelInfo> voxels(max);
00469 G4VoxelComparator comp(voxels);
00470 set<G4int, G4VoxelComparator> voxelSet(comp);
00471 vector<G4int> mergings;
00472
00473 for (G4int j = 0; j < max; ++j)
00474 {
00475 G4VoxelInfo &voxel = voxels[j];
00476 voxel.count = candidatesCount[j];
00477 voxel.previous = j - 1;
00478 voxel.next = j + 1;
00479 voxels[j] = voxel;
00480 }
00481
00482 for (G4int j = 0; j < max - 1; ++j) { voxelSet.insert(j); }
00483
00484
00485 G4int count = 0;
00486 while (true)
00487 {
00488 G4double reduction = reductionRatio[k];
00489 if (reduction == 0)
00490 break;
00491 G4int currentCount = voxelSet.size();
00492 if (currentCount <= 2)
00493 break;
00494 G4double currentRatio = 1 - (G4double) count / max;
00495 if (currentRatio <= reduction && currentCount <= 1000)
00496 break;
00497 const G4int pos = *voxelSet.begin();
00498 mergings.push_back(pos);
00499
00500 G4VoxelInfo &voxel = voxels[pos];
00501 G4VoxelInfo &nextVoxel = voxels[voxel.next];
00502
00503 voxelSet.erase(pos);
00504 if (voxel.next != max - 1) { voxelSet.erase(voxel.next); }
00505 if (voxel.previous != -1) { voxelSet.erase(voxel.previous); }
00506
00507 nextVoxel.count += voxel.count;
00508 voxel.count = 0;
00509 nextVoxel.previous = voxel.previous;
00510
00511 if (voxel.next != max - 1)
00512 voxelSet.insert(voxel.next);
00513
00514 if (voxel.previous != -1)
00515 {
00516 voxels[voxel.previous].next = voxel.next;
00517 voxelSet.insert(voxel.previous);
00518 }
00519 count++;
00520 }
00521
00522 if (mergings.size())
00523 {
00524 sort(mergings.begin(), mergings.end());
00525
00526 vector<G4double> &boundary = boundaries[k];
00527 vector<G4double> reducedBoundary(boundary.size() - mergings.size());
00528 G4int skip = mergings[0] + 1, cur = 0, i = 0;
00529 max = boundary.size();
00530 for (G4int j = 0; j < max; ++j)
00531 {
00532 if (j != skip)
00533 {
00534 reducedBoundary[cur++] = boundary[j];
00535 }
00536 else
00537 {
00538 if (++i < (G4int)mergings.size()) { skip = mergings[i] + 1; }
00539 }
00540 }
00541 boundaries[k] = reducedBoundary;
00542 }
00543 }
00544 }
00545
00546
00547 void G4SurfaceVoxelizer::BuildReduceVoxels2(std::vector<G4double> boundaries[],
00548 G4ThreeVector reductionRatio)
00549 {
00550 for (G4int k = 0; k <= 2; ++k)
00551 {
00552 vector<G4int> &candidatesCount = fCandidatesCounts[k];
00553 G4int max = candidatesCount.size();
00554 G4int total = 0;
00555 for (G4int i = 0; i < max; i++) total += candidatesCount[i];
00556
00557 G4double reduction = reductionRatio[k];
00558 if (reduction == 0)
00559 break;
00560
00561 G4int destination = (G4int) (reduction * max) + 1;
00562 if (destination > 1000) destination = 1000;
00563 if (destination < 2) destination = 2;
00564 G4double average = ((G4double)total / max) / reduction;
00565
00566 vector<G4int> mergings;
00567
00568 vector<G4double> &boundary = boundaries[k];
00569 vector<G4double> reducedBoundary(destination);
00570
00571 G4int sum = 0, cur = 0;
00572 for (G4int i = 0; i < max; i++)
00573 {
00574 sum += candidatesCount[i];
00575 if (sum > average * (cur + 1) || i == 0)
00576 {
00577 G4double val = boundary[i];
00578 reducedBoundary[cur] = val;
00579 cur++;
00580 if (cur == destination)
00581 break;
00582 }
00583 }
00584 reducedBoundary[destination-1] = boundary[max];
00585 boundaries[k] = reducedBoundary;
00586 }
00587 }
00588
00589
00590 void G4SurfaceVoxelizer::CreateMiniVoxels(std::vector<G4double> boundaries[],
00591 G4SurfBits bitmasks[])
00592 {
00593 vector<G4int> voxel(3), maxVoxels(3);
00594 for (G4int i = 0; i <= 2; ++i) maxVoxels[i] = boundaries[i].size();
00595
00596 G4ThreeVector point;
00597 for (voxel[2] = 0; voxel[2] < maxVoxels[2] - 1; ++voxel[2])
00598 {
00599 for (voxel[1] = 0; voxel[1] < maxVoxels[1] - 1; ++voxel[1])
00600 {
00601 for (voxel[0] = 0; voxel[0] < maxVoxels[0] - 1; ++voxel[0])
00602 {
00603 vector<G4int> candidates;
00604 if (GetCandidatesVoxelArray(voxel, bitmasks, candidates, 0))
00605 {
00606
00607 G4VoxelBox box;
00608 for (G4int i = 0; i <= 2; ++i)
00609 {
00610 G4int index = voxel[i];
00611 const vector<G4double> &boundary = boundaries[i];
00612 G4double hlen = 0.5 * (boundary[index+1] - boundary[index]);
00613 box.hlen[i] = hlen;
00614 box.pos[i] = boundary[index] + hlen;
00615 }
00616 fVoxelBoxes.push_back(box);
00617 vector<G4int>(candidates).swap(candidates);
00618 fVoxelBoxesCandidates.push_back(candidates);
00619 }
00620 }
00621 }
00622 }
00623 }
00624
00625
00626 void G4SurfaceVoxelizer::Voxelize(std::vector<G4VFacet *> &facets)
00627 {
00628 G4int maxVoxels = fMaxVoxels;
00629 G4ThreeVector reductionRatio = fReductionRatio;
00630
00631 G4int size = facets.size();
00632 if (size < 10)
00633 {
00634 for (G4int i = 0; i < (G4int) facets.size(); i++)
00635 {
00636 if (facets[i]->GetNumberOfVertices() > 3) size++;
00637 }
00638 }
00639
00640 if ((size >= 10 || maxVoxels > 0) && maxVoxels != 0 && maxVoxels != 1)
00641 {
00642 #ifdef G4SPECSDEBUG
00643 G4cout << "Building voxel limits..." << G4endl;
00644 #endif
00645
00646 BuildVoxelLimits(facets);
00647
00648 #ifdef G4SPECSDEBUG
00649 G4cout << "Building boundaries..." << G4endl;
00650 #endif
00651
00652 BuildBoundaries();
00653
00654 #ifdef G4SPECSDEBUG
00655 G4cout << "Building bitmasks..." << G4endl;
00656 #endif
00657
00658 BuildBitmasks(fBoundaries, 0, true);
00659
00660 if (maxVoxels < 0 && reductionRatio == G4ThreeVector())
00661 {
00662 maxVoxels = fTotalCandidates;
00663 if (fTotalCandidates > 1000000) maxVoxels = 1000000;
00664 }
00665
00666 SetReductionRatio(maxVoxels, reductionRatio);
00667
00668 fCountOfVoxels = CountVoxels(fBoundaries);
00669
00670 #ifdef G4SPECSDEBUG
00671 G4cout << "Total number of voxels: " << fCountOfVoxels << G4endl;
00672 #endif
00673
00674 BuildReduceVoxels2(fBoundaries, reductionRatio);
00675
00676 fCountOfVoxels = CountVoxels(fBoundaries);
00677
00678 #ifdef G4SPECSDEBUG
00679 G4cout << "Total number of voxels after reduction: "
00680 << fCountOfVoxels << G4endl;
00681 #endif
00682
00683 #ifdef G4SPECSDEBUG
00684 G4cout << "Building bitmasks..." << G4endl;
00685 #endif
00686
00687 BuildBitmasks(fBoundaries, fBitmasks);
00688
00689 G4ThreeVector reductionRatioMini;
00690
00691 G4SurfBits bitmasksMini[3];
00692
00693
00694
00695 vector<G4double> miniBoundaries[3];
00696
00697 for (G4int i = 0; i <= 2; ++i) { miniBoundaries[i] = fBoundaries[i]; }
00698
00699 G4int voxelsCountMini = (fCountOfVoxels >= 1000)
00700 ? 100 : fCountOfVoxels / 10;
00701
00702
00703
00704
00705 SetReductionRatio(voxelsCountMini, reductionRatioMini);
00706
00707 #ifdef G4SPECSDEBUG
00708 G4cout << "Building reduced voxels..." << G4endl;
00709 #endif
00710
00711 BuildReduceVoxels(miniBoundaries, reductionRatioMini);
00712
00713 #ifdef G4SPECSDEBUG
00714 G4int total = CountVoxels(miniBoundaries);
00715 G4cout << "Total number of mini voxels: " << total << G4endl;
00716 #endif
00717
00718 #ifdef G4SPECSDEBUG
00719 G4cout << "Building mini bitmasks..." << G4endl;
00720 #endif
00721
00722 BuildBitmasks(miniBoundaries, bitmasksMini);
00723
00724 #ifdef G4SPECSDEBUG
00725 G4cout << "Creating Mini Voxels..." << G4endl;
00726 #endif
00727
00728 CreateMiniVoxels(miniBoundaries, bitmasksMini);
00729
00730 #ifdef G4SPECSDEBUG
00731 G4cout << "Building bounding box..." << G4endl;
00732 #endif
00733
00734 BuildBoundingBox();
00735
00736 #ifdef G4SPECSDEBUG
00737 G4cout << "Building empty..." << G4endl;
00738 #endif
00739
00740 BuildEmpty();
00741
00742 #ifdef G4SPECSDEBUG
00743 G4cout << "Deallocating unnecessary fields during runtime..." << G4endl;
00744 #endif
00745
00746
00747 fBoxes.resize(0);
00748 for (G4int i = 0; i < 3; i++)
00749 {
00750 fCandidatesCounts[i].resize(0);
00751 fBitmasks[i].Clear();
00752 }
00753 }
00754 }
00755
00756
00757
00758 void G4SurfaceVoxelizer::GetCandidatesVoxel(std::vector<G4int> &voxels)
00759 {
00760
00761
00762
00763 G4cout << " Candidates in voxel [" << voxels[0] << " ; " << voxels[1]
00764 << " ; " << voxels[2] << "]: ";
00765 vector<G4int> candidates;
00766 G4int count = GetCandidatesVoxelArray(voxels, candidates);
00767 G4cout << "[ ";
00768 for (G4int i = 0; i < count; ++i) G4cout << candidates[i];
00769 G4cout << "] " << G4endl;
00770 }
00771
00772
00773 void G4SurfaceVoxelizer::FindComponentsFastest(unsigned int mask,
00774 std::vector<G4int> &list, G4int i)
00775 {
00776 for (G4int byte = 0; byte < (G4int) (sizeof(unsigned int)); byte++)
00777 {
00778 if (G4int maskByte = mask & 0xFF)
00779 {
00780 for (G4int bit = 0; bit < 8; bit++)
00781 {
00782 if (maskByte & 1)
00783 { list.push_back(8*(sizeof(unsigned int)*i+ byte) + bit); }
00784 if (!(maskByte >>= 1)) break;
00785 }
00786 }
00787 mask >>= 8;
00788 }
00789 }
00790
00791
00792 G4int G4SurfaceVoxelizer::GetCandidatesVoxelArray(const G4ThreeVector &point,
00793 std::vector<G4int> &list, G4SurfBits *crossed) const
00794 {
00795
00796
00797 list.clear();
00798
00799 for (G4int i = 0; i <= 2; ++i)
00800 {
00801 if(point[i] < fBoundaries[i].front() || point[i] >= fBoundaries[i].back())
00802 return 0;
00803 }
00804
00805 if (fTotalCandidates == 1)
00806 {
00807 list.push_back(0);
00808 return 1;
00809 }
00810 else
00811 {
00812 if (fNPerSlice == 1)
00813 {
00814 unsigned int mask;
00815 G4int slice = BinarySearch(fBoundaries[0], point.x());
00816 if (!(mask = ((unsigned int *) fBitmasks[0].fAllBits)[slice]
00817 )) return 0;
00818 slice = BinarySearch(fBoundaries[1], point.y());
00819 if (!(mask &= ((unsigned int *) fBitmasks[1].fAllBits)[slice]
00820 )) return 0;
00821 slice = BinarySearch(fBoundaries[2], point.z());
00822 if (!(mask &= ((unsigned int *) fBitmasks[2].fAllBits)[slice]
00823 )) return 0;
00824 if (crossed && (!(mask &= ~((unsigned int *)crossed->fAllBits)[0])))
00825 return 0;
00826
00827 FindComponentsFastest(mask, list, 0);
00828 }
00829 else
00830 {
00831 unsigned int *masks[3], mask;
00832 for (G4int i = 0; i <= 2; ++i)
00833 {
00834 G4int slice = BinarySearch(fBoundaries[i], point[i]);
00835 masks[i] = ((unsigned int *) fBitmasks[i].fAllBits) + slice*fNPerSlice;
00836 }
00837 unsigned int *maskCrossed =
00838 crossed ? (unsigned int *)crossed->fAllBits : 0;
00839
00840 for (G4int i = 0 ; i < fNPerSlice; ++i)
00841 {
00842
00843
00844
00845 if (!(mask = masks[0][i])) continue;
00846 if (!(mask &= masks[1][i])) continue;
00847 if (!(mask &= masks[2][i])) continue;
00848 if (maskCrossed && !(mask &= ~maskCrossed[i])) continue;
00849
00850 FindComponentsFastest(mask, list, i);
00851 }
00852 }
00853 }
00854 return list.size();
00855 }
00856
00857
00858 G4int
00859 G4SurfaceVoxelizer::GetCandidatesVoxelArray(const std::vector<G4int> &voxels,
00860 const G4SurfBits bitmasks[],
00861 std::vector<G4int> &list,
00862 G4SurfBits *crossed) const
00863 {
00864 list.clear();
00865
00866 if (fTotalCandidates == 1)
00867 {
00868 list.push_back(0);
00869 return 1;
00870 }
00871 else
00872 {
00873 if (fNPerSlice == 1)
00874 {
00875 unsigned int mask;
00876 if (!(mask = ((unsigned int *) bitmasks[0].fAllBits)[voxels[0]]
00877 )) return 0;
00878 if (!(mask &= ((unsigned int *) bitmasks[1].fAllBits)[voxels[1]]
00879 )) return 0;
00880 if (!(mask &= ((unsigned int *) bitmasks[2].fAllBits)[voxels[2]]
00881 )) return 0;
00882 if (crossed && (!(mask &= ~((unsigned int *)crossed->fAllBits)[0])))
00883 return 0;
00884
00885 FindComponentsFastest(mask, list, 0);
00886 }
00887 else
00888 {
00889 unsigned int *masks[3], mask;
00890 for (G4int i = 0; i <= 2; ++i)
00891 {
00892 masks[i] = ((unsigned int *) bitmasks[i].fAllBits)
00893 + voxels[i]*fNPerSlice;
00894 }
00895 unsigned int *maskCrossed =
00896 crossed ? (unsigned int *)crossed->fAllBits : 0;
00897
00898 for (G4int i = 0 ; i < fNPerSlice; ++i)
00899 {
00900
00901
00902
00903 if (!(mask = masks[0][i])) continue;
00904 if (!(mask &= masks[1][i])) continue;
00905 if (!(mask &= masks[2][i])) continue;
00906 if (maskCrossed && !(mask &= ~maskCrossed[i])) continue;
00907
00908 FindComponentsFastest(mask, list, i);
00909 }
00910 }
00911 }
00912 return list.size();
00913 }
00914
00915
00916 G4int
00917 G4SurfaceVoxelizer::GetCandidatesVoxelArray(const std::vector<G4int> &voxels,
00918 std::vector<G4int> &list, G4SurfBits *crossed) const
00919 {
00920
00921
00922 return GetCandidatesVoxelArray(voxels, fBitmasks, list, crossed);
00923 }
00924
00925
00926 G4bool G4SurfaceVoxelizer::Contains(const G4ThreeVector &point) const
00927 {
00928 for (G4int i = 0; i < 3; ++i)
00929 {
00930 if (point[i] < fBoundaries[i].front() || point[i] > fBoundaries[i].back())
00931 return false;
00932 }
00933 return true;
00934 }
00935
00936
00937 G4double
00938 G4SurfaceVoxelizer::DistanceToFirst(const G4ThreeVector &point,
00939 const G4ThreeVector &direction) const
00940 {
00941 G4ThreeVector pointShifted = point - fBoundingBoxCenter;
00942 G4double shift = fBoundingBox.DistanceToIn(pointShifted, direction);
00943 return shift;
00944 }
00945
00946
00947 G4double
00948 G4SurfaceVoxelizer::DistanceToBoundingBox(const G4ThreeVector &point) const
00949 {
00950 G4ThreeVector pointShifted = point - fBoundingBoxCenter;
00951 G4double shift = MinDistanceToBox(pointShifted, fBoundingBoxSize);
00952 return shift;
00953 }
00954
00955
00956 G4double
00957 G4SurfaceVoxelizer::MinDistanceToBox (const G4ThreeVector &aPoint,
00958 const G4ThreeVector &f)
00959 {
00960
00961
00962
00963
00964 G4double safe, safx, safy, safz;
00965 safe = safx = -f.x() + std::abs(aPoint.x());
00966 safy = -f.y() + std::abs(aPoint.y());
00967 if ( safy > safe ) safe = safy;
00968 safz = -f.z() + std::abs(aPoint.z());
00969 if ( safz > safe ) safe = safz;
00970 if (safe < 0.0) return 0.0;
00971
00972 G4double safsq = 0.0;
00973 G4int count = 0;
00974 if ( safx > 0 ) { safsq += safx*safx; count++; }
00975 if ( safy > 0 ) { safsq += safy*safy; count++; }
00976 if ( safz > 0 ) { safsq += safz*safz; count++; }
00977 if (count == 1) return safe;
00978 return std::sqrt(safsq);
00979 }
00980
00981
00982 G4double
00983 G4SurfaceVoxelizer::DistanceToNext(const G4ThreeVector &point,
00984 const G4ThreeVector &direction,
00985 const std::vector<G4int> &curVoxel) const
00986 {
00987 G4double shift = kInfinity;
00988
00989 for (G4int i = 0; i <= 2; ++i)
00990 {
00991
00992
00993 const vector<G4double> &boundary = fBoundaries[i];
00994 G4int cur = curVoxel[i];
00995 if(direction[i] >= 1e-10)
00996 {
00997 if (boundary[cur] - point[i] < fTolerance)
00998 if (++cur >= (G4int) boundary.size())
00999 continue;
01000 }
01001 else
01002 {
01003 if(direction[i] <= -1e-10)
01004 {
01005 if (point[i] - boundary[cur] < fTolerance)
01006 if (--cur < 0)
01007 continue;
01008 }
01009 else
01010 continue;
01011 }
01012 G4double dif = boundary[cur] - point[i];
01013 G4double distance = dif / direction[i];
01014
01015 if (shift > distance)
01016 shift = distance;
01017 }
01018
01019 return shift;
01020 }
01021
01022
01023 G4bool
01024 G4SurfaceVoxelizer::UpdateCurrentVoxel(const G4ThreeVector &point,
01025 const G4ThreeVector &direction,
01026 std::vector<G4int> &curVoxel) const
01027 {
01028 for (G4int i = 0; i <= 2; ++i)
01029 {
01030 G4int index = curVoxel[i];
01031 const vector<G4double> &boundary = fBoundaries[i];
01032
01033 if (direction[i] > 0)
01034 {
01035 if (point[i] >= boundary[++index])
01036 if (++curVoxel[i] >= (G4int) boundary.size() - 1)
01037 return false;
01038 }
01039 else
01040 {
01041 if (point[i] < boundary[index])
01042 if (--curVoxel[i] < 0)
01043 return false;
01044 }
01045 #ifdef G4SPECSDEBUG
01046 G4int indexOK = BinarySearch(boundary, point[i]);
01047 if (curVoxel[i] != indexOK)
01048 curVoxel[i] = indexOK;
01049 #endif
01050 }
01051 return true;
01052 }
01053
01054
01055 void G4SurfaceVoxelizer::SetMaxVoxels(G4int max)
01056 {
01057 fMaxVoxels = max;
01058 fReductionRatio.set(0,0,0);
01059 }
01060
01061
01062 void G4SurfaceVoxelizer::SetMaxVoxels(const G4ThreeVector &ratioOfReduction)
01063 {
01064 fMaxVoxels = -1;
01065 fReductionRatio = ratioOfReduction;
01066 }
01067
01068
01069 G4int G4SurfaceVoxelizer::SetDefaultVoxelsCount(G4int count)
01070 {
01071 G4int res = fDefaultVoxelsCount;
01072 fDefaultVoxelsCount = count;
01073 return res;
01074 }
01075
01076
01077 G4int G4SurfaceVoxelizer::GetDefaultVoxelsCount()
01078 {
01079 return fDefaultVoxelsCount;
01080 }
01081
01082
01083 G4int G4SurfaceVoxelizer::AllocatedMemory()
01084 {
01085 G4int size = fEmpty.GetNbytes();
01086 size += fBoxes.capacity() * sizeof(G4VoxelBox);
01087 size += sizeof(G4double) * (fBoundaries[0].capacity()
01088 + fBoundaries[1].capacity() + fBoundaries[2].capacity());
01089 size += sizeof(G4int) * (fCandidatesCounts[0].capacity()
01090 + fCandidatesCounts[1].capacity() + fCandidatesCounts[2].capacity());
01091 size += fBitmasks[0].GetNbytes() + fBitmasks[1].GetNbytes()
01092 + fBitmasks[2].GetNbytes();
01093
01094 G4int csize = fCandidates.size();
01095 for (G4int i = 0; i < csize; i++)
01096 {
01097 size += sizeof(vector<G4int>) + fCandidates[i].capacity() * sizeof(G4int);
01098 }
01099
01100 return size;
01101 }
01102