Geogram  Version 1.9.1
A programming library of geometric algorithms
kd_tree.h
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2000-2022 Inria
3  * All rights reserved.
4  *
5  * Redistribution and use in source and binary forms, with or without
6  * modification, are permitted provided that the following conditions are met:
7  *
8  * * Redistributions of source code must retain the above copyright notice,
9  * this list of conditions and the following disclaimer.
10  * * Redistributions in binary form must reproduce the above copyright notice,
11  * this list of conditions and the following disclaimer in the documentation
12  * and/or other materials provided with the distribution.
13  * * Neither the name of the ALICE Project-Team nor the names of its
14  * contributors may be used to endorse or promote products derived from this
15  * software without specific prior written permission.
16  *
17  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27  * POSSIBILITY OF SUCH DAMAGE.
28  *
29  * Contact: Bruno Levy
30  *
31  * https://www.inria.fr/fr/bruno-levy
32  *
33  * Inria,
34  * Domaine de Voluceau,
35  * 78150 Le Chesnay - Rocquencourt
36  * FRANCE
37  *
38  */
39 
40 #ifndef GEOGRAM_POINTS_KD_TREE
41 #define GEOGRAM_POINTS_KD_TREE
42 
43 #include <geogram/basic/common.h>
45 #include <algorithm>
46 
52 namespace GEO {
53 
57  class GEOGRAM_API KdTree : public NearestNeighborSearch {
58  public:
64 
66  void set_points(index_t nb_points, const double* points) override;
67 
69  bool stride_supported() const override;
70 
72  void set_points(
73  index_t nb_points, const double* points, index_t stride
74  ) override;
75 
78  index_t nb_neighbors,
79  const double* query_point,
80  index_t* neighbors,
81  double* neighbors_sq_dist
82  ) const override;
83 
86  index_t nb_neighbors,
87  const double* query_point,
88  index_t* neighbors,
89  double* neighbors_sq_dist,
91  ) const override;
92 
95  index_t nb_neighbors,
96  index_t query_point,
97  index_t* neighbors,
98  double* neighbors_sq_dist
99  ) const override;
100 
101  /**********************************************************************/
102 
109 
130  index_t nb_neighbors_in,
131  index_t* user_neighbors_in,
132  double* user_neighbors_sq_dist_in,
133  index_t* work_neighbors_in,
134  double* work_neighbors_sq_dist_in
135  ) :
136  nb_neighbors(0),
137  nb_neighbors_max(nb_neighbors_in),
138  neighbors(work_neighbors_in),
139  neighbors_sq_dist(work_neighbors_sq_dist_in),
140  user_neighbors(user_neighbors_in),
141  user_neighbors_sq_dist(user_neighbors_sq_dist_in),
142  nb_visited(0)
143  {
144  // Yes, '<=' because we got space for n+1 neigbors
145  // in the work arrays.
146  for(index_t i = 0; i <= nb_neighbors; ++i) {
147  neighbors[i] = index_t(-1);
148  neighbors_sq_dist[i] = Numeric::max_float64();
149  }
150  }
151 
156  double furthest_neighbor_sq_dist() const {
157  return
158  nb_neighbors == nb_neighbors_max ?
159  neighbors_sq_dist[nb_neighbors - 1] :
161  ;
162  }
163 
173  void insert(
174  index_t neighbor, double sq_dist
175  ) {
177  sq_dist <= furthest_neighbor_sq_dist()
178  );
179 
180  int i;
181  for(i=int(nb_neighbors); i>0; --i) {
182  if(neighbors_sq_dist[i - 1] < sq_dist) {
183  break;
184  }
185  neighbors[i] = neighbors[i - 1];
186  neighbors_sq_dist[i] = neighbors_sq_dist[i - 1];
187  }
188 
189  neighbors[i] = neighbor;
190  neighbors_sq_dist[i] = sq_dist;
191 
192  if(nb_neighbors < nb_neighbors_max) {
193  ++nb_neighbors;
194  }
195  }
196 
204  void copy_from_user() {
205  for(index_t i=0; i<nb_neighbors_max; ++i) {
206  neighbors[i] = user_neighbors[i];
207  neighbors_sq_dist[i] = user_neighbors_sq_dist[i];
208  }
209  neighbors[nb_neighbors_max] = index_t(-1);
210  neighbors_sq_dist[nb_neighbors_max] = Numeric::max_float64();
211  nb_neighbors = nb_neighbors_max;
212  }
213 
220  void copy_to_user() {
221  for(index_t i=0; i<nb_neighbors_max; ++i) {
222  user_neighbors[i] = neighbors[i];
223  user_neighbors_sq_dist[i] = neighbors_sq_dist[i];
224  }
225  }
226 
229 
232 
238 
244 
250 
257 
262  size_t nb_visited;
263  };
264 
296  index_t node_index, index_t b, index_t e,
297  double* bbox_min, double* bbox_max,
298  double bbox_dist, const double* query_point,
299  NearestNeighbors& neighbors
300  ) const;
301 
318  double* bbox_min, double* bbox_max,
319  double& box_dist, const double* query_point
320  ) const;
321 
326  index_t root() const {
327  return root_;
328  }
329 
330  protected:
334  static const index_t MAX_LEAF_SIZE = 16;
335 
340  virtual index_t build_tree() = 0 ;
341 
361  virtual void get_node(
362  index_t n, index_t b, index_t e,
363  index_t& left_child, index_t& right_child,
364  coord_index_t& splitting_coord,
365  index_t& m,
366  double& splitting_val
367  ) const = 0;
368 
369 
370 
385  index_t node_index, index_t b, index_t e,
386  const double* query_point,
387  NearestNeighbors& neighbors
388  ) const;
389 
399  index_t b, index_t e, coord_index_t coord,
400  double& minval, double& maxval
401  ) const {
402  minval = Numeric::max_float64();
403  maxval = Numeric::min_float64();
404  for(index_t i = b; i < e; ++i) {
405  double val = point_ptr(point_index_[i])[coord];
406  minval = std::min(minval, val);
407  maxval = std::max(maxval, val);
408  }
409  }
410 
419  double spread(index_t b, index_t e, coord_index_t coord) const {
420  double minval,maxval;
421  get_minmax(b,e,coord,minval,maxval);
422  return maxval - minval;
423  }
424 
428  ~KdTree() override;
429 
430  protected:
431  vector<index_t> point_index_;
432  vector<double> bbox_min_;
433  vector<double> bbox_max_;
434  index_t root_;
435  };
436 
437  /*********************************************************************/
438 
447  class GEOGRAM_API BalancedKdTree : public KdTree {
448  public:
454 
455  protected:
459  ~BalancedKdTree() override;
460 
469  index_t node_id, index_t b, index_t e
470  ) {
471  if(e - b <= MAX_LEAF_SIZE) {
472  return node_id;
473  }
474  index_t m = b + (e - b) / 2;
475  return std::max(
476  max_node_index(2 * node_id, b, m),
477  max_node_index(2 * node_id + 1, m, e)
478  );
479  }
480 
488 
498  index_t node_index, index_t b, index_t e
499  ) {
500  if(e - b <= MAX_LEAF_SIZE) {
501  return;
502  }
503  index_t m = split_kd_node(node_index, b, e);
504  create_kd_tree_recursive(2 * node_index, b, m);
505  create_kd_tree_recursive(2 * node_index + 1, m, e);
506  }
507 
519  index_t node_index, index_t b, index_t e
520  );
521 
523  index_t build_tree() override;
524 
526  void get_node(
527  index_t n, index_t b, index_t e,
528  index_t& left_child, index_t& right_child,
529  coord_index_t& splitting_coord,
530  index_t& m,
531  double& splitting_val
532  ) const override;
533 
534  protected:
535 
540 
545 
549  index_t m0_, m1_, m2_, m3_, m4_, m5_, m6_, m7_, m8_;
550  };
551 
552  /*********************************************************************/
553 
570  class GEOGRAM_API AdaptiveKdTree : public KdTree {
571  public:
577 
578  protected:
580  index_t build_tree() override;
581 
583  void get_node(
584  index_t n, index_t b, index_t e,
585  index_t& left_child, index_t& right_child,
586  coord_index_t& splitting_coord,
587  index_t& m,
588  double& splitting_val
589  ) const override;
590 
603  index_t b, index_t e,
604  double* bbox_min, double* bbox_max
605  );
606 
622  virtual void split_kd_node(
623  index_t b, index_t e,
624  double* bbox_min, double* bbox_max,
625  index_t& m, coord_index_t& cut_dim, double& cut_val
626  );
627 
642  virtual void plane_split(
643  index_t b, index_t e, coord_index_t coord, double val,
644  index_t& br1, index_t& br2
645  );
646 
653  double point_coord(int index, coord_index_t coord) {
654  geo_debug_assert(index >= 0);
655  geo_debug_assert(index_t(index) < nb_points());
656  geo_debug_assert(coord < dimension());
657  index_t direct_index = point_index_[index_t(index)];
658  geo_debug_assert(direct_index < nb_points());
659  return (points_ + direct_index * stride_)[coord];
660  }
661 
662 
667  index_t nb_nodes() const {
668  return splitting_coord_.size();
669  }
670 
675  virtual index_t new_node();
676 
677  protected:
682 
687 
695 
701  };
702 
703  /*********************************************************************/
704 }
705 
706 #endif
#define geo_debug_assert(x)
Verifies that a condition is met.
Definition: assert.h:196
Implements NearestNeighborSearch using an Adaptive Kd-tree.
Definition: kd_tree.h:570
virtual void split_kd_node(index_t b, index_t e, double *bbox_min, double *bbox_max, index_t &m, coord_index_t &cut_dim, double &cut_val)
Computes and stores the splitting coordinate and splitting value of the node node_index,...
vector< index_t > node_m_
One per node, node splitting index.
Definition: kd_tree.h:694
double point_coord(int index, coord_index_t coord)
Gets a point coordinate by index and coordinate.
Definition: kd_tree.h:653
virtual index_t create_kd_tree_recursive(index_t b, index_t e, double *bbox_min, double *bbox_max)
Creates the subtree under a node.
index_t build_tree() override
Builds the tree.
vector< double > splitting_val_
One per node, splitting coordinate value.
Definition: kd_tree.h:686
virtual index_t new_node()
Creates a new node.
vector< index_t > node_right_child_
One per node, right child index.
Definition: kd_tree.h:700
AdaptiveKdTree(coord_index_t dim)
Creates a new BalancedKdTree.
vector< coord_index_t > splitting_coord_
One per node, splitting coordinate.
Definition: kd_tree.h:681
virtual void plane_split(index_t b, index_t e, coord_index_t coord, double val, index_t &br1, index_t &br2)
Reorders the points in a sequence in such a way that the specified coordinate in the beginning of the...
index_t nb_nodes() const
Gets the number of nodes.
Definition: kd_tree.h:667
void get_node(index_t n, index_t b, index_t e, index_t &left_child, index_t &right_child, coord_index_t &splitting_coord, index_t &m, double &splitting_val) const override
Gets all the attributes of a node.
Implements NearestNeighborSearch using a balanced Kd-tree.
Definition: kd_tree.h:447
void get_node(index_t n, index_t b, index_t e, index_t &left_child, index_t &right_child, coord_index_t &splitting_coord, index_t &m, double &splitting_val) const override
Gets all the attributes of a node.
void create_kd_tree_recursive(index_t node_index, index_t b, index_t e)
Creates the subtree under a node.
Definition: kd_tree.h:497
BalancedKdTree(coord_index_t dim)
Creates a new BalancedKdTree.
vector< coord_index_t > splitting_coord_
One per node, splitting coordinate.
Definition: kd_tree.h:539
~BalancedKdTree() override
BalancedKdTree destructor.
index_t split_kd_node(index_t node_index, index_t b, index_t e)
Computes and stores the splitting coordinate and splitting value of the node node_index,...
coord_index_t best_splitting_coord(index_t b, index_t e)
Computes the coordinate along which a point sequence will be split.
index_t build_tree() override
Builds the tree.
vector< double > splitting_val_
One per node, splitting coordinate value.
Definition: kd_tree.h:544
index_t m0_
Indices for multithreaded tree construction.
Definition: kd_tree.h:549
static index_t max_node_index(index_t node_id, index_t b, index_t e)
Returns the maximum node index in subtree.
Definition: kd_tree.h:468
Base class for all Kd-tree implementations.
Definition: kd_tree.h:57
index_t root() const
Gets the root node.
Definition: kd_tree.h:326
void get_nearest_neighbors(index_t nb_neighbors, index_t query_point, index_t *neighbors, double *neighbors_sq_dist) const override
Finds the nearest neighbors of a point given by coordinates.
bool stride_supported() const override
Tests whether the stride variant of set_points() is supported.
virtual index_t build_tree()=0
Builds the tree.
void get_minmax(index_t b, index_t e, coord_index_t coord, double &minval, double &maxval) const
Computes the minimum and maximum point coordinates along a coordinate.
Definition: kd_tree.h:398
virtual void get_nearest_neighbors_recursive(index_t node_index, index_t b, index_t e, double *bbox_min, double *bbox_max, double bbox_dist, const double *query_point, NearestNeighbors &neighbors) const
The recursive function to implement KdTree traversal and nearest neighbors computation.
void set_points(index_t nb_points, const double *points, index_t stride) override
Sets the points and create the search data structure.
void get_nearest_neighbors(index_t nb_neighbors, const double *query_point, index_t *neighbors, double *neighbors_sq_dist, KeepInitialValues) const override
Finds the nearest neighbors of a point given by coordinates.
KdTree(coord_index_t dim)
KdTree constructor.
void set_points(index_t nb_points, const double *points) override
Sets the points and create the search data structure.
double spread(index_t b, index_t e, coord_index_t coord) const
Computes the extent of a point sequence along a given coordinate.
Definition: kd_tree.h:419
void init_bbox_and_bbox_dist_for_traversal(double *bbox_min, double *bbox_max, double &box_dist, const double *query_point) const
Initializes bounding box and box distance for Kd-Tree traversal.
virtual void get_nearest_neighbors_leaf(index_t node_index, index_t b, index_t e, const double *query_point, NearestNeighbors &neighbors) const
The recursive function to implement KdTree traversal and nearest neighbors computation in a leaf.
void get_nearest_neighbors(index_t nb_neighbors, const double *query_point, index_t *neighbors, double *neighbors_sq_dist) const override
Finds the nearest neighbors of a point given by coordinates.
virtual void get_node(index_t n, index_t b, index_t e, index_t &left_child, index_t &right_child, coord_index_t &splitting_coord, index_t &m, double &splitting_val) const =0
Gets all the attributes of a node.
~KdTree() override
KdTree destructor.
Abstract interface for nearest neighbor search algorithms.
Definition: nn_search.h:69
Common include file, providing basic definitions. Should be included before anything else by all head...
float64 max_float64()
Gets 64 bits float maximum positive value.
Definition: numeric.h:172
float64 min_float64()
Gets 64 bits float minimum negative value.
Definition: numeric.h:179
Global Vorpaline namespace.
Definition: basic.h:55
geo_index_t index_t
The type for storing and manipulating indices.
Definition: numeric.h:329
geo_coord_index_t coord_index_t
The type for storing coordinate indices, and iterating on the coordinates of a point.
Definition: numeric.h:363
Abstract interface for nearest neighbor searching.
The context for traversing a KdTree.
Definition: kd_tree.h:108
void copy_from_user()
Copies the user neighbors and distances into the work zone and initializes nb_neighbors to max_nb_nei...
Definition: kd_tree.h:204
double furthest_neighbor_sq_dist() const
Gets the squared distance to the furthest neighbor.
Definition: kd_tree.h:156
index_t nb_neighbors
Current number of neighbors.
Definition: kd_tree.h:228
index_t * neighbors
Internal array of neighbors.
Definition: kd_tree.h:237
index_t nb_neighbors_max
Maximum number of neighbors.
Definition: kd_tree.h:231
double * user_neighbors_sq_dist
User-provided array of neighbors squared distances.
Definition: kd_tree.h:256
void insert(index_t neighbor, double sq_dist)
Inserts a new neighbor.
Definition: kd_tree.h:173
size_t nb_visited
Number of points visited during traversal.
Definition: kd_tree.h:262
double * neighbors_sq_dist
Internal squared distance to neigbors.
Definition: kd_tree.h:243
index_t * user_neighbors
User-provided array of neighbors.
Definition: kd_tree.h:249
NearestNeighbors(index_t nb_neighbors_in, index_t *user_neighbors_in, double *user_neighbors_sq_dist_in, index_t *work_neighbors_in, double *work_neighbors_sq_dist_in)
Creates a new NearestNeighbors.
Definition: kd_tree.h:129
void copy_to_user()
Copies the found nearest neighbors from the work zone to the user neighbors and squared distance arra...
Definition: kd_tree.h:220
A structure to discriminate between the two versions of get_nearest_neighbors()
Definition: nn_search.h:141