2// ********************************************************************
3// * License and Disclaimer *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
27// Author: HoangTRAN, 20/2/2019
29#define OCTREE G4Octree<Iterator,Extractor,Point>
30#define OCTREE_TEMPLATE typename Iterator, class Extractor,typename Point
32//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
34template <OCTREE_TEMPLATE>
35G4ThreadLocal G4Allocator<OCTREE>* OCTREE::fgAllocator = nullptr;
37//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39template <OCTREE_TEMPLATE>
41 : functor_(Extractor())
46//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48template <OCTREE_TEMPLATE>
49OCTREE::G4Octree(Iterator begin, Iterator end)
50 : G4Octree(begin,end,Extractor())
56//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58template <OCTREE_TEMPLATE>
59OCTREE::G4Octree(Iterator begin, Iterator end, Extractor f)
60 : functor_(f),head_(nullptr),size_(0)
62 std::vector<std::pair<Iterator,Point>> v;
63 for(auto it=begin;it!=end;++it)
65 v.push_back(std::pair<Iterator,Point>(it,functor_(it)));
71//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73template <OCTREE_TEMPLATE>
74OCTREE::G4Octree(OCTREE::tree_type&& rhs)
75 : functor_(rhs.functor_)
80//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82template <OCTREE_TEMPLATE>
83void OCTREE::swap(OCTREE::tree_type& rhs)
85 std::swap(head_, rhs.head_);
86 std::swap(functor_, rhs.functor_);
87 std::swap(size_, rhs.size_);
90//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92template <OCTREE_TEMPLATE>
93typename OCTREE::tree_type& OCTREE::operator=(typename OCTREE::tree_type rhs)
99//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101template <OCTREE_TEMPLATE>
102typename OCTREE::tree_type& OCTREE::operator=(typename OCTREE::tree_type&& rhs)
107//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109template <OCTREE_TEMPLATE>
115template <OCTREE_TEMPLATE>
116size_t OCTREE::size() const
121template <OCTREE_TEMPLATE>
122OCTREE::Node::Node(const NodeVector& input_values)
124 G4DNABoundingBox(InnerIterator(input_values.begin()),
125 InnerIterator(input_values.end())),
129//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
131template <OCTREE_TEMPLATE>
133 const NodeVector& input_values,
134 const G4DNABoundingBox& box,
135 size_t current_depth):
140 if (current_depth > max_depth)
142 init_max_depth_leaf(input_values);
144 else if (input_values.size() <= max_per_node)
146 init_leaf(input_values);
150 init_internal(input_values, current_depth);
154//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155template <OCTREE_TEMPLATE>
158 if (fNodeType == NodeTypes::INTERNAL)
160 childNodeArray& children = *static_cast<childNodeArray*>(fpValue);
161 for (size_t i = 0; i < 8; ++i)
163 if (children[i] != nullptr)
166 children[i] = nullptr;
171 else if (fNodeType == NodeTypes::LEAF)
173 auto toDelete = static_cast<LeafValues*>(fpValue);
175 delete static_cast<LeafValues*>(fpValue);
177 else if (fNodeType == NodeTypes::MAX_DEPTH_LEAF)
179 delete static_cast<NodeVector*>(fpValue);
184//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
186template <OCTREE_TEMPLATE>
187void OCTREE::Node::init_max_depth_leaf(
188 const NodeVector& input_values)
190 fpValue = new NodeVector(input_values);
191 fNodeType = NodeTypes::MAX_DEPTH_LEAF;
193//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
194template <OCTREE_TEMPLATE>
195void OCTREE::Node::init_leaf(const NodeVector& input_values)
197 std::array<std::pair<Iterator, Point>, max_per_node> a;
198 std::copy(input_values.begin(), input_values.end(), a.begin());
199 fpValue = new LeafValues{a, input_values.size()};
200 fNodeType = NodeTypes::LEAF;
202//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
203template <OCTREE_TEMPLATE>
204void OCTREE::Node::init_internal(
205 const NodeVector& input_values,
206 size_t current_depth)
208 std::array<NodeVector, 8> childVectors;
209 std::array<G4DNABoundingBox, 8> boxes = fBigVolume.partition();
210 std::array<Node*, 8> children;
211 for (size_t child = 0; child < 8; ++child)
213 NodeVector& childVector = childVectors[child];
214 childVector.reserve(input_values.size()/8);
215 std::copy_if(input_values.begin(),input_values.end(),
216 std::back_inserter(childVector),
217 [&boxes, child](const std::pair<Iterator, Point>& element)
220 const Point& p = element.second;
221 return boxes[child].contains(p);
224 children[child] = childVector.empty()
226 : new Node(childVector, boxes[child], ++current_depth);
228 fpValue = new std::array<Node*, 8>(children);
229 fNodeType = NodeTypes::INTERNAL;
231//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
233template <OCTREE_TEMPLATE>
234template <typename OutPutContainer>
235G4bool OCTREE::Node::radiusNeighbors(const Point& query,
237 OutPutContainer& resultIndices) const
239 G4bool success = false;
240 G4double distance = 0;
241 if (fNodeType == NodeTypes::INTERNAL)
243 childNodeArray& children = *static_cast<childNodeArray*>(fpValue);
244 for (auto eachChild : children)
246 if (eachChild == nullptr)
250 if(!eachChild->fBigVolume.overlap(query,radius))
254 success = eachChild->radiusNeighbors(query, radius, resultIndices) || success;
257 else if (fNodeType == NodeTypes::LEAF)
259 if(fpValue != nullptr)
261 LeafValues& children = *static_cast<LeafValues*>(fpValue);
262 for (size_t i = 0; i < children.size_; ++i)
264 distance = (query - std::get<1>(children.values_[i])).mag();
268 if( distance < radius )//TODO: find another solution for this using boundingbox
270 resultIndices.push_back(std::make_pair(std::get<0>(children.values_[i]),distance));
277 else if (fNodeType == NodeTypes::MAX_DEPTH_LEAF)
279 NodeVector& children = *static_cast<NodeVector*>(fpValue);
280 for (auto child : children)
282 const Point& point = std::get<1>(child);
283 //if (this->fBigVolume.contains(query, point, radius))
284 distance = (query - point).mag();
289 if( distance < radius )
293 throw std::runtime_error("distance == 0 => OCTREE::Node::radiusNeighbors : find itself");
295 Iterator resultIndex = std::get<0>(child);
296 resultIndices.push_back(std::make_pair(resultIndex,distance));
303 throw std::runtime_error("fNodeType is not set : find itself");
308//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
310template <OCTREE_TEMPLATE>
311template <typename OutPutContainer>
312void OCTREE::radiusNeighbors(const Point& query,
313 const G4double& radius,
314 OutPutContainer& resultIndices) const
316 resultIndices.clear();
317 if (head_ == nullptr)
321 head_->radiusNeighbors(query, radius, resultIndices);
324//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
326template <OCTREE_TEMPLATE>
327void* OCTREE::operator new(size_t)
331 fgAllocator = new G4Allocator<OCTREE>;
333 return (void *) fgAllocator->MallocSingle();
336template <OCTREE_TEMPLATE>
337void OCTREE::operator delete(void *a)
339 fgAllocator->FreeSingle((OCTREE*)a);