moleculardna/src/OctreeNode.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 #include "OctreeNode.hh"
0028
0029 #include "G4VPhysicalVolume.hh"
0030
0031 #include <cassert>
0032 #include <utility>
0033
0034 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0035
0036 OctreeNode::OctreeNode(const G4ThreeVector& position, const G4ThreeVector& halfLengths,
0037                        G4int maxContents, OctreeNode* parent)
0038   : fPosition(position),
0039     fHalfLengths(halfLengths),
0040     fMaxContents(maxContents),
0041     fParent(parent),
0042     fChildren(std::array<OctreeNode*, 8>())
0043 {
0044   fHalfLengthsMag = fHalfLengths.mag();
0045 }
0046
0047 OctreeNode::~OctreeNode()
0048 {
0049   // Delete children
0050   for (auto& it : fChildren) {
0051     delete it;
0052   }
0053 }
0054
0055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0056
0057 const OctreeNode* OctreeNode::GetChildFromPosition(const G4ThreeVector& pos) const
0058 {
0059   G4ThreeVector localPosition = pos - fPosition;
0060
0061   // Get the right child
0062   // The scheme is defined by quadrant as follows:
0063   // Zero is treated as positive
0064   // X Y Z Index | X Y Z Index
0065   // + + +   0   | - + +   4
0066   // + + -   1   | - + -   5
0067   // + - +   2   | - - +   6
0068   // + - -   3   | - - -   7
0069   int bit2 = (localPosition.getX() < 0) ? 4 : 0;
0070   int bit1 = (localPosition.getY() < 0) ? 2 : 0;
0071   int bit0 = (localPosition.getZ() < 0) ? 1 : 0;
0072   return fChildren[bit0 + bit1 + bit2];
0073 }
0074
0075 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0076
0077 // Non-const version of above function. Ugly but useful in this case
0078 OctreeNode* OctreeNode::GetChildFromPosition(const G4ThreeVector& pos)
0079 {
0080   G4ThreeVector localPosition = pos - fPosition;
0081
0082   // Get the right child
0083   // The scheme is defined by quadrant as follows:
0084   // Zero is treated as positive
0085   // X Y Z Index | X Y Z Index
0086   // + + +   0   | - + +   4
0087   // + + -   1   | - + -   5
0088   // + - +   2   | - - +   6
0089   // + - -   3   | - - -   7
0090   int bit2 = (localPosition.getX() < 0) ? 4 : 0;
0091   int bit1 = (localPosition.getY() < 0) ? 2 : 0;
0092   int bit0 = (localPosition.getZ() < 0) ? 1 : 0;
0093   int idx = bit0 + bit1 + bit2;
0094
0095   assert(idx < (int)fChildren.size());
0096   return fChildren[idx];
0097 }
0098
0099 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0100
0101 void OctreeNode::AddPhysicalVolume(G4VPhysicalVolume* pv)
0102 {
0103   // Consider adding test for the validity of the PV position
0104   if (this->HasChildren()) {
0105     OctreeNode* child = this->GetChildFromPosition(pv->GetTranslation());
0106     child->AddPhysicalVolume(pv);
0107   }
0108   else {
0109     // if there are fMaxContents elements in the bin, we need to split
0110     // the bin before adding a new element.
0111     if ((G4int)fContents.size() == fMaxContents) {
0112       this->Split();
0113       this->AddPhysicalVolume(pv);
0114     }
0115     else {
0116       fContents.push_back(pv);
0117     }
0118   }
0119 }
0120
0121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0122
0123 void OctreeNode::Split()
0124 {
0125   G4ThreeVector hl = 0.5 * fHalfLengths;
0126   G4double xp = fPosition.getX();
0127   G4double yp = fPosition.getY();
0128   G4double zp = fPosition.getZ();
0129   G4double xl = hl.getX();
0130   G4double yl = hl.getY();
0131   G4double zl = hl.getZ();
0132   G4ThreeVector newpos;
0133
0134   // Construct children
0135   newpos = G4ThreeVector(xp + xl, yp + yl, zp + zl);
0136   fChildren[0] = new OctreeNode(newpos, hl, fMaxContents, this);
0137   newpos = G4ThreeVector(xp + xl, yp + yl, zp - zl);
0138   fChildren[1] = new OctreeNode(newpos, hl, fMaxContents, this);
0139   newpos = G4ThreeVector(xp + xl, yp - yl, zp + zl);
0140   fChildren[2] = new OctreeNode(newpos, hl, fMaxContents, this);
0141   newpos = G4ThreeVector(xp + xl, yp - yl, zp - zl);
0142   fChildren[3] = new OctreeNode(newpos, hl, fMaxContents, this);
0143   newpos = G4ThreeVector(xp - xl, yp + yl, zp + zl);
0144   fChildren[4] = new OctreeNode(newpos, hl, fMaxContents, this);
0145   newpos = G4ThreeVector(xp - xl, yp + yl, zp - zl);
0146   fChildren[5] = new OctreeNode(newpos, hl, fMaxContents, this);
0147   newpos = G4ThreeVector(xp - xl, yp - yl, zp + zl);
0148   fChildren[6] = new OctreeNode(newpos, hl, fMaxContents, this);
0149   newpos = G4ThreeVector(xp - xl, yp - yl, zp - zl);
0150   fChildren[7] = new OctreeNode(newpos, hl, fMaxContents, this);
0151
0152   // Distribute contents to children
0153
0154   for (int i = fContents.size() - 1; i >= 0; --i) {
0155     G4VPhysicalVolume* pv = fContents[i];
0156     OctreeNode* child = this->GetChildFromPosition(pv->GetTranslation());
0157     assert(child != this);
0158     child->AddPhysicalVolume(pv);
0159   }
0160   fContents.clear();
0161 }
0162
0163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0164
0165 std::vector<G4VPhysicalVolume*> OctreeNode::GetContents() const
0166 {
0167   if (this->HasChildren())  // if has sub-nodes
0168   {
0169     std::vector<G4VPhysicalVolume*> vec;
0170     std::vector<G4VPhysicalVolume*> childCont;
0171     for (auto it = fChildren.begin(); it != fChildren.end(); ++it) {
0172       childCont = (*it)->GetContents();
0173       for (auto jt = childCont.begin(); jt != childCont.end(); ++jt) {
0174         vec.push_back(*jt);
0175       }
0176     }
0177     return vec;
0178   }
0179   else  // if no sub-nodes
0180   {
0181     return fContents;
0182   }
0183 }
0184
0185 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0186
0187 // Search for Octree nodes in a volume
0188 const std::vector<G4VPhysicalVolume*> OctreeNode::SearchOctree(const G4ThreeVector& pos,
0189                                                                G4double _rad) const
0190 {
0191   // Need to search based on absolute position of each volume rather than
0192   // relative  position, and maintain a list of candidate nodes and
0193   // final nodes within the radius
0194
0195   std::vector<const OctreeNode*> nodes;
0196   std::vector<const OctreeNode*> candidates;
0197   nodes.reserve(512);
0198   candidates.reserve(512);
0199
0200   if (this->HasChildren()) {
0201     for (auto it = this->GetChildren().begin(); it != this->GetChildren().end(); it++) {
0202       candidates.push_back(*it);
0203     }
0204   }
0205   else {
0206     candidates.push_back(this);
0207   }
0208
0209   const OctreeNode* aNode;
0210   G4double d;
0211   while (!candidates.empty()) {
0212     aNode = candidates.back();
0213     candidates.pop_back();
0214     d = (aNode->GetPosition() - pos).mag() - aNode->GetHalfLengthsMag();
0215     // if node within circle
0216     if (d < _rad) {
0217       if (aNode->HasChildren()) {
0218         for (auto it = aNode->GetChildren().begin(); it != aNode->GetChildren().end(); ++it) {
0219           candidates.push_back(*it);
0220         }
0221       }
0222       else {
0223         nodes.push_back(aNode);
0224       }
0225     }
0226   }
0227
0228   // G4cout << "Found " << nodes.size() << " nodes" << G4endl;
0229
0230   // Get the physical volumes
0231   G4double r2 = _rad * _rad;
0232   std::vector<G4VPhysicalVolume*> pvols;
0233   for (auto it = nodes.begin(); it != nodes.end(); ++it) {
0234     std::vector<G4VPhysicalVolume*> conts = (*it)->GetContents();
0235     for (auto jt = conts.begin(); jt != conts.end(); ++jt) {
0236       if ((pos - (*jt)->GetTranslation()).mag2() < r2) {
0237         pvols.push_back((*jt));
0238       }
0239     }
0240   }
0241   // G4cout << "Found " << pvols.size() << " candidates" << G4endl;
0242   return pvols;
0243 }
0244
0245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0246
0247 // Search for Octree nodes in a volume
0248 void OctreeNode::SearchOctree(const G4ThreeVector& pos, std::vector<G4VPhysicalVolume*>& out,
0249                               G4double _rad) const
0250 {
0251   // Need to search based on absolute position of each volume rather than
0252   // relative  position, and maintain a list of candidate nodes and
0253   // final nodes within the radius
0254
0255   std::vector<const OctreeNode*> nodes;
0256   std::vector<const OctreeNode*> candidates;
0257   nodes.reserve(512);
0258   candidates.reserve(512);
0259
0260   if (this->HasChildren()) {
0261     for (auto it = this->GetChildren().begin(); it != this->GetChildren().end(); ++it) {
0262       candidates.push_back(*it);
0263     }
0264   }
0265   else {
0266     candidates.push_back(this);
0267   }
0268
0269   const OctreeNode* aNode;
0270   G4double d;
0271   while (!candidates.empty()) {
0272     aNode = candidates.back();
0273     candidates.pop_back();
0274     d = (aNode->GetPosition() - pos).mag() - aNode->GetHalfLengthsMag();
0275     // if node within circle
0276     if (d < _rad) {
0277       if (aNode->HasChildren()) {
0278         for (auto it = aNode->GetChildren().begin(); it != aNode->GetChildren().end(); it++) {
0279           candidates.push_back(*it);
0280         }
0281       }
0282       else {
0283         nodes.push_back(aNode);
0284       }
0285     }
0286   }
0287
0288   // Get the physical volumes
0289   G4double r2 = _rad * _rad;
0290   for (auto it = nodes.begin(); it != nodes.end(); ++it) {
0291     std::vector<G4VPhysicalVolume*> conts = (*it)->GetContents();
0292     for (auto jt = conts.begin(); jt != conts.end(); ++jt) {
0293       if ((pos - (*jt)->GetTranslation()).mag2() < r2) {
0294         out.push_back((*jt));
0295       }
0296     }
0297   }
0298 }
0299
0300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0301
0302 const std::vector<G4VPhysicalVolume*> OctreeNode::SearchOctree(const G4ThreeVector& pos) const
0303 {
0304   const OctreeNode* child = GetChildFromPosition(pos);
0305   return child->GetContents();
0306 }
0307
0308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0309
0310 G4int OctreeNode::GetNumberOfTerminalNodes()
0311 {
0312   if (!this->HasChildren()) {
0313     return 1;
0314   }
0315   // sum the number of nodes of each child
0316   G4int n = 0;
0317   for (auto it = this->GetChildren().begin(); it != this->GetChildren().end(); ++it) {
0318     n = n + (*it)->GetNumberOfTerminalNodes();
0319   }
0320
0321   return n;
0322 }
0323 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......