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 /* 26 /* 27 * G4KDTree.cc 27 * G4KDTree.cc 28 * 28 * 29 * Created on: 22 oct. 2013 29 * Created on: 22 oct. 2013 30 * Author: kara 30 * Author: kara 31 */ 31 */ 32 32 33 #include "globals.hh" 33 #include "globals.hh" 34 #include <cstdio> << 34 #include <stdio.h> >> 35 #include <stdlib.h> 35 #include <cmath> 36 #include <cmath> 36 #include "G4KDTree.hh" 37 #include "G4KDTree.hh" 37 #include "G4KDMap.hh" 38 #include "G4KDMap.hh" 38 #include "G4KDNode.hh" 39 #include "G4KDNode.hh" 39 #include "G4KDTreeResult.hh" 40 #include "G4KDTreeResult.hh" 40 #include <list> 41 #include <list> 41 #include <iostream> 42 #include <iostream> 42 43 43 using namespace std; 44 using namespace std; 44 45 45 G4Allocator<G4KDTree>*& G4KDTree::fgAllocator( 46 G4Allocator<G4KDTree>*& G4KDTree::fgAllocator() 46 { 47 { 47 G4ThreadLocalStatic G4Allocator<G4KDTree>* _ << 48 G4ThreadLocalStatic G4Allocator<G4KDTree>* _instance = nullptr; 48 return _instance; << 49 return _instance; 49 } 50 } 50 51 51 //____________________________________________ 52 //______________________________________________________________________ 52 // KDTree methods 53 // KDTree methods 53 G4KDTree::G4KDTree(size_t k) << 54 G4KDTree::G4KDTree(size_t k) : 54 : fDim(k) << 55 fKDMap(new G4KDMap(k)) 55 ,fKDMap(new G4KDMap(k)) << 56 { 56 {} << 57 fDim = k; >> 58 fRoot = 0; >> 59 fRect = 0; >> 60 fNbNodes = 0; >> 61 fNbActiveNodes = 0; >> 62 } 57 63 58 G4KDTree::~G4KDTree() 64 G4KDTree::~G4KDTree() 59 { 65 { 60 if(fRoot != nullptr){ << 66 if (fRoot) __Clear_Rec(fRoot); 61 __Clear_Rec(fRoot); << 67 fRoot = 0; 62 fRoot = nullptr; << 63 } << 64 68 65 if(fRect != nullptr){ << 69 if (fRect) >> 70 { 66 delete fRect; 71 delete fRect; 67 fRect = nullptr; << 72 fRect = 0; 68 } 73 } 69 74 70 if(fKDMap != nullptr){ << 75 if (fKDMap) delete fKDMap; 71 delete fKDMap; << 72 fKDMap = nullptr; << 73 } << 74 } 76 } 75 77 76 void* G4KDTree::operator new(size_t) 78 void* G4KDTree::operator new(size_t) 77 { 79 { 78 if(fgAllocator() == nullptr){ << 80 if (!fgAllocator()) fgAllocator() = new G4Allocator<G4KDTree>; 79 fgAllocator() = new G4Allocator<G4KDTree>; << 81 return (void *) fgAllocator()->MallocSingle(); 80 } << 81 return (void*) fgAllocator()->MallocSingle() << 82 } 82 } 83 83 84 void G4KDTree::operator delete(void* aNode) << 84 void G4KDTree::operator delete(void *aNode) 85 { 85 { 86 fgAllocator()->FreeSingle((G4KDTree*) aNode) 86 fgAllocator()->FreeSingle((G4KDTree*) aNode); 87 } 87 } 88 88 89 void G4KDTree::Print(std::ostream& out) const 89 void G4KDTree::Print(std::ostream& out) const 90 { 90 { 91 if(fRoot != nullptr){ << 91 if (fRoot) fRoot->Print(out); 92 fRoot->Print(out); << 93 } << 94 } 92 } 95 93 96 void G4KDTree::Clear() 94 void G4KDTree::Clear() 97 { 95 { 98 __Clear_Rec(fRoot); 96 __Clear_Rec(fRoot); 99 fRoot = nullptr; << 97 fRoot = 0; 100 fNbNodes = 0; 98 fNbNodes = 0; 101 99 102 if(fRect != nullptr) << 100 if (fRect) 103 { 101 { 104 delete fRect; 102 delete fRect; 105 fRect = nullptr; << 103 fRect = 0; 106 } 104 } 107 } 105 } 108 106 109 void G4KDTree::__Clear_Rec(G4KDNode_Base* node 107 void G4KDTree::__Clear_Rec(G4KDNode_Base* node) 110 { 108 { 111 if(node == nullptr) << 109 if (!node) return; 112 { << 113 return; << 114 } << 115 110 116 if(node->GetLeft() != nullptr) << 111 if (node->GetLeft()) __Clear_Rec(node->GetLeft()); 117 { << 112 if (node->GetRight()) __Clear_Rec(node->GetRight()); 118 __Clear_Rec(node->GetLeft()); << 119 } << 120 if(node->GetRight() != nullptr) << 121 { << 122 __Clear_Rec(node->GetRight()); << 123 } << 124 113 125 delete node; 114 delete node; 126 } 115 } 127 116 128 void G4KDTree::__InsertMap(G4KDNode_Base* node << 117 void G4KDTree::__InsertMap(G4KDNode_Base *node) >> 118 { >> 119 fKDMap->Insert(node); >> 120 } 129 121 130 void G4KDTree::Build() 122 void G4KDTree::Build() 131 { 123 { 132 size_t Nnodes = fKDMap->GetSize(); 124 size_t Nnodes = fKDMap->GetSize(); 133 125 134 G4cout << "********************" << G4endl; 126 G4cout << "********************" << G4endl; 135 G4cout << "template<typename PointT> G4KDTre 127 G4cout << "template<typename PointT> G4KDTree<PointT>::Build" << G4endl; 136 G4cout << "Map size = " << Nnodes << G4endl; 128 G4cout << "Map size = " << Nnodes << G4endl; 137 129 138 G4KDNode_Base* root = fKDMap->PopOutMiddle(0 130 G4KDNode_Base* root = fKDMap->PopOutMiddle(0); 139 131 140 if(root == nullptr) << 132 if(root == 0) return; 141 { << 142 return; << 143 } << 144 133 145 fRoot = root; 134 fRoot = root; 146 fNbActiveNodes++; 135 fNbActiveNodes++; 147 fRect = new HyperRect(fDim); 136 fRect = new HyperRect(fDim); 148 fRect->SetMinMax(*fRoot, *fRoot); 137 fRect->SetMinMax(*fRoot, *fRoot); 149 138 150 Nnodes--; 139 Nnodes--; 151 140 152 G4KDNode_Base* parent = fRoot; 141 G4KDNode_Base* parent = fRoot; 153 142 154 for(size_t n = 0; n < Nnodes; n += fDim) << 143 for (size_t n = 0; n < Nnodes; n += fDim) 155 { 144 { 156 for(size_t dim = 0; dim < fDim; dim++) << 145 for (size_t dim = 0; dim < fDim; dim++) 157 { 146 { 158 G4KDNode_Base* node = fKDMap->PopOutMidd 147 G4KDNode_Base* node = fKDMap->PopOutMiddle(dim); 159 if(node != nullptr) << 148 if (node) 160 { 149 { 161 parent->Insert(node); 150 parent->Insert(node); 162 fNbActiveNodes++; 151 fNbActiveNodes++; 163 fRect->Extend(*node); 152 fRect->Extend(*node); 164 parent = node; 153 parent = node; 165 } 154 } 166 } 155 } 167 } 156 } 168 } 157 } 169 158 170 G4KDTreeResultHandle G4KDTree::Nearest(G4KDNod 159 G4KDTreeResultHandle G4KDTree::Nearest(G4KDNode_Base* node) 171 { 160 { 172 if(fRect == nullptr) << 161 // G4cout << "Nearest(node)" << G4endl ; >> 162 if (!fRect) 173 { 163 { 174 return nullptr; << 164 G4cout << "Tree empty" << G4endl; >> 165 return 0; 175 } 166 } 176 167 177 std::vector<G4KDNode_Base*> result; << 168 std::vector<G4KDNode_Base* > result; 178 G4double dist_sq = DBL_MAX; << 169 double dist_sq = DBL_MAX; 179 170 180 /* Duplicate the bounding hyperrectangle, we 171 /* Duplicate the bounding hyperrectangle, we will work on the copy */ 181 auto newrect = new HyperRect(*fRect); << 172 HyperRect *newrect = new HyperRect(*fRect); 182 173 183 /* Search for the nearest neighbour recursiv 174 /* Search for the nearest neighbour recursively */ 184 G4int nbresult = 0; << 175 int nbresult = 0; 185 176 186 __NearestToNode(node, fRoot, *node, result, 177 __NearestToNode(node, fRoot, *node, result, &dist_sq, newrect, nbresult); 187 178 188 /* Free the copy of the hyperrect */ 179 /* Free the copy of the hyperrect */ 189 delete newrect; 180 delete newrect; 190 181 191 /* Store the result */ 182 /* Store the result */ 192 if(!result.empty()) << 183 if (!result.empty()) 193 { 184 { 194 G4KDTreeResultHandle rset(new G4KDTreeResu 185 G4KDTreeResultHandle rset(new G4KDTreeResult(this)); 195 G4int j = 0; << 186 int j = 0; 196 while(j < nbresult) << 187 while (j<nbresult) 197 { 188 { 198 rset->Insert(dist_sq, result[j]); 189 rset->Insert(dist_sq, result[j]); 199 j++; 190 j++; 200 } 191 } 201 rset->Rewind(); 192 rset->Rewind(); 202 193 203 return rset; 194 return rset; 204 } 195 } 205 << 196 else 206 return nullptr; << 197 { >> 198 return 0; >> 199 } 207 } 200 } 208 201 209 G4KDTreeResultHandle G4KDTree::NearestInRange( 202 G4KDTreeResultHandle G4KDTree::NearestInRange(G4KDNode_Base* node, 210 << 203 const double& range) 211 { 204 { 212 if(node == nullptr) << 205 if (!node) return 0; 213 { << 206 int ret(-1); 214 return nullptr; << 215 } << 216 G4int ret(-1); << 217 207 218 auto* rset = new G4KDTreeResult(this); << 208 G4KDTreeResult *rset = new G4KDTreeResult(this); 219 209 220 const G4double range_sq = sqr(range); << 210 const double range_sq = sqr(range); 221 211 222 if((ret = __NearestInRange(fRoot, *node, ran << 212 if ((ret = __NearestInRange(fRoot, *node, range_sq, range, *rset, 0, node)) == -1) 223 -1) << 224 { 213 { 225 delete rset; 214 delete rset; 226 return nullptr; << 215 return 0; 227 } 216 } 228 rset->Sort(); 217 rset->Sort(); 229 rset->Rewind(); 218 rset->Rewind(); 230 return rset; 219 return rset; 231 } 220 } 232 221