Geant4-11
G4KDTree.icc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
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. *
10// * *
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. *
17// * *
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// ********************************************************************
25//
26//
27// Author: Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr)
28//
29// History:
30// -----------
31// 10 Oct 2011 M.Karamitros created
32//
33// -------------------------------------------------------------------
34/*
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/
39 *
40 * Redistribution and use in source and binary forms, with or without
41 * modification, are permitted provided that the following conditions are
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 */
67
68template<typename PointT>
69 G4KDNode_Base* G4KDTree::InsertMap(PointT* point)
70 {
71 G4KDNode<PointT>* node = new G4KDNode<PointT>(this, point, 0);
72 this->__InsertMap(node);
73 return node;
74 }
75
76template<typename PointT>
77 G4KDNode_Base* G4KDTree::Insert(PointT* pos)
78 {
79 G4KDNode_Base* node = 0;
80 if (!fRoot)
81 {
82 fRoot = new G4KDNode<PointT>(this, pos, 0);
83 node = fRoot;
84 fNbNodes = 0;
85 fNbNodes++;
86 fNbActiveNodes++;
87 }
88 else
89 {
90 if ((node = fRoot->Insert<PointT>(pos)))
91 {
92 fNbNodes++;
93 fNbActiveNodes++;
94 }
95 }
96
97 if (fRect == 0)
98 {
99 fRect = new HyperRect(fDim);
100 fRect->SetMinMax(*pos, *pos);
101 }
102 else
103 {
104 fRect->Extend(*pos);
105 }
106
107 return node;
108 }
109
110template<typename PointT>
111 G4KDNode_Base* G4KDTree::Insert(const PointT& pos)
112 {
113 G4KDNode_Base* node = 0;
114 if (!fRoot)
115 {
116 fRoot = new G4KDNodeCopy<PointT>(this, pos, 0);
117 node = fRoot;
118 fNbNodes = 0;
119 fNbNodes++;
120 fNbActiveNodes++;
121 }
122 else
123 {
124 if ((node = fRoot->Insert<PointT>(pos)))
125 {
126 fNbNodes++;
127 fNbActiveNodes++;
128 }
129 }
130
131 if (fRect == 0)
132 {
133 fRect = new HyperRect(fDim);
134 fRect->SetMinMax(pos, pos);
135 }
136 else
137 {
138 fRect->Extend(pos);
139 }
140
141 return node;
142 }
143
144//__________________________________________________________________
145template<typename Position>
146 int G4KDTree::__NearestInRange(G4KDNode_Base* node,
147 const Position& pos,
148 const double& range_sq,
149 const double& range,
150 G4KDTreeResult& list,
151 int ordered,
152 G4KDNode_Base *source_node)
153 {
154 if (!node) return 0;
155
156 double dist_sq(DBL_MAX), dx(DBL_MAX);
157 int ret(-1), added_res(0);
158
159 if (node->IsValid() && node != source_node)
160 {
161 bool do_break = false;
162 dist_sq = 0;
163 for (size_t i = 0; i < fDim; i++)
164 {
165 dist_sq += sqr((*node)[i] - pos[i]);
166 if (dist_sq > range_sq)
167 {
168 do_break = true;
169 break;
170 }
171 }
172 if (!do_break && dist_sq <= range_sq)
173 {
174 list.Insert(dist_sq, node);
175 added_res = 1;
176 }
177 }
178
179 dx = pos[node->GetAxis()] - (*node)[node->GetAxis()];
180
181 ret = __NearestInRange(dx <= 0.0 ? node->GetLeft() : node->GetRight(), pos,
182 range_sq, range, list, ordered, source_node);
183 if (ret >= 0 && std::fabs(dx) <= range)
184 {
185 added_res += ret;
186 ret = __NearestInRange(dx <= 0.0 ? node->GetRight() : node->GetLeft(),
187 pos, range_sq, range, list, ordered, source_node);
188 }
189
190 if (ret == -1)
191 {
192 return -1;
193 }
194 added_res += ret;
195
196 return added_res;
197 }
198
199//__________________________________________________________________
200template<typename Position>
201 void G4KDTree::__NearestToPosition(G4KDNode_Base *node,
202 const Position &pos,
203 G4KDNode_Base *&result,
204 double *result_dist_sq,
205 HyperRect* rect)
206 {
207 int dir = node->GetAxis();
208 double dummy(0.), dist_sq(-1.);
209 G4KDNode_Base* nearer_subtree(0), *farther_subtree(0);
210 double *nearer_hyperrect_coord(0), *farther_hyperrect_coord(0);
211
212 /* Decide whether to go left or right in the tree */
213 dummy = pos[dir] - (*node)[dir];
214 if (dummy <= 0)
215 {
216 nearer_subtree = node->GetLeft();
217 farther_subtree = node->GetRight();
218
219 nearer_hyperrect_coord = rect->GetMax() + dir;
220 farther_hyperrect_coord = rect->GetMin() + dir;
221 }
222 else
223 {
224 nearer_subtree = node->GetRight();
225 farther_subtree = node->GetLeft();
226 nearer_hyperrect_coord = rect->GetMin() + dir;
227 farther_hyperrect_coord = rect->GetMax() + dir;
228 }
229
230 if (nearer_subtree)
231 {
232 /* Slice the hyperrect to get the hyperrect of the nearer subtree */
233 dummy = *nearer_hyperrect_coord;
234 *nearer_hyperrect_coord = (*node)[dir];
235 /* Recurse down into nearer subtree */
236 __NearestToPosition(nearer_subtree, pos, result, result_dist_sq, rect);
237 /* Undo the slice */
238 *nearer_hyperrect_coord = dummy;
239 }
240
241 /* Check the distance of the point at the current node, compare it
242 * with our best so far */
243 if (node->IsValid()) // TODO
244 {
245 dist_sq = 0;
246 bool do_break = false;
247 for (size_t i = 0; i < fDim; i++)
248 {
249 dist_sq += sqr((*node)[i] - pos[i]);
250 if (dist_sq > *result_dist_sq)
251 {
252 do_break = true;
253 break;
254 }
255 }
256 if (!do_break && dist_sq < *result_dist_sq)
257 {
258 result = node;
259 *result_dist_sq = dist_sq;
260 }
261 }
262
263 if (farther_subtree)
264 {
265 /* Get the hyperrect of the farther subtree */
266 dummy = *farther_hyperrect_coord;
267 *farther_hyperrect_coord = (*node)[dir];
268 /* Check if we have to recurse down by calculating the closest
269 * point of the hyperrect and see if it's closer than our
270 * minimum distance in result_dist_sq. */
271 if (rect->CompareDistSqr(pos, result_dist_sq))
272 {
273 /* Recurse down into farther subtree */
274 __NearestToPosition(farther_subtree, pos, result, result_dist_sq, rect);
275 }
276 /* Undo the slice on the hyperrect */
277 *farther_hyperrect_coord = dummy;
278 }
279 }
280
281template<typename Position>
282 G4KDTreeResultHandle G4KDTree::Nearest(const Position& pos)
283 {
284 // G4cout << "Nearest(pos)" << G4endl ;
285
286 if (!fRect) return 0;
287
288 G4KDNode_Base *result(0);
289 double dist_sq = DBL_MAX;
290
291 /* Duplicate the bounding hyperrectangle, we will work on the copy */
292 HyperRect* newrect = new HyperRect(*fRect);
293
294 /* Our first estimate is the root node */
295 /* Search for the nearest neighbour recursively */
296 __NearestToPosition(fRoot, pos, result, &dist_sq, newrect);
297
298 /* Free the copy of the hyperrect */
299 delete newrect;
300
301 /* Store the result */
302 if (result)
303 {
304 G4KDTreeResultHandle rset = new G4KDTreeResult(this);
305 rset->Insert(dist_sq, result);
306 rset->Rewind();
307 return rset;
308 }
309 else
310 {
311 return 0;
312 }
313 }
314
315//__________________________________________________________________
316template<typename Position>
317 void G4KDTree::__NearestToNode(G4KDNode_Base* source_node,
318 G4KDNode_Base* node,
319 const Position& pos,
320 std::vector<G4KDNode_Base*>& result,
321 double *result_dist_sq,
322 HyperRect* rect,
323 int& nbresult)
324 {
325 int dir = node->GetAxis();
326 double dummy, dist_sq;
327 G4KDNode_Base *nearer_subtree(0), *farther_subtree(0);
328 double *nearer_hyperrect_coord(0), *farther_hyperrect_coord(0);
329
330 /* Decide whether to go left or right in the tree */
331 dummy = pos[dir] - (*node)[dir];
332 if (dummy <= 0)
333 {
334 nearer_subtree = node->GetLeft();
335 farther_subtree = node->GetRight();
336 nearer_hyperrect_coord = rect->GetMax() + dir;
337 farther_hyperrect_coord = rect->GetMin() + dir;
338 }
339 else
340 {
341 nearer_subtree = node->GetRight();
342 farther_subtree = node->GetLeft();
343 nearer_hyperrect_coord = rect->GetMin() + dir;
344 farther_hyperrect_coord = rect->GetMax() + dir;
345 }
346
347 if (nearer_subtree)
348 {
349 /* Slice the hyperrect to get the hyperrect of the nearer subtree */
350 dummy = *nearer_hyperrect_coord;
351 *nearer_hyperrect_coord = (*node)[dir];
352 /* Recurse down into nearer subtree */
353 __NearestToNode(source_node, nearer_subtree, pos, result, result_dist_sq,
354 rect, nbresult);
355 /* Undo the slice */
356 *nearer_hyperrect_coord = dummy;
357 }
358
359 /* Check the distance of the point at the current node, compare it
360 * with our best so far */
361 if (node->IsValid() && node != source_node)
362 {
363 dist_sq = 0;
364 bool do_break = false;
365 for (size_t i = 0; i < fDim; i++)
366 {
367 dist_sq += sqr((*node)[i] - pos[i]);
368 if (dist_sq > *result_dist_sq)
369 {
370 do_break = true;
371 break;
372 }
373 }
374 if (!do_break)
375 {
376 if (dist_sq < *result_dist_sq)
377 {
378 result.clear();
379 nbresult = 1;
380 result.push_back(node);
381 *result_dist_sq = dist_sq;
382 }
383 else if (dist_sq == *result_dist_sq)
384 {
385 result.push_back(node);
386 nbresult++;
387 }
388 }
389 }
390
391 if (farther_subtree)
392 {
393 /* Get the hyperrect of the farther subtree */
394 dummy = *farther_hyperrect_coord;
395 *farther_hyperrect_coord = (*node)[dir];
396 /* Check if we have to recurse down by calculating the closest
397 * point of the hyperrect and see if it's closer than our
398 * minimum distance in result_dist_sq. */
399 // if (hyperrect_dist_sq(rect, pos) < *result_dist_sq)
400 if (rect->CompareDistSqr(pos, result_dist_sq))
401 {
402 /* Recurse down into farther subtree */
403 __NearestToNode(source_node, farther_subtree, pos, result,
404 result_dist_sq, rect, nbresult);
405 }
406 /* Undo the slice on the hyperrect */
407 *farther_hyperrect_coord = dummy;
408 }
409 }
410
411template<typename Position>
412 G4KDTreeResultHandle G4KDTree::NearestInRange(const Position& pos,
413 const double& range)
414 {
415 int ret(-1);
416
417 const double range_sq = sqr(range);
418
419 G4KDTreeResultHandle rset = new G4KDTreeResult(this);
420 if ((ret = __NearestInRange(fRoot, pos, range_sq, range, *(rset()), 0)) == -1)
421 {
422 rset = 0;
423 return rset;
424 }
425 rset->Sort();
426 rset->Rewind();
427 return rset;
428 }
Generated by doxygen 1.9.3