最近点算法

  1. type Point = (f64, f64);
  2. use std::cmp::Ordering;
  3. fn point_cmp((a1, a2): &Point, (b1, b2): &Point) -> Ordering {
  4. let acmp = f64_cmp(a1, b1);
  5. match acmp {
  6. Ordering::Equal => f64_cmp(a2, b2),
  7. _ => acmp,
  8. }
  9. }
  10. fn f64_cmp(a: &f64, b: &f64) -> Ordering {
  11. a.partial_cmp(b).unwrap()
  12. }
  13. /// returns the two closest points
  14. /// or None if there are zero or one point
  15. pub fn closest_points(points: &[Point]) -> Option<(Point, Point)> {
  16. let mut points: Vec<Point> = points.to_vec();
  17. points.sort_by(point_cmp);
  18. closest_points_aux(&points, 0, points.len())
  19. }
  20. fn dist((x1, y1): &Point, (x2, y2): &Point) -> f64 {
  21. let dx = *x1 - *x2;
  22. let dy = *y1 - *y2;
  23. (dx * dx + dy * dy).sqrt()
  24. }
  25. fn closest_points_aux(
  26. points: &[Point],
  27. mut start: usize,
  28. mut end: usize,
  29. ) -> Option<(Point, Point)> {
  30. let n = end - start;
  31. if n <= 1 {
  32. return None;
  33. }
  34. if n <= 3 {
  35. // bruteforce
  36. let mut min = dist(&points[0], &points[1]);
  37. let mut pair = (points[0], points[1]);
  38. for i in 0..n {
  39. for j in (i + 1)..n {
  40. let new = dist(&points[i], &points[j]);
  41. if new < min {
  42. min = new;
  43. pair = (points[i], points[j]);
  44. }
  45. }
  46. }
  47. return Some(pair);
  48. }
  49. let mid = (start + end) / 2;
  50. let left = closest_points_aux(points, start, mid);
  51. let right = closest_points_aux(points, mid, end);
  52. let (mut min_dist, mut pair) = match (left, right) {
  53. (Some((l1, l2)), Some((r1, r2))) => {
  54. let dl = dist(&l1, &l2);
  55. let dr = dist(&r1, &r2);
  56. if dl < dr {
  57. (dl, (l1, l2))
  58. } else {
  59. (dr, (r1, r2))
  60. }
  61. }
  62. (Some((a, b)), None) => (dist(&a, &b), (a, b)),
  63. (None, Some((a, b))) => (dist(&a, &b), (a, b)),
  64. (None, None) => unreachable!(),
  65. };
  66. let mid_x = points[mid].0;
  67. while points[start].0 < mid_x - min_dist {
  68. start += 1;
  69. }
  70. while points[end - 1].0 > mid_x + min_dist {
  71. end -= 1;
  72. }
  73. let mut mids: Vec<&Point> = points[start..end].iter().collect();
  74. mids.sort_by(|a, b| f64_cmp(&a.1, &b.1));
  75. for (i, e) in mids.iter().enumerate() {
  76. for k in 1..8 {
  77. if i + k >= mids.len() {
  78. break;
  79. }
  80. let new = dist(e, mids[i + k]);
  81. if new < min_dist {
  82. min_dist = new;
  83. pair = (**e, *mids[i + k]);
  84. }
  85. }
  86. }
  87. Some(pair)
  88. }
  89. #[cfg(test)]
  90. mod tests {
  91. use super::closest_points;
  92. use super::Point;
  93. fn eq(p1: Option<(Point, Point)>, p2: Option<(Point, Point)>) -> bool {
  94. match (p1, p2) {
  95. (None, None) => true,
  96. (Some((p1, p2)), Some((p3, p4))) => (p1 == p3 && p2 == p4) || (p1 == p4 && p2 == p3),
  97. _ => false,
  98. }
  99. }
  100. macro_rules! assert_display {
  101. ($left: expr, $right: expr) => {
  102. assert!(
  103. eq($left, $right),
  104. "assertion failed: `(left == right)`\nleft: `{:?}`,\nright: `{:?}`",
  105. $left,
  106. $right
  107. )
  108. };
  109. }
  110. #[test]
  111. fn zero_points() {
  112. let vals: [Point; 0] = [];
  113. assert_display!(closest_points(&vals), None::<(Point, Point)>);
  114. }
  115. #[test]
  116. fn one_points() {
  117. let vals = [(0., 0.)];
  118. assert_display!(closest_points(&vals), None::<(Point, Point)>);
  119. }
  120. #[test]
  121. fn two_points() {
  122. let vals = [(0., 0.), (1., 1.)];
  123. assert_display!(closest_points(&vals), Some(((0., 0.), (1., 1.))));
  124. }
  125. #[test]
  126. fn three_points() {
  127. let vals = [(0., 0.), (1., 1.), (3., 3.)];
  128. assert_display!(closest_points(&vals), Some(((0., 0.), (1., 1.))));
  129. }
  130. #[test]
  131. fn list_1() {
  132. let vals = [
  133. (0., 0.),
  134. (2., 1.),
  135. (5., 2.),
  136. (2., 3.),
  137. (4., 0.),
  138. (0., 4.),
  139. (5., 6.),
  140. (4., 4.),
  141. (7., 3.),
  142. (-1., 2.),
  143. (2., 6.),
  144. ];
  145. assert_display!(closest_points(&vals), Some(((2., 1.), (2., 3.))));
  146. }
  147. #[test]
  148. fn list_2() {
  149. let vals = [
  150. (1., 3.),
  151. (4., 6.),
  152. (8., 8.),
  153. (7., 5.),
  154. (5., 3.),
  155. (10., 3.),
  156. (7., 1.),
  157. (8., 3.),
  158. (4., 9.),
  159. (4., 12.),
  160. (4., 15.),
  161. (7., 14.),
  162. (8., 12.),
  163. (6., 10.),
  164. (4., 14.),
  165. (2., 7.),
  166. (3., 8.),
  167. (5., 8.),
  168. (6., 7.),
  169. (8., 10.),
  170. (6., 12.),
  171. ];
  172. assert_display!(closest_points(&vals), Some(((4., 14.), (4., 15.))));
  173. }
  174. #[test]
  175. fn vertical_points() {
  176. let vals = [
  177. (0., 0.),
  178. (0., 50.),
  179. (0., -25.),
  180. (0., 40.),
  181. (0., 42.),
  182. (0., 100.),
  183. (0., 17.),
  184. (0., 29.),
  185. (0., -50.),
  186. (0., 37.),
  187. (0., 34.),
  188. (0., 8.),
  189. (0., 3.),
  190. (0., 46.),
  191. ];
  192. assert_display!(closest_points(&vals), Some(((0., 40.), (0., 42.))));
  193. }
  194. }