00001 //---------------------------------------------------------------------- 00002 // File: kd_pr_search.cpp 00003 // Programmer: Sunil Arya and David Mount 00004 // Description: Priority search for kd-trees 00005 // Last modified: 01/04/05 (Version 1.0) 00006 //---------------------------------------------------------------------- 00007 // Copyright (c) 1997-2005 University of Maryland and Sunil Arya and 00008 // David Mount. All Rights Reserved. 00009 // 00010 // This software and related documentation is part of the Approximate 00011 // Nearest Neighbor Library (ANN). This software is provided under 00012 // the provisions of the Lesser GNU Public License (LGPL). See the 00013 // file ../ReadMe.txt for further information. 00014 // 00015 // The University of Maryland (U.M.) and the authors make no 00016 // representations about the suitability or fitness of this software for 00017 // any purpose. It is provided "as is" without express or implied 00018 // warranty. 00019 //---------------------------------------------------------------------- 00020 // History: 00021 // Revision 0.1 03/04/98 00022 // Initial release 00023 //---------------------------------------------------------------------- 00024 00025 #include "kd_pr_search.h" // kd priority search declarations 00026 00027 //---------------------------------------------------------------------- 00028 // Approximate nearest neighbor searching by priority search. 00029 // The kd-tree is searched for an approximate nearest neighbor. 00030 // The point is returned through one of the arguments, and the 00031 // distance returned is the SQUARED distance to this point. 00032 // 00033 // The method used for searching the kd-tree is called priority 00034 // search. (It is described in Arya and Mount, ``Algorithms for 00035 // fast vector quantization,'' Proc. of DCC '93: Data Compression 00036 // Conference}, eds. J. A. Storer and M. Cohn, IEEE Press, 1993, 00037 // 381--390.) 00038 // 00039 // The cell of the kd-tree containing the query point is located, 00040 // and cells are visited in increasing order of distance from the 00041 // query point. This is done by placing each subtree which has 00042 // NOT been visited in a priority queue, according to the closest 00043 // distance of the corresponding enclosing rectangle from the 00044 // query point. The search stops when the distance to the nearest 00045 // remaining rectangle exceeds the distance to the nearest point 00046 // seen by a factor of more than 1/(1+eps). (Implying that any 00047 // point found subsequently in the search cannot be closer by more 00048 // than this factor.) 00049 // 00050 // The main entry point is annkPriSearch() which sets things up and 00051 // then call the recursive routine ann_pri_search(). This is a 00052 // recursive routine which performs the processing for one node in 00053 // the kd-tree. There are two versions of this virtual procedure, 00054 // one for splitting nodes and one for leaves. When a splitting node 00055 // is visited, we determine which child to continue the search on 00056 // (the closer one), and insert the other child into the priority 00057 // queue. When a leaf is visited, we compute the distances to the 00058 // points in the buckets, and update information on the closest 00059 // points. 00060 // 00061 // Some trickery is used to incrementally update the distance from 00062 // a kd-tree rectangle to the query point. This comes about from 00063 // the fact that which each successive split, only one component 00064 // (along the dimension that is split) of the squared distance to 00065 // the child rectangle is different from the squared distance to 00066 // the parent rectangle. 00067 //---------------------------------------------------------------------- 00068 00069 //---------------------------------------------------------------------- 00070 // To keep argument lists short, a number of global variables 00071 // are maintained which are common to all the recursive calls. 00072 // These are given below. 00073 //---------------------------------------------------------------------- 00074 00075 double ANNprEps; // the error bound 00076 int ANNprDim; // dimension of space 00077 ANNpoint ANNprQ; // query point 00078 double ANNprMaxErr; // max tolerable squared error 00079 ANNpointArray ANNprPts; // the points 00080 ANNpr_queue *ANNprBoxPQ; // priority queue for boxes 00081 ANNmin_k *ANNprPointMK; // set of k closest points 00082 00083 //---------------------------------------------------------------------- 00084 // annkPriSearch - priority search for k nearest neighbors 00085 //---------------------------------------------------------------------- 00086 00087 void ANNkd_tree::annkPriSearch( 00088 ANNpoint q, // query point 00089 int k, // number of near neighbors to return 00090 ANNidxArray nn_idx, // nearest neighbor indices (returned) 00091 ANNdistArray dd, // dist to near neighbors (returned) 00092 double eps) // error bound (ignored) 00093 { 00094 // max tolerable squared error 00095 ANNprMaxErr = ANN_POW(1.0 + eps); 00096 ANN_FLOP(2) // increment floating ops 00097 00098 ANNprDim = dim; // copy arguments to static equivs 00099 ANNprQ = q; 00100 ANNprPts = pts; 00101 ANNptsVisited = 0; // initialize count of points visited 00102 00103 ANNprPointMK = new ANNmin_k(k); // create set for closest k points 00104 00105 // distance to root box 00106 ANNdist box_dist = annBoxDistance(q, 00107 bnd_box_lo, bnd_box_hi, dim); 00108 00109 ANNprBoxPQ = new ANNpr_queue(n_pts);// create priority queue for boxes 00110 ANNprBoxPQ->insert(box_dist, root); // insert root in priority queue 00111 00112 while (ANNprBoxPQ->non_empty() && 00113 (!(ANNmaxPtsVisited != 0 && ANNptsVisited > ANNmaxPtsVisited))) { 00114 ANNkd_ptr np; // next box from prior queue 00115 00116 // extract closest box from queue 00117 ANNprBoxPQ->extr_min(box_dist, (void *&) np); 00118 00119 ANN_FLOP(2) // increment floating ops 00120 if (box_dist*ANNprMaxErr >= ANNprPointMK->max_key()) 00121 break; 00122 00123 np->ann_pri_search(box_dist); // search this subtree. 00124 } 00125 00126 for (int i = 0; i < k; i++) { // extract the k-th closest points 00127 dd[i] = ANNprPointMK->ith_smallest_key(i); 00128 nn_idx[i] = ANNprPointMK->ith_smallest_info(i); 00129 } 00130 00131 delete ANNprPointMK; // deallocate closest point set 00132 delete ANNprBoxPQ; // deallocate priority queue 00133 } 00134 00135 //---------------------------------------------------------------------- 00136 // kd_split::ann_pri_search - search a splitting node 00137 //---------------------------------------------------------------------- 00138 00139 void ANNkd_split::ann_pri_search(ANNdist box_dist) 00140 { 00141 ANNdist new_dist; // distance to child visited later 00142 // distance to cutting plane 00143 ANNcoord cut_diff = ANNprQ[cut_dim] - cut_val; 00144 00145 if (cut_diff < 0) { // left of cutting plane 00146 ANNcoord box_diff = cd_bnds[ANN_LO] - ANNprQ[cut_dim]; 00147 if (box_diff < 0) // within bounds - ignore 00148 box_diff = 0; 00149 // distance to further box 00150 new_dist = (ANNdist) ANN_SUM(box_dist, 00151 ANN_DIFF(ANN_POW(box_diff), ANN_POW(cut_diff))); 00152 00153 if (child[ANN_HI] != KD_TRIVIAL)// enqueue if not trivial 00154 ANNprBoxPQ->insert(new_dist, child[ANN_HI]); 00155 // continue with closer child 00156 child[ANN_LO]->ann_pri_search(box_dist); 00157 } 00158 else { // right of cutting plane 00159 ANNcoord box_diff = ANNprQ[cut_dim] - cd_bnds[ANN_HI]; 00160 if (box_diff < 0) // within bounds - ignore 00161 box_diff = 0; 00162 // distance to further box 00163 new_dist = (ANNdist) ANN_SUM(box_dist, 00164 ANN_DIFF(ANN_POW(box_diff), ANN_POW(cut_diff))); 00165 00166 if (child[ANN_LO] != KD_TRIVIAL)// enqueue if not trivial 00167 ANNprBoxPQ->insert(new_dist, child[ANN_LO]); 00168 // continue with closer child 00169 child[ANN_HI]->ann_pri_search(box_dist); 00170 } 00171 ANN_SPL(1) // one more splitting node visited 00172 ANN_FLOP(8) // increment floating ops 00173 } 00174 00175 //---------------------------------------------------------------------- 00176 // kd_leaf::ann_pri_search - search points in a leaf node 00177 // 00178 // This is virtually identical to the ann_search for standard search. 00179 //---------------------------------------------------------------------- 00180 00181 void ANNkd_leaf::ann_pri_search(ANNdist box_dist) 00182 { 00183 register ANNdist dist; // distance to data point 00184 register ANNcoord* pp; // data coordinate pointer 00185 register ANNcoord* qq; // query coordinate pointer 00186 register ANNdist min_dist; // distance to k-th closest point 00187 register ANNcoord t; 00188 register int d; 00189 00190 min_dist = ANNprPointMK->max_key(); // k-th smallest distance so far 00191 00192 for (int i = 0; i < n_pts; i++) { // check points in bucket 00193 00194 pp = ANNprPts[bkt[i]]; // first coord of next data point 00195 qq = ANNprQ; // first coord of query point 00196 dist = 0; 00197 00198 for(d = 0; d < ANNprDim; d++) { 00199 ANN_COORD(1) // one more coordinate hit 00200 ANN_FLOP(4) // increment floating ops 00201 00202 t = *(qq++) - *(pp++); // compute length and adv coordinate 00203 // exceeds dist to k-th smallest? 00204 if( (dist = ANN_SUM(dist, ANN_POW(t))) > min_dist) { 00205 break; 00206 } 00207 } 00208 00209 if (d >= ANNprDim && // among the k best? 00210 (ANN_ALLOW_SELF_MATCH || dist!=0)) { // and no self-match problem 00211 // add it to the list 00212 ANNprPointMK->insert(dist, bkt[i]); 00213 min_dist = ANNprPointMK->max_key(); 00214 } 00215 } 00216 ANN_LEAF(1) // one more leaf node visited 00217 ANN_PTS(n_pts) // increment points visited 00218 ANNptsVisited += n_pts; // increment number of points visited 00219 }
1.6.1