include/Geant4/G4Octree.icc

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 // Author: HoangTRAN, 20/2/2019
0028
0029 #define OCTREE G4Octree<Iterator,Extractor,Point>
0030 #define OCTREE_TEMPLATE typename Iterator, class Extractor,typename Point
0031
0032 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0033
0034 template <OCTREE_TEMPLATE>
0035 G4ThreadLocal G4Allocator<OCTREE>* OCTREE::fgAllocator = nullptr;
0036
0037 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0038
0039 template <OCTREE_TEMPLATE>
0040 OCTREE::G4Octree()
0041     : functor_(Extractor())
0042     , head_(nullptr)
0043     , size_(0)
0044 {}
0045
0046 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0047
0048 template <OCTREE_TEMPLATE>
0049 OCTREE::G4Octree(Iterator begin, Iterator end)
0050     : G4Octree(begin,end,Extractor())
0051 {
0052     head_ = nullptr;
0053     size_ = 0;
0054 }
0055
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0057
0058 template <OCTREE_TEMPLATE>
0059 OCTREE::G4Octree(Iterator begin, Iterator end, Extractor f)
0060     : functor_(f),head_(nullptr),size_(0)
0061 {
0062     std::vector<std::pair<Iterator,Point>> v;
0063     for(auto it=begin;it!=end;++it)
0064     {
0065         v.push_back(std::pair<Iterator,Point>(it,functor_(it)));
0066     }
0067     size_ = v.size();
0068     head_ = new Node(v);
0069 }
0070
0071 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0072
0073 template <OCTREE_TEMPLATE>
0074 OCTREE::G4Octree(OCTREE::tree_type&& rhs)
0075     : functor_(rhs.functor_)
0076     , head_(rhs.head_)
0077     , size_(rhs.size_)
0078 {}
0079
0080 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0081
0082 template <OCTREE_TEMPLATE>
0083 void OCTREE::swap(OCTREE::tree_type& rhs)
0084 {
0085     std::swap(head_, rhs.head_);
0086     std::swap(functor_, rhs.functor_);
0087     std::swap(size_, rhs.size_);
0088 }
0089
0090 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0091
0092 template <OCTREE_TEMPLATE>
0093 typename OCTREE::tree_type& OCTREE::operator=(typename OCTREE::tree_type rhs)
0094 {
0095     swap(rhs);
0096     return *this;
0097 }
0098
0099 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0100
0101 template <OCTREE_TEMPLATE>
0102 typename OCTREE::tree_type& OCTREE::operator=(typename OCTREE::tree_type&& rhs)
0103 {
0104     swap(rhs);
0105     return *this;
0106 }
0107 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0108
0109 template <OCTREE_TEMPLATE>
0110 OCTREE::~G4Octree()
0111 {
0112     delete head_;
0113 }
0114
0115 template <OCTREE_TEMPLATE>
0116 size_t OCTREE::size() const
0117 {
0118     return size_;
0119 }
0120
0121 template <OCTREE_TEMPLATE>
0122 OCTREE::Node::Node(const NodeVector& input_values)
0123     : Node(input_values,
0124          G4DNABoundingBox(InnerIterator(input_values.begin()),
0125                         InnerIterator(input_values.end())),
0126             0)
0127 {}
0128
0129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0130
0131 template <OCTREE_TEMPLATE>
0132 OCTREE::Node::Node(
0133     const NodeVector& input_values,
0134     const G4DNABoundingBox& box,
0135     size_t current_depth):
0136     fpValue(nullptr),
0137     fBigVolume(box),
0138     fNodeType(DEFAULT)
0139 {
0140     if (current_depth > max_depth)
0141     {
0142         init_max_depth_leaf(input_values);
0143     }
0144     else if (input_values.size() <= max_per_node)
0145     {
0146         init_leaf(input_values);
0147     }
0148     else
0149     {
0150         init_internal(input_values, current_depth);
0151     }
0152 }
0153
0154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0155 template <OCTREE_TEMPLATE>
0156 OCTREE::Node::~Node()
0157 {
0158     if (fNodeType == NodeTypes::INTERNAL)
0159     {
0160         childNodeArray& children = *static_cast<childNodeArray*>(fpValue);
0161         for (size_t i = 0; i < 8; ++i)
0162         {
0163             if (children[i] != nullptr)
0164             {
0165                 delete children[i];
0166                 children[i] = nullptr;
0167             }
0168         }
0169         delete &children;
0170     }
0171     else if (fNodeType == NodeTypes::LEAF)
0172     {
0173         auto toDelete = static_cast<LeafValues*>(fpValue);
0174         toDelete->size_ = 0;
0175         delete static_cast<LeafValues*>(fpValue);
0176     }
0177     else if (fNodeType == NodeTypes::MAX_DEPTH_LEAF)
0178     {
0179         delete static_cast<NodeVector*>(fpValue);
0180     }
0181     fpValue = nullptr;
0182 }
0183
0184 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0185
0186 template <OCTREE_TEMPLATE>
0187 void OCTREE::Node::init_max_depth_leaf(
0188     const NodeVector& input_values)
0189 {
0190     fpValue = new NodeVector(input_values);
0191     fNodeType = NodeTypes::MAX_DEPTH_LEAF;
0192 }
0193 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0194 template <OCTREE_TEMPLATE>
0195 void OCTREE::Node::init_leaf(const NodeVector& input_values)
0196 {
0197     std::array<std::pair<Iterator, Point>, max_per_node> a;
0198     std::copy(input_values.begin(), input_values.end(), a.begin());
0199     fpValue = new LeafValues{a, input_values.size()};
0200     fNodeType = NodeTypes::LEAF;
0201 }
0202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0203 template <OCTREE_TEMPLATE>
0204 void OCTREE::Node::init_internal(
0205     const NodeVector& input_values,
0206     size_t current_depth)
0207 {
0208     std::array<NodeVector, 8> childVectors;
0209     std::array<G4DNABoundingBox, 8> boxes = fBigVolume.partition();
0210     std::array<Node*, 8> children{};
0211     for (size_t child = 0; child < 8; ++child)
0212     {
0213         NodeVector& childVector = childVectors[child];
0214         childVector.reserve(input_values.size()/8);
0215         std::copy_if(input_values.begin(),input_values.end(),
0216                     std::back_inserter(childVector),
0217                 [&boxes, child](const std::pair<Iterator, Point>& element)
0218                 -> G4bool
0219         {
0220             const Point& p = element.second;
0221             return boxes[child].contains(p);
0222         }
0223         );
0224         children[child] = childVector.empty()
0225             ? nullptr
0226             : new Node(childVector, boxes[child], ++current_depth);
0227     }
0228     fpValue = new std::array<Node*, 8>(children);
0229     fNodeType = NodeTypes::INTERNAL;
0230 }
0231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0232
0233 template <OCTREE_TEMPLATE>
0234 template <typename OutPutContainer>
0235 G4bool OCTREE::Node::radiusNeighbors(const Point& query,
0236                                     G4double radius,
0237                                     OutPutContainer& resultIndices) const
0238 {
0239     G4bool success = false;
0240     G4double distance = 0;
0241     if (fNodeType == NodeTypes::INTERNAL)
0242     {
0243         childNodeArray& children = *static_cast<childNodeArray*>(fpValue);
0244         for (auto eachChild : children)
0245         {
0246             if (eachChild == nullptr)
0247             {
0248                 continue;
0249             }
0250             if(!eachChild->fBigVolume.overlap(query,radius))
0251             {
0252                 continue;
0253             }
0254             success = eachChild->radiusNeighbors(query, radius, resultIndices) || success;
0255         }
0256     }
0257     else if (fNodeType == NodeTypes::LEAF)
0258     {
0259         if(fpValue != nullptr)
0260         {
0261             LeafValues& children = *static_cast<LeafValues*>(fpValue);
0262             for (size_t i = 0; i < children.size_; ++i)
0263             {
0264                 distance = (query - std::get<1>(children.values_[i])).mag();
0265
0266                 if(distance != 0)
0267                 {
0268                     if( distance < radius )//TODO: find another solution for this using boundingbox
0269                     {
0270                         resultIndices.push_back(std::make_pair(std::get<0>(children.values_[i]),distance));
0271                         success = true;
0272                     }
0273                 }
0274             }
0275         }
0276     }
0277     else if (fNodeType == NodeTypes::MAX_DEPTH_LEAF)
0278     {
0279         NodeVector& children = *static_cast<NodeVector*>(fpValue);
0280         for (auto & child : children)
0281         {
0282             const Point& point = std::get<1>(child);
0283             //if (this->fBigVolume.contains(query, point, radius))
0284             distance = (query - point).mag();
0285             if( distance == 0. )
0286             {
0287                 continue;
0288             }
0289             if( distance < radius )
0290             {
0291                 if(distance == 0)
0292                 {
0293                     throw std::runtime_error("distance == 0 => OCTREE::Node::radiusNeighbors : find itself");
0294                 }
0295                 Iterator resultIndex = std::get<0>(child);
0296                 resultIndices.push_back(std::make_pair(resultIndex,distance));
0297                 success = true;
0298             }
0299         }
0300     }
0301     else
0302     {
0303         throw std::runtime_error("fNodeType is not set : find itself");
0304     }
0305     return success;
0306 }
0307
0308 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0309
0310 template <OCTREE_TEMPLATE>
0311 template <typename OutPutContainer>
0312 void OCTREE::radiusNeighbors(const Point& query,
0313                              const G4double& radius,
0314                              OutPutContainer& resultIndices) const
0315 {
0316     resultIndices.clear();
0317     if (head_ == nullptr)
0318     {
0319          return;
0320     }
0321     head_->radiusNeighbors(query, radius, resultIndices);
0322 }
0323
0324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0325
0326 template <OCTREE_TEMPLATE>
0327 void* OCTREE::operator new(size_t)
0328 {
0329     if (!fgAllocator)
0330     {
0331         fgAllocator = new G4Allocator<OCTREE>;
0332     }
0333     return (void *) fgAllocator->MallocSingle();
0334 }
0335
0336 template <OCTREE_TEMPLATE>
0337 void OCTREE::operator delete(void *a)
0338 {
0339     fgAllocator->FreeSingle((OCTREE*)a);
0340 }