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

并行计算

开发平台:

Visual C++

  1. //----------------------------------------------------------------------
  2. // File: kd_util.cpp
  3. // Programmer: Sunil Arya and David Mount
  4. // Description: Common utilities 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_util.h" // kd-utility declarations
  25. #include "../ANNperf.h" // performance evaluation
  26. //----------------------------------------------------------------------
  27. // The following routines are utility functions for manipulating
  28. // points sets, used in determining splitting planes for kd-tree
  29. // construction.
  30. //----------------------------------------------------------------------
  31. //----------------------------------------------------------------------
  32. // NOTE: Virtually all point indexing is done through an index (i.e.
  33. // permutation) array pidx.  Consequently, a reference to the d-th
  34. // coordinate of the i-th point is pa[pidx[i]][d].  The macro PA(i,d)
  35. // is a shorthand for this.
  36. //----------------------------------------------------------------------
  37. // standard 2-d indirect indexing
  38. #define PA(i,d) (pa[pidx[(i)]][(d)])
  39. // accessing a single point
  40. #define PP(i) (pa[pidx[(i)]])
  41. //----------------------------------------------------------------------
  42. // annAspectRatio
  43. // Compute the aspect ratio (ratio of longest to shortest side)
  44. // of a rectangle.
  45. //----------------------------------------------------------------------
  46. double annAspectRatio(
  47. int dim, // dimension
  48. const ANNorthRect &bnd_box) // bounding cube
  49. {
  50. ANNcoord length = bnd_box.hi[0] - bnd_box.lo[0];
  51. ANNcoord min_length = length; // min side length
  52. ANNcoord max_length = length; // max side length
  53. for (int d = 0; d < dim; d++) {
  54. length = bnd_box.hi[d] - bnd_box.lo[d];
  55. if (length < min_length) min_length = length;
  56. if (length > max_length) max_length = length;
  57. }
  58. return max_length/min_length;
  59. }
  60. //----------------------------------------------------------------------
  61. // annEnclRect, annEnclCube
  62. // These utilities compute the smallest rectangle and cube enclosing
  63. // a set of points, respectively.
  64. //----------------------------------------------------------------------
  65. void annEnclRect(
  66. ANNpointArray pa, // point array
  67. ANNidxArray pidx, // point indices
  68. int n, // number of points
  69. int dim, // dimension
  70. ANNorthRect &bnds) // bounding cube (returned)
  71. {
  72. for (int d = 0; d < dim; d++) { // find smallest enclosing rectangle
  73. ANNcoord lo_bnd = PA(0,d); // lower bound on dimension d
  74. ANNcoord hi_bnd = PA(0,d); // upper bound on dimension d
  75. for (int i = 0; i < n; i++) {
  76. if (PA(i,d) < lo_bnd) lo_bnd = PA(i,d);
  77. else if (PA(i,d) > hi_bnd) hi_bnd = PA(i,d);
  78. }
  79. bnds.lo[d] = lo_bnd;
  80. bnds.hi[d] = hi_bnd;
  81. }
  82. }
  83. void annEnclCube( // compute smallest enclosing cube
  84. ANNpointArray pa, // point array
  85. ANNidxArray pidx, // point indices
  86. int n, // number of points
  87. int dim, // dimension
  88. ANNorthRect &bnds) // bounding cube (returned)
  89. {
  90. int d;
  91. // compute smallest enclosing rect
  92. annEnclRect(pa, pidx, n, dim, bnds);
  93. ANNcoord max_len = 0; // max length of any side
  94. for (d = 0; d < dim; d++) { // determine max side length
  95. ANNcoord len = bnds.hi[d] - bnds.lo[d];
  96. if (len > max_len) { // update max_len if longest
  97. max_len = len;
  98. }
  99. }
  100. for (d = 0; d < dim; d++) { // grow sides to match max
  101. ANNcoord len = bnds.hi[d] - bnds.lo[d];
  102. ANNcoord half_diff = (max_len - len) / 2;
  103. bnds.lo[d] -= half_diff;
  104. bnds.hi[d] += half_diff;
  105. }
  106. }
  107. //----------------------------------------------------------------------
  108. // annBoxDistance - utility routine which computes distance from point to
  109. // box (Note: most distances to boxes are computed using incremental
  110. // distance updates, not this function.)
  111. //----------------------------------------------------------------------
  112. ANNdist annBoxDistance( // compute distance from point to box
  113. const ANNpoint q, // the point
  114. const ANNpoint lo, // low point of box
  115. const ANNpoint hi, // high point of box
  116. int dim) // dimension of space
  117. {
  118. register ANNdist dist = 0.0; // sum of squared distances
  119. register ANNdist t;
  120. for (register int d = 0; d < dim; d++) {
  121. if (q[d] < lo[d]) { // q is left of box
  122. t = ANNdist(lo[d]) - ANNdist(q[d]);
  123. dist = ANN_SUM(dist, ANN_POW(t));
  124. }
  125. else if (q[d] > hi[d]) { // q is right of box
  126. t = ANNdist(q[d]) - ANNdist(hi[d]);
  127. dist = ANN_SUM(dist, ANN_POW(t));
  128. }
  129. }
  130. ANN_FLOP(4*dim) // increment floating op count
  131. return dist;
  132. }
  133. //----------------------------------------------------------------------
  134. // annSpread - find spread along given dimension
  135. // annMinMax - find min and max coordinates along given dimension
  136. // annMaxSpread - find dimension of max spread
  137. //----------------------------------------------------------------------
  138. ANNcoord annSpread( // compute point spread along dimension
  139. ANNpointArray pa, // point array
  140. ANNidxArray pidx, // point indices
  141. int n, // number of points
  142. int d) // dimension to check
  143. {
  144. ANNcoord min = PA(0,d); // compute max and min coords
  145. ANNcoord max = PA(0,d);
  146. for (int i = 1; i < n; i++) {
  147. ANNcoord c = PA(i,d);
  148. if (c < min) min = c;
  149. else if (c > max) max = c;
  150. }
  151. return (max - min); // total spread is difference
  152. }
  153. void annMinMax( // compute min and max coordinates along dim
  154. ANNpointArray pa, // point array
  155. ANNidxArray pidx, // point indices
  156. int n, // number of points
  157. int d, // dimension to check
  158. ANNcoord &min, // minimum value (returned)
  159. ANNcoord &max) // maximum value (returned)
  160. {
  161. min = PA(0,d); // compute max and min coords
  162. max = PA(0,d);
  163. for (int i = 1; i < n; i++) {
  164. ANNcoord c = PA(i,d);
  165. if (c < min) min = c;
  166. else if (c > max) max = c;
  167. }
  168. }
  169. int annMaxSpread( // compute dimension of max spread
  170. ANNpointArray pa, // point array
  171. ANNidxArray pidx, // point indices
  172. int n, // number of points
  173. int dim) // dimension of space
  174. {
  175. int max_dim = 0; // dimension of max spread
  176. ANNcoord max_spr = 0; // amount of max spread
  177. if (n == 0) return max_dim; // no points, who cares?
  178. for (int d = 0; d < dim; d++) { // compute spread along each dim
  179. ANNcoord spr = annSpread(pa, pidx, n, d);
  180. if (spr > max_spr) { // bigger than current max
  181. max_spr = spr;
  182. max_dim = d;
  183. }
  184. }
  185. return max_dim;
  186. }
  187. //----------------------------------------------------------------------
  188. // annMedianSplit - split point array about its median
  189. // Splits a subarray of points pa[0..n] about an element of given
  190. // rank (median: n_lo = n/2) with respect to dimension d.  It places
  191. // the element of rank n_lo-1 correctly (because our splitting rule
  192. // takes the mean of these two).  On exit, the array is permuted so
  193. // that:
  194. //
  195. // pa[0..n_lo-2][d] <= pa[n_lo-1][d] <= pa[n_lo][d] <= pa[n_lo+1..n-1][d].
  196. //
  197. // The mean of pa[n_lo-1][d] and pa[n_lo][d] is returned as the
  198. // splitting value.
  199. //
  200. // All indexing is done indirectly through the index array pidx.
  201. //
  202. // This function uses the well known selection algorithm due to
  203. // C.A.R. Hoare.
  204. //----------------------------------------------------------------------
  205. // swap two points in pa array
  206. #define PASWAP(a,b) { int tmp = pidx[a]; pidx[a] = pidx[b]; pidx[b] = tmp; }
  207. void annMedianSplit(
  208. ANNpointArray pa, // points to split
  209. ANNidxArray pidx, // point indices
  210. int n, // number of points
  211. int d, // dimension along which to split
  212. ANNcoord &cv, // cutting value
  213. int n_lo) // split into n_lo and n-n_lo
  214. {
  215. int l = 0; // left end of current subarray
  216. int r = n-1; // right end of current subarray
  217. while (l < r) {
  218. register int i = (r+l)/2; // select middle as pivot
  219. register int k;
  220. if (PA(i,d) > PA(r,d)) // make sure last > pivot
  221. PASWAP(i,r)
  222. PASWAP(l,i); // move pivot to first position
  223. ANNcoord c = PA(l,d); // pivot value
  224. i = l;
  225. k = r;
  226. for(;;) { // pivot about c
  227. while (PA(++i,d) < c) ;
  228. while (PA(--k,d) > c) ;
  229. if (i < k) PASWAP(i,k) else break;
  230. }
  231. PASWAP(l,k); // pivot winds up in location k
  232. if (k > n_lo)    r = k-1; // recurse on proper subarray
  233. else if (k < n_lo) l = k+1;
  234. else break; // got the median exactly
  235. }
  236. if (n_lo > 0) { // search for next smaller item
  237. ANNcoord c = PA(0,d); // candidate for max
  238. int k = 0; // candidate's index
  239. for (int i = 1; i < n_lo; i++) {
  240. if (PA(i,d) > c) {
  241. c = PA(i,d);
  242. k = i;
  243. }
  244. }
  245. PASWAP(n_lo-1, k); // max among pa[0..n_lo-1] to pa[n_lo-1]
  246. }
  247. // cut value is midpoint value
  248. cv = (PA(n_lo-1,d) + PA(n_lo,d))/2.0;
  249. }
  250. //----------------------------------------------------------------------
  251. // annPlaneSplit - split point array about a cutting plane
  252. // Split the points in an array about a given plane along a
  253. // given cutting dimension.  On exit, br1 and br2 are set so
  254. // that:
  255. //
  256. // pa[ 0 ..br1-1] <  cv
  257. // pa[br1..br2-1] == cv
  258. // pa[br2.. n -1] >  cv
  259. //
  260. // All indexing is done indirectly through the index array pidx.
  261. //
  262. //----------------------------------------------------------------------
  263. void annPlaneSplit( // split points by a plane
  264. ANNpointArray pa, // points to split
  265. ANNidxArray pidx, // point indices
  266. int n, // number of points
  267. int d, // dimension along which to split
  268. ANNcoord cv, // cutting value
  269. int &br1, // first break (values < cv)
  270. int &br2) // second break (values == cv)
  271. {
  272. int l = 0;
  273. int r = n-1;
  274. for(;;) { // partition pa[0..n-1] about cv
  275. while (l < n && PA(l,d) < cv) l++;
  276. while (r >= 0 && PA(r,d) >= cv) r--;
  277. if (l > r) break;
  278. PASWAP(l,r);
  279. l++; r--;
  280. }
  281. br1 = l; // now: pa[0..br1-1] < cv <= pa[br1..n-1]
  282. r = n-1;
  283. for(;;) { // partition pa[br1..n-1] about cv
  284. while (l < n && PA(l,d) <= cv) l++;
  285. while (r >= br1 && PA(r,d) > cv) r--;
  286. if (l > r) break;
  287. PASWAP(l,r);
  288. l++; r--;
  289. }
  290. br2 = l; // now: pa[br1..br2-1] == cv < pa[br2..n-1]
  291. }
  292. //----------------------------------------------------------------------
  293. // annBoxSplit - split point array about a orthogonal rectangle
  294. // Split the points in an array about a given orthogonal
  295. // rectangle.  On exit, n_in is set to the number of points
  296. // that are inside (or on the boundary of) the rectangle.
  297. //
  298. // All indexing is done indirectly through the index array pidx.
  299. //
  300. //----------------------------------------------------------------------
  301. void annBoxSplit( // split points by a box
  302. ANNpointArray pa, // points to split
  303. ANNidxArray pidx, // point indices
  304. int n, // number of points
  305. int dim, // dimension of space
  306. ANNorthRect &box, // the box
  307. int &n_in) // number of points inside (returned)
  308. {
  309. int l = 0;
  310. int r = n-1;
  311. for(;;) { // partition pa[0..n-1] about box
  312. while (l < n && box.inside(dim, PP(l))) l++;
  313. while (r >= 0 && !box.inside(dim, PP(r))) r--;
  314. if (l > r) break;
  315. PASWAP(l,r);
  316. l++; r--;
  317. }
  318. n_in = l; // now: pa[0..n_in-1] inside and rest outside
  319. }
  320. //----------------------------------------------------------------------
  321. // annSplitBalance - compute balance factor for a given plane split
  322. // Balance factor is defined as the number of points lying
  323. // below the splitting value minus n/2 (median).  Thus, a
  324. // median split has balance 0, left of this is negative and
  325. // right of this is positive.  (The points are unchanged.)
  326. //----------------------------------------------------------------------
  327. int annSplitBalance( // determine balance factor of a split
  328. ANNpointArray pa, // points to split
  329. ANNidxArray pidx, // point indices
  330. int n, // number of points
  331. int d, // dimension along which to split
  332. ANNcoord cv) // cutting value
  333. {
  334. int n_lo = 0;
  335. for(int i = 0; i < n; i++) { // count number less than cv
  336. if (PA(i,d) < cv) n_lo++;
  337. }
  338. return n_lo - n/2;
  339. }
  340. //----------------------------------------------------------------------
  341. // annBox2Bnds - convert bounding box to list of bounds
  342. // Given two boxes, an inner box enclosed within a bounding
  343. // box, this routine determines all the sides for which the
  344. // inner box is strictly contained with the bounding box,
  345. // and adds an appropriate entry to a list of bounds.  Then
  346. // we allocate storage for the final list of bounds, and return
  347. // the resulting list and its size.
  348. //----------------------------------------------------------------------
  349. void annBox2Bnds( // convert inner box to bounds
  350. const ANNorthRect &inner_box, // inner box
  351. const ANNorthRect &bnd_box, // enclosing box
  352. int dim, // dimension of space
  353. int &n_bnds, // number of bounds (returned)
  354. ANNorthHSArray &bnds) // bounds array (returned)
  355. {
  356. int i;
  357. n_bnds = 0; // count number of bounds
  358. for (i = 0; i < dim; i++) {
  359. if (inner_box.lo[i] > bnd_box.lo[i]) // low bound is inside
  360. n_bnds++;
  361. if (inner_box.hi[i] < bnd_box.hi[i]) // high bound is inside
  362. n_bnds++;
  363. }
  364. bnds = new ANNorthHalfSpace[n_bnds]; // allocate appropriate size
  365. int j = 0;
  366. for (i = 0; i < dim; i++) { // fill the array
  367. if (inner_box.lo[i] > bnd_box.lo[i]) {
  368. bnds[j].cd = i;
  369. bnds[j].cv = inner_box.lo[i];
  370. bnds[j].sd = +1;
  371. j++;
  372. }
  373. if (inner_box.hi[i] < bnd_box.hi[i]) {
  374. bnds[j].cd = i;
  375. bnds[j].cv = inner_box.hi[i];
  376. bnds[j].sd = -1;
  377. j++;
  378. }
  379. }
  380. }
  381. //----------------------------------------------------------------------
  382. // annBnds2Box - convert list of bounds to bounding box
  383. // Given an enclosing box and a list of bounds, this routine
  384. // computes the corresponding inner box.  It is assumed that
  385. // the box points have been allocated already.
  386. //----------------------------------------------------------------------
  387. void annBnds2Box(
  388. const ANNorthRect &bnd_box, // enclosing box
  389. int dim, // dimension of space
  390. int n_bnds, // number of bounds
  391. ANNorthHSArray bnds, // bounds array
  392. ANNorthRect &inner_box) // inner box (returned)
  393. {
  394. annAssignRect(dim, inner_box, bnd_box); // copy bounding box to inner
  395. for (int i = 0; i < n_bnds; i++) {
  396. bnds[i].project(inner_box.lo); // project each endpoint
  397. bnds[i].project(inner_box.hi);
  398. }
  399. }