| 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 //
>> 27 // Author: Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr)
>> 28 //
>> 29 // History:
>> 30 // -----------
>> 31 // 10 Oct 2011 M.Karamitros created
>> 32 //
>> 33 // -------------------------------------------------------------------
26 /* 34 /*
27 * G4KDTree.cc << 35 * Based on ``kdtree'', a library for working with kd-trees.
>> 36 * Copyright (C) 2007-2009 John Tsiombikas <nuclear@siggraph.org>
>> 37 * The original open-source version of this code
>> 38 * may be found at http://code.google.com/p/kdtree/
28 * 39 *
29 * Created on: 22 oct. 2013 << 40 * Redistribution and use in source and binary forms, with or without
30 * Author: kara << 41 * modification, are permitted provided that the following conditions are
31 */ << 42 * met:
>> 43 * 1. Redistributions of source code must retain the above copyright
>> 44 * notice, this
>> 45 * list of conditions and the following disclaimer.
>> 46 * 2. Redistributions in binary form must reproduce the above copyright
>> 47 * notice,
>> 48 * this list of conditions and the following disclaimer in the
>> 49 * documentation
>> 50 * and/or other materials provided with the distribution.
>> 51 * 3. The name of the author may not be used to endorse or promote products
>> 52 * derived from this software without specific prior written permission.
>> 53 *
>> 54 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
>> 55 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
>> 56 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
>> 57 * IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
>> 58 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
>> 59 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
>> 60 * OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
>> 61 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
>> 62 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
>> 63 */
>> 64 /* single nearest neighbor search written by Tamas Nepusz
>> 65 * <tamas@cs.rhul.ac.uk>
>> 66 */
32 67
33 #include "globals.hh" 68 #include "globals.hh"
34 #include <cstdio> << 69 #include <stdio.h>
35 #include <cmath> << 70 #include <stdlib.h>
>> 71 #include <math.h>
36 #include "G4KDTree.hh" 72 #include "G4KDTree.hh"
37 #include "G4KDMap.hh" <<
38 #include "G4KDNode.hh" 73 #include "G4KDNode.hh"
39 #include "G4KDTreeResult.hh" 74 #include "G4KDTreeResult.hh"
40 #include <list> 75 #include <list>
41 #include <iostream> 76 #include <iostream>
42 77
43 using namespace std; 78 using namespace std;
44 79
45 G4Allocator<G4KDTree>*& G4KDTree::fgAllocator( << 80 //______________________________________________________________________
>> 81 struct HyperRect
46 { 82 {
47 G4ThreadLocalStatic G4Allocator<G4KDTree>* _ << 83 public:
48 return _instance; << 84 HyperRect(int dim, const double *min, const double *max)
49 } << 85 {
>> 86 fDim = dim;
>> 87 fMin = new double[fDim];
>> 88 fMax = new double[fDim];
>> 89 size_t size = fDim * sizeof(double);
>> 90 memcpy(fMin, min, size);
>> 91 memcpy(fMax, max, size);
>> 92 }
>> 93
>> 94
>> 95 ~HyperRect()
>> 96 {
>> 97 delete[] fMin;
>> 98 delete[] fMax;
>> 99 }
>> 100
>> 101 HyperRect(const HyperRect& rect)
>> 102 {
>> 103 fDim = rect.fDim;
>> 104 fMin = new double[fDim];
>> 105 fMax = new double[fDim];
>> 106 size_t size = fDim * sizeof(double);
>> 107 memcpy(fMin, rect.fMin, size);
>> 108 memcpy(fMax, rect.fMax, size);
>> 109 }
>> 110
>> 111 void Extend(const double *pos)
>> 112 {
>> 113 int i;
>> 114
>> 115 for (i=0; i < fDim; i++)
>> 116 {
>> 117 if (pos[i] < fMin[i])
>> 118 {
>> 119 fMin[i] = pos[i];
>> 120 }
>> 121 if (pos[i] > fMax[i])
>> 122 {
>> 123 fMax[i] = pos[i];
>> 124 }
>> 125 }
>> 126 }
>> 127
>> 128 bool CompareDistSqr(const double *pos, const double* bestmatch)
>> 129 {
>> 130 double result = 0;
>> 131
>> 132 for (int i=0; i < fDim; i++)
>> 133 {
>> 134 if (pos[i] < fMin[i])
>> 135 {
>> 136 result += sqr(fMin[i] - pos[i]);
>> 137 }
>> 138 else if (pos[i] > fMax[i])
>> 139 {
>> 140 result += sqr(fMax[i] - pos[i]);
>> 141 }
>> 142
>> 143 if(result >= *bestmatch) return false ;
>> 144 }
>> 145
>> 146 return true ;
>> 147 }
>> 148
>> 149 int GetDim(){return fDim;}
>> 150 double* GetMin(){return fMin;}
>> 151 double* GetMax(){return fMax;}
>> 152
>> 153 protected:
>> 154 int fDim;
>> 155 double *fMin, *fMax; /* minimum/maximum coords */
>> 156
>> 157 private:
>> 158 // should not be used
>> 159 HyperRect& operator=(const HyperRect& rhs)
>> 160 {
>> 161 if(this == &rhs) return *this;
>> 162 return *this;
>> 163 }
>> 164 };
50 165
51 //____________________________________________ 166 //______________________________________________________________________
52 // KDTree methods 167 // KDTree methods
53 G4KDTree::G4KDTree(size_t k) << 168 G4KDTree::G4KDTree (int k)
54 : fDim(k) << 169 {
55 ,fKDMap(new G4KDMap(k)) << 170 fDim = k;
56 {} << 171 fRoot = 0;
>> 172 fDestr = 0;
>> 173 fRect = 0;
>> 174 fNbNodes = 0;
>> 175 }
57 176
58 G4KDTree::~G4KDTree() << 177 G4KDTree::~G4KDTree ()
>> 178 {
>> 179 if(fRoot) __Clear_Rec(fRoot);
>> 180 fRoot = 0;
>> 181
>> 182 if (fRect)
>> 183 {
>> 184 delete fRect;
>> 185 fRect = 0;
>> 186 }
>> 187 }
>> 188
>> 189 void G4KDTree::Clear()
59 { 190 {
60 if(fRoot != nullptr){ <<
61 __Clear_Rec(fRoot); 191 __Clear_Rec(fRoot);
62 fRoot = nullptr; << 192 fRoot = 0;
63 } << 193 fNbNodes = 0;
>> 194
>> 195 if (fRect)
>> 196 {
>> 197 delete fRect;
>> 198 fRect = 0;
>> 199 }
>> 200 }
64 201
65 if(fRect != nullptr){ << 202 void G4KDTree::__Clear_Rec(G4KDNode *node)
66 delete fRect; << 203 {
67 fRect = nullptr; << 204 if(!node) return;
68 } <<
69 205
70 if(fKDMap != nullptr){ << 206 if(node->GetLeft()) __Clear_Rec(node->GetLeft());
71 delete fKDMap; << 207 if(node->GetRight()) __Clear_Rec(node->GetRight());
72 fKDMap = nullptr; << 208
73 } << 209 if(fDestr)
>> 210 {
>> 211 if(node->GetData())
>> 212 {
>> 213 fDestr(node->GetData());
>> 214 node->SetData(0);
>> 215 }
>> 216 }
>> 217 delete node;
74 } 218 }
75 219
76 void* G4KDTree::operator new(size_t) << 220 G4KDNode* G4KDTree::Insert(const double *pos, void *data)
77 { 221 {
78 if(fgAllocator() == nullptr){ << 222 G4KDNode* node = 0 ;
79 fgAllocator() = new G4Allocator<G4KDTree>; << 223 if(!fRoot)
80 } << 224 {
81 return (void*) fgAllocator()->MallocSingle() << 225 fRoot = new G4KDNode(this,pos,data,0, 0);
>> 226 node = fRoot;
>> 227 fNbNodes = 0;
>> 228 fNbNodes++;
>> 229 }
>> 230 else
>> 231 {
>> 232 if((node=fRoot->Insert(pos, data)))
>> 233 {
>> 234 fNbNodes++;
>> 235 }
>> 236 }
>> 237
>> 238 if (fRect == 0)
>> 239 {
>> 240 fRect = new HyperRect(fDim,pos,pos);
>> 241 }
>> 242 else
>> 243 {
>> 244 fRect->Extend(pos);
>> 245 }
>> 246
>> 247 return node;
82 } 248 }
83 249
84 void G4KDTree::operator delete(void* aNode) << 250 G4KDNode* G4KDTree::Insert(const double& x, const double& y, const double& z, void *data)
85 { 251 {
86 fgAllocator()->FreeSingle((G4KDTree*) aNode) << 252 double buf[3];
>> 253 buf[0] = x;
>> 254 buf[1] = y;
>> 255 buf[2] = z;
>> 256 return Insert(buf, data);
87 } 257 }
88 258
89 void G4KDTree::Print(std::ostream& out) const << 259 //__________________________________________________________________
>> 260 int G4KDTree::__NearestInRange(G4KDNode *node, const double *pos, const double& range_sq,
>> 261 const double& range, G4KDTreeResult& list, int ordered, G4KDNode *source_node)
90 { 262 {
91 if(fRoot != nullptr){ << 263 if(!node) return 0;
92 fRoot->Print(out); << 264
93 } << 265 double dist_sq(DBL_MAX), dx(DBL_MAX);
>> 266 int ret(-1), added_res(0);
>> 267
>> 268 if(node-> GetData() && node != source_node)
>> 269 {
>> 270 bool do_break = false ;
>> 271 dist_sq = 0;
>> 272 for(int i=0; i<fDim; i++)
>> 273 {
>> 274 dist_sq += sqr(node->GetPosition()[i] - pos[i]);
>> 275 if(dist_sq > range_sq)
>> 276 {
>> 277 do_break = true;
>> 278 break;
>> 279 }
>> 280 }
>> 281 if(!do_break && dist_sq <= range_sq)
>> 282 {
>> 283 list.Insert(dist_sq, node);
>> 284 added_res = 1;
>> 285 }
>> 286 }
>> 287
>> 288 dx = pos[node->GetAxis()] - node->GetPosition()[node->GetAxis()];
>> 289
>> 290 ret = __NearestInRange(dx <= 0.0 ? node->GetLeft() : node->GetRight(), pos, range_sq, range, list, ordered, source_node);
>> 291 if(ret >= 0 && fabs(dx) <= range)
>> 292 {
>> 293 added_res += ret;
>> 294 ret = __NearestInRange(dx <= 0.0 ? node->GetRight() : node->GetLeft(), pos, range_sq, range, list, ordered, source_node);
>> 295 }
>> 296
>> 297 if(ret == -1)
>> 298 {
>> 299 return -1;
>> 300 }
>> 301 added_res += ret;
>> 302
>> 303 return added_res;
94 } 304 }
95 305
96 void G4KDTree::Clear() << 306 //__________________________________________________________________
>> 307 void G4KDTree::__NearestToPosition(G4KDNode *node, const double *pos, G4KDNode *&result,
>> 308 double *result_dist_sq, HyperRect* rect)
97 { 309 {
98 __Clear_Rec(fRoot); << 310 int dir = node->GetAxis();
99 fRoot = nullptr; << 311 int i;
100 fNbNodes = 0; << 312 double dummy(0.), dist_sq(-1.);
>> 313 G4KDNode *nearer_subtree(0), *farther_subtree (0);
>> 314 double *nearer_hyperrect_coord(0),*farther_hyperrect_coord(0);
>> 315
>> 316 /* Decide whether to go left or right in the tree */
>> 317 dummy = pos[dir] - node->GetPosition()[dir];
>> 318 if (dummy <= 0)
>> 319 {
>> 320 nearer_subtree = node->GetLeft();
>> 321 farther_subtree = node->GetRight();
101 322
102 if(fRect != nullptr) << 323 nearer_hyperrect_coord = rect->GetMax() + dir;
103 { << 324 farther_hyperrect_coord = rect->GetMin() + dir;
104 delete fRect; << 325 }
105 fRect = nullptr; << 326 else
106 } << 327 {
>> 328 nearer_subtree = node->GetRight();
>> 329 farther_subtree = node->GetLeft();
>> 330 nearer_hyperrect_coord = rect->GetMin() + dir;
>> 331 farther_hyperrect_coord = rect->GetMax() + dir;
>> 332 }
>> 333
>> 334 if (nearer_subtree)
>> 335 {
>> 336 /* Slice the hyperrect to get the hyperrect of the nearer subtree */
>> 337 dummy = *nearer_hyperrect_coord;
>> 338 *nearer_hyperrect_coord = node->GetPosition()[dir];
>> 339 /* Recurse down into nearer subtree */
>> 340 __NearestToPosition(nearer_subtree, pos, result, result_dist_sq, rect);
>> 341 /* Undo the slice */
>> 342 *nearer_hyperrect_coord = dummy;
>> 343 }
>> 344
>> 345 /* Check the distance of the point at the current node, compare it
>> 346 * with our best so far */
>> 347 if(node->GetData())
>> 348 {
>> 349 dist_sq = 0;
>> 350 bool do_break = false ;
>> 351 for(i=0; i < fDim; i++)
>> 352 {
>> 353 dist_sq += sqr(node->GetPosition()[i] - pos[i]);
>> 354 if(dist_sq > *result_dist_sq)
>> 355 {
>> 356 do_break = true;
>> 357 break ;
>> 358 }
>> 359 }
>> 360 if (!do_break && dist_sq < *result_dist_sq)
>> 361 {
>> 362 result = node;
>> 363 *result_dist_sq = dist_sq;
>> 364 }
>> 365 }
>> 366
>> 367 if (farther_subtree)
>> 368 {
>> 369 /* Get the hyperrect of the farther subtree */
>> 370 dummy = *farther_hyperrect_coord;
>> 371 *farther_hyperrect_coord = node->GetPosition()[dir];
>> 372 /* Check if we have to recurse down by calculating the closest
>> 373 * point of the hyperrect and see if it's closer than our
>> 374 * minimum distance in result_dist_sq. */
>> 375 if (rect->CompareDistSqr(pos,result_dist_sq))
>> 376 {
>> 377 /* Recurse down into farther subtree */
>> 378 __NearestToPosition(farther_subtree, pos, result, result_dist_sq, rect);
>> 379 }
>> 380 /* Undo the slice on the hyperrect */
>> 381 *farther_hyperrect_coord = dummy;
>> 382 }
107 } 383 }
108 384
109 void G4KDTree::__Clear_Rec(G4KDNode_Base* node << 385 G4KDTreeResultHandle G4KDTree::Nearest(const double *pos)
110 { 386 {
111 if(node == nullptr) << 387 // G4cout << "Nearest(pos)" << G4endl ;
112 { << 388
113 return; << 389 if (!fRect) return 0;
114 } << 390
>> 391 G4KDNode *result(0);
>> 392 double dist_sq = DBL_MAX;
>> 393
>> 394 /* Duplicate the bounding hyperrectangle, we will work on the copy */
>> 395 HyperRect *newrect = new HyperRect(*fRect);
>> 396
>> 397 /* Our first guesstimate is the root node */
>> 398 /* Search for the nearest neighbour recursively */
>> 399 __NearestToPosition(fRoot, pos, result, &dist_sq, newrect);
115 400
116 if(node->GetLeft() != nullptr) << 401 /* Free the copy of the hyperrect */
117 { << 402 delete newrect;
118 __Clear_Rec(node->GetLeft()); <<
119 } <<
120 if(node->GetRight() != nullptr) <<
121 { <<
122 __Clear_Rec(node->GetRight()); <<
123 } <<
124 403
125 delete node; << 404 /* Store the result */
>> 405 if (result)
>> 406 {
>> 407 G4KDTreeResultHandle rset = new G4KDTreeResult(this);
>> 408 rset->Insert(dist_sq, result);
>> 409 rset -> Rewind();
>> 410 return rset;
>> 411 }
>> 412 else
>> 413 {
>> 414 return 0;
>> 415 }
126 } 416 }
127 417
128 void G4KDTree::__InsertMap(G4KDNode_Base* node << 418 //__________________________________________________________________
>> 419 void G4KDTree::__NearestToNode(G4KDNode *source_node, G4KDNode *node,
>> 420 const double *pos, std::vector<G4KDNode*>& result, double *result_dist_sq,
>> 421 HyperRect* rect, int& nbresult)
>> 422 {
>> 423 int dir = node->GetAxis();
>> 424 double dummy, dist_sq;
>> 425 G4KDNode *nearer_subtree (0), *farther_subtree (0);
>> 426 double *nearer_hyperrect_coord (0), *farther_hyperrect_coord (0);
>> 427
>> 428 /* Decide whether to go left or right in the tree */
>> 429 dummy = pos[dir] - node->GetPosition()[dir];
>> 430 if (dummy <= 0)
>> 431 {
>> 432 nearer_subtree = node->GetLeft();
>> 433 farther_subtree = node->GetRight();
>> 434 nearer_hyperrect_coord = rect->GetMax() + dir;
>> 435 farther_hyperrect_coord = rect->GetMin() + dir;
>> 436 }
>> 437 else
>> 438 {
>> 439 nearer_subtree = node->GetRight();
>> 440 farther_subtree = node->GetLeft();
>> 441 nearer_hyperrect_coord = rect->GetMin() + dir;
>> 442 farther_hyperrect_coord = rect->GetMax() + dir;
>> 443 }
>> 444
>> 445 if (nearer_subtree)
>> 446 {
>> 447 /* Slice the hyperrect to get the hyperrect of the nearer subtree */
>> 448 dummy = *nearer_hyperrect_coord;
>> 449 *nearer_hyperrect_coord = node->GetPosition()[dir];
>> 450 /* Recurse down into nearer subtree */
>> 451 __NearestToNode(source_node, nearer_subtree, pos, result, result_dist_sq, rect, nbresult);
>> 452 /* Undo the slice */
>> 453 *nearer_hyperrect_coord = dummy;
>> 454 }
>> 455
>> 456 /* Check the distance of the point at the current node, compare it
>> 457 * with our best so far */
>> 458 if(node->GetData() && node != source_node)
>> 459 {
>> 460 dist_sq = 0;
>> 461 bool do_break = false;
>> 462 for(int i=0; i < fDim; i++)
>> 463 {
>> 464 dist_sq += sqr(node->GetPosition()[i] - pos[i]);
>> 465 if(dist_sq > *result_dist_sq)
>> 466 {
>> 467 do_break = true;
>> 468 break ;
>> 469 }
>> 470 }
>> 471 if(!do_break)
>> 472 {
>> 473 if (dist_sq < *result_dist_sq)
>> 474 {
>> 475 result.clear();
>> 476 nbresult = 1 ;
>> 477 result.push_back(node);
>> 478 *result_dist_sq = dist_sq;
>> 479 }
>> 480 else if(dist_sq == *result_dist_sq)
>> 481 {
>> 482 result.push_back(node);
>> 483 nbresult++;
>> 484 }
>> 485 }
>> 486 }
129 487
130 void G4KDTree::Build() << 488 if (farther_subtree)
>> 489 {
>> 490 /* Get the hyperrect of the farther subtree */
>> 491 dummy = *farther_hyperrect_coord;
>> 492 *farther_hyperrect_coord = node->GetPosition()[dir];
>> 493 /* Check if we have to recurse down by calculating the closest
>> 494 * point of the hyperrect and see if it's closer than our
>> 495 * minimum distance in result_dist_sq. */
>> 496 // if (hyperrect_dist_sq(rect, pos) < *result_dist_sq)
>> 497 if (rect->CompareDistSqr(pos,result_dist_sq))
>> 498 {
>> 499 /* Recurse down into farther subtree */
>> 500 __NearestToNode(source_node, farther_subtree, pos, result, result_dist_sq, rect, nbresult);
>> 501 }
>> 502 /* Undo the slice on the hyperrect */
>> 503 *farther_hyperrect_coord = dummy;
>> 504 }
>> 505 }
>> 506
>> 507 G4KDTreeResultHandle G4KDTree::Nearest(G4KDNode* node)
131 { 508 {
132 size_t Nnodes = fKDMap->GetSize(); << 509 // G4cout << "Nearest(node)" << G4endl ;
>> 510 if (!fRect)
>> 511 {
>> 512 G4cout << "Tree empty" << G4endl ;
>> 513 return 0;
>> 514 }
133 515
134 G4cout << "********************" << G4endl; << 516 const double* pos = node->GetPosition();
135 G4cout << "template<typename PointT> G4KDTre << 517 std::vector<G4KDNode*> result;
136 G4cout << "Map size = " << Nnodes << G4endl; << 518 double dist_sq = DBL_MAX;
137 519
138 G4KDNode_Base* root = fKDMap->PopOutMiddle(0 << 520 /* Duplicate the bounding hyperrectangle, we will work on the copy */
>> 521 HyperRect *newrect = new HyperRect(*fRect);
139 522
140 if(root == nullptr) << 523 /* Search for the nearest neighbour recursively */
141 { << 524 int nbresult = 0 ;
142 return; <<
143 } <<
144 525
145 fRoot = root; << 526 __NearestToNode(node, fRoot, pos, result, &dist_sq, newrect, nbresult);
146 fNbActiveNodes++; <<
147 fRect = new HyperRect(fDim); <<
148 fRect->SetMinMax(*fRoot, *fRoot); <<
149 527
150 Nnodes--; << 528 /* Free the copy of the hyperrect */
>> 529 delete newrect;
151 530
152 G4KDNode_Base* parent = fRoot; << 531 /* Store the result */
>> 532 if (!result.empty())
>> 533 {
>> 534 G4KDTreeResultHandle rset(new G4KDTreeResult(this));
>> 535 int j = 0 ;
>> 536 while (j<nbresult)
>> 537 {
>> 538 rset->Insert(dist_sq, result[j]);
>> 539 j++;
>> 540 }
>> 541 rset->Rewind();
153 542
154 for(size_t n = 0; n < Nnodes; n += fDim) << 543 return rset;
155 { << 544 }
156 for(size_t dim = 0; dim < fDim; dim++) << 545 else
157 { 546 {
158 G4KDNode_Base* node = fKDMap->PopOutMidd << 547 return 0;
159 if(node != nullptr) <<
160 { <<
161 parent->Insert(node); <<
162 fNbActiveNodes++; <<
163 fRect->Extend(*node); <<
164 parent = node; <<
165 } <<
166 } 548 }
167 } <<
168 } 549 }
169 550
170 G4KDTreeResultHandle G4KDTree::Nearest(G4KDNod << 551 G4KDTreeResultHandle G4KDTree::Nearest(const double& x, const double& y, const double& z)
>> 552 {
>> 553 double pos[3];
>> 554 pos[0] = x;
>> 555 pos[1] = y;
>> 556 pos[2] = z;
>> 557 return Nearest(pos);
>> 558 }
>> 559
>> 560 G4KDTreeResultHandle G4KDTree::NearestInRange(const double *pos, const double& range)
171 { 561 {
172 if(fRect == nullptr) << 562 int ret(-1);
173 { <<
174 return nullptr; <<
175 } <<
176 563
177 std::vector<G4KDNode_Base*> result; << 564 const double range_sq = sqr(range) ;
178 G4double dist_sq = DBL_MAX; << 565
>> 566 G4KDTreeResultHandle rset = new G4KDTreeResult(this);
>> 567 if((ret = __NearestInRange(fRoot, pos, range_sq, range, *(rset()), 0)) == -1)
>> 568 {
>> 569 rset = 0;
>> 570 return rset;
>> 571 }
>> 572 rset->Sort();
>> 573 rset->Rewind();
>> 574 return rset;
>> 575 }
179 576
180 /* Duplicate the bounding hyperrectangle, we << 577 G4KDTreeResultHandle G4KDTree::NearestInRange(const double& x,
181 auto newrect = new HyperRect(*fRect); << 578 const double& y,
>> 579 const double& z,
>> 580 const double& range)
>> 581 {
>> 582 double buf[3];
>> 583 buf[0] = x;
>> 584 buf[1] = y;
>> 585 buf[2] = z;
>> 586 return NearestInRange(buf, range);
>> 587 }
182 588
183 /* Search for the nearest neighbour recursiv << 589 G4KDTreeResultHandle G4KDTree::NearestInRange( G4KDNode* node, const double& range)
184 G4int nbresult = 0; << 590 {
>> 591 if(!node) return 0 ;
>> 592 int ret(-1);
185 593
186 __NearestToNode(node, fRoot, *node, result, << 594 G4KDTreeResult *rset = new G4KDTreeResult(this);
187 595
188 /* Free the copy of the hyperrect */ << 596 const double range_sq = sqr(range) ;
189 delete newrect; <<
190 597
191 /* Store the result */ << 598 if((ret = __NearestInRange(fRoot, node->GetPosition(), range_sq, range, *rset, 0, node)) == -1)
192 if(!result.empty()) <<
193 { <<
194 G4KDTreeResultHandle rset(new G4KDTreeResu <<
195 G4int j = 0; <<
196 while(j < nbresult) <<
197 { 599 {
198 rset->Insert(dist_sq, result[j]); << 600 delete rset;
199 j++; << 601 return 0;
200 } 602 }
>> 603 rset->Sort();
201 rset->Rewind(); 604 rset->Rewind();
202 <<
203 return rset; 605 return rset;
204 } <<
205 <<
206 return nullptr; <<
207 } <<
208 <<
209 G4KDTreeResultHandle G4KDTree::NearestInRange( <<
210 <<
211 { <<
212 if(node == nullptr) <<
213 { <<
214 return nullptr; <<
215 } <<
216 G4int ret(-1); <<
217 <<
218 auto* rset = new G4KDTreeResult(this); <<
219 <<
220 const G4double range_sq = sqr(range); <<
221 <<
222 if((ret = __NearestInRange(fRoot, *node, ran <<
223 -1) <<
224 { <<
225 delete rset; <<
226 return nullptr; <<
227 } <<
228 rset->Sort(); <<
229 rset->Rewind(); <<
230 return rset; <<
231 } 606 }
232 607