kd_pr_search.cpp
上传用户:chinafayin
上传日期:2022-04-05
资源大小:153k
文件大小:9k
源码类别:

并行计算

开发平台:

Visual C++

  1. //----------------------------------------------------------------------
  2. // File: kd_pr_search.cpp
  3. // Programmer: Sunil Arya and David Mount
  4. // Description: Priority search for kd-trees
  5. // Last modified: 01/04/05 (Version 1.0)
  6. //----------------------------------------------------------------------
  7. // Copyright (c) 1997-2005 University of Maryland and Sunil Arya and
  8. // David Mount.  All Rights Reserved.
  9. // 
  10. // This software and related documentation is part of the Approximate
  11. // Nearest Neighbor Library (ANN).  This software is provided under
  12. // the provisions of the Lesser GNU Public License (LGPL).  See the
  13. // file ../ReadMe.txt for further information.
  14. // 
  15. // The University of Maryland (U.M.) and the authors make no
  16. // representations about the suitability or fitness of this software for
  17. // any purpose.  It is provided "as is" without express or implied
  18. // warranty.
  19. //----------------------------------------------------------------------
  20. // History:
  21. // Revision 0.1  03/04/98
  22. // Initial release
  23. //----------------------------------------------------------------------
  24. #include "kd_pr_search.h" // kd priority search declarations
  25. //----------------------------------------------------------------------
  26. // Approximate nearest neighbor searching by priority search.
  27. // The kd-tree is searched for an approximate nearest neighbor.
  28. // The point is returned through one of the arguments, and the
  29. // distance returned is the SQUARED distance to this point.
  30. //
  31. // The method used for searching the kd-tree is called priority
  32. // search.  (It is described in Arya and Mount, ``Algorithms for
  33. // fast vector quantization,'' Proc. of DCC '93: Data Compression
  34. // Conference}, eds. J. A. Storer and M. Cohn, IEEE Press, 1993,
  35. // 381--390.)
  36. //
  37. // The cell of the kd-tree containing the query point is located,
  38. // and cells are visited in increasing order of distance from the
  39. // query point.  This is done by placing each subtree which has
  40. // NOT been visited in a priority queue, according to the closest
  41. // distance of the corresponding enclosing rectangle from the
  42. // query point.  The search stops when the distance to the nearest
  43. // remaining rectangle exceeds the distance to the nearest point
  44. // seen by a factor of more than 1/(1+eps). (Implying that any
  45. // point found subsequently in the search cannot be closer by more
  46. // than this factor.)
  47. //
  48. // The main entry point is annkPriSearch() which sets things up and
  49. // then call the recursive routine ann_pri_search().  This is a
  50. // recursive routine which performs the processing for one node in
  51. // the kd-tree.  There are two versions of this virtual procedure,
  52. // one for splitting nodes and one for leaves. When a splitting node
  53. // is visited, we determine which child to continue the search on
  54. // (the closer one), and insert the other child into the priority
  55. // queue.  When a leaf is visited, we compute the distances to the
  56. // points in the buckets, and update information on the closest
  57. // points.
  58. //
  59. // Some trickery is used to incrementally update the distance from
  60. // a kd-tree rectangle to the query point.  This comes about from
  61. // the fact that which each successive split, only one component
  62. // (along the dimension that is split) of the squared distance to
  63. // the child rectangle is different from the squared distance to
  64. // the parent rectangle.
  65. //----------------------------------------------------------------------
  66. //----------------------------------------------------------------------
  67. // To keep argument lists short, a number of global variables
  68. // are maintained which are common to all the recursive calls.
  69. // These are given below.
  70. //----------------------------------------------------------------------
  71. ANN_THREAD_LOCAL double ANNprEps; // the error bound
  72. ANN_THREAD_LOCAL int ANNprDim; // dimension of space
  73. ANN_THREAD_LOCAL ANNpoint ANNprQ; // query point
  74. ANN_THREAD_LOCAL double ANNprMaxErr; // max tolerable squared error
  75. ANN_THREAD_LOCAL ANNpointArray ANNprPts; // the points
  76. ANN_THREAD_LOCAL ANNpr_queue *ANNprBoxPQ; // priority queue for boxes
  77. ANN_THREAD_LOCAL ANNmin_k *ANNprPointMK; // set of k closest points
  78. //----------------------------------------------------------------------
  79. // annkPriSearch - priority search for k nearest neighbors
  80. //----------------------------------------------------------------------
  81. void ANNkd_tree::annkPriSearch(
  82. ANNpoint q, // query point
  83. int k, // number of near neighbors to return
  84. ANNidxArray nn_idx, // nearest neighbor indices (returned)
  85. ANNdistArray dd, // dist to near neighbors (returned)
  86. double eps) // error bound (ignored)
  87. {
  88. // max tolerable squared error
  89. ANNprMaxErr = ANN_POW(1.0 + eps);
  90. ANN_FLOP(2) // increment floating ops
  91. ANNprDim = dim; // copy arguments to static equivs
  92. ANNprQ = q;
  93. ANNprPts = pts;
  94. ANNptsVisited = 0; // initialize count of points visited
  95. ANNprPointMK = new ANNmin_k(k); // create set for closest k points
  96. // distance to root box
  97. ANNdist box_dist = annBoxDistance(q,
  98. bnd_box_lo, bnd_box_hi, dim);
  99. ANNprBoxPQ = new ANNpr_queue(n_pts);// create priority queue for boxes
  100. ANNprBoxPQ->insert(box_dist, root); // insert root in priority queue
  101. while (ANNprBoxPQ->non_empty() &&
  102. (!(ANNmaxPtsVisited != 0 && ANNptsVisited > ANNmaxPtsVisited))) {
  103. ANNkd_ptr np; // next box from prior queue
  104. // extract closest box from queue
  105. ANNprBoxPQ->extr_min(box_dist, (void *&) np);
  106. ANN_FLOP(2) // increment floating ops
  107. if (box_dist*ANNprMaxErr >= ANNprPointMK->max_key())
  108. break;
  109. np->ann_pri_search(box_dist); // search this subtree.
  110. }
  111. for (int i = 0; i < k; i++) { // extract the k-th closest points
  112. dd[i] = ANNprPointMK->ith_smallest_key(i);
  113. nn_idx[i] = ANNprPointMK->ith_smallest_info(i);
  114. }
  115. delete ANNprPointMK; // deallocate closest point set
  116. delete ANNprBoxPQ; // deallocate priority queue
  117. }
  118. //----------------------------------------------------------------------
  119. // kd_split::ann_pri_search - search a splitting node
  120. //----------------------------------------------------------------------
  121. void ANNkd_split::ann_pri_search(ANNdist box_dist)
  122. {
  123. ANNdist new_dist; // distance to child visited later
  124. // distance to cutting plane
  125. ANNcoord cut_diff = ANNprQ[cut_dim] - cut_val;
  126. if (cut_diff < 0) { // left of cutting plane
  127. ANNcoord box_diff = cd_bnds[ANN_LO] - ANNprQ[cut_dim];
  128. if (box_diff < 0) // within bounds - ignore
  129. box_diff = 0;
  130. // distance to further box
  131. new_dist = (ANNdist) ANN_SUM(box_dist,
  132. ANN_DIFF(ANN_POW(box_diff), ANN_POW(cut_diff)));
  133. if (child[ANN_HI] != KD_TRIVIAL)// enqueue if not trivial
  134. ANNprBoxPQ->insert(new_dist, child[ANN_HI]);
  135. // continue with closer child
  136. child[ANN_LO]->ann_pri_search(box_dist);
  137. }
  138. else { // right of cutting plane
  139. ANNcoord box_diff = ANNprQ[cut_dim] - cd_bnds[ANN_HI];
  140. if (box_diff < 0) // within bounds - ignore
  141. box_diff = 0;
  142. // distance to further box
  143. new_dist = (ANNdist) ANN_SUM(box_dist,
  144. ANN_DIFF(ANN_POW(box_diff), ANN_POW(cut_diff)));
  145. if (child[ANN_LO] != KD_TRIVIAL)// enqueue if not trivial
  146. ANNprBoxPQ->insert(new_dist, child[ANN_LO]);
  147. // continue with closer child
  148. child[ANN_HI]->ann_pri_search(box_dist);
  149. }
  150. ANN_SPL(1) // one more splitting node visited
  151. ANN_FLOP(8) // increment floating ops
  152. }
  153. //----------------------------------------------------------------------
  154. // kd_leaf::ann_pri_search - search points in a leaf node
  155. //
  156. // This is virtually identical to the ann_search for standard search.
  157. //----------------------------------------------------------------------
  158. void ANNkd_leaf::ann_pri_search(ANNdist)
  159. {
  160. register ANNdist dist; // distance to data point
  161. register ANNcoord* pp; // data coordinate pointer
  162. register ANNcoord* qq; // query coordinate pointer
  163. register ANNdist min_dist; // distance to k-th closest point
  164. register ANNcoord t;
  165. register int d;
  166. min_dist = ANNprPointMK->max_key(); // k-th smallest distance so far
  167. for (int i = 0; i < n_pts; i++) { // check points in bucket
  168. pp = ANNprPts[bkt[i]]; // first coord of next data point
  169. qq = ANNprQ; // first coord of query point
  170. dist = 0;
  171. for(d = 0; d < ANNprDim; d++) {
  172. ANN_COORD(1) // one more coordinate hit
  173. ANN_FLOP(4) // increment floating ops
  174. t = *(qq++) - *(pp++); // compute length and adv coordinate
  175. // exceeds dist to k-th smallest?
  176. if( (dist = ANN_SUM(dist, ANN_POW(t))) > min_dist) {
  177. break;
  178. }
  179. }
  180. if (d >= ANNprDim && // among the k best?
  181.    (ANN_ALLOW_SELF_MATCH || dist!=0)) { // and no self-match problem
  182. // add it to the list
  183. ANNprPointMK->insert(dist, bkt[i]);
  184. min_dist = ANNprPointMK->max_key();
  185. }
  186. }
  187. ANN_LEAF(1) // one more leaf node visited
  188. ANN_PTS(n_pts) // increment points visited
  189. ANNptsVisited += n_pts; // increment number of points visited
  190. }