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