最近点算法

1. type Point = (f64, f64);
2. use std::cmp::Ordering;

4. fn point_cmp((a1, a2): &Point, (b1, b2): &Point) -> Ordering {
5. let acmp = f64_cmp(a1, b1);
6. match acmp {
7. Ordering::Equal => f64_cmp(a2, b2),
8. _ => acmp,
9. }
10. }

12. fn f64_cmp(a: &f64, b: &f64) -> Ordering {
13. a.partial_cmp(b).unwrap()
14. }

16. /// returns the two closest points
17. /// or None if there are zero or one point
18. pub fn closest_points(points: &[Point]) -> Option<(Point, Point)> {
19. let mut points: Vec<Point> = points.to_vec();
20. points.sort_by(point_cmp);

22. closest_points_aux(&points, 0, points.len())
23. }

25. fn dist((x1, y1): &Point, (x2, y2): &Point) -> f64 {
26. let dx = *x1 - *x2;
27. let dy = *y1 - *y2;

29. (dx * dx + dy * dy).sqrt()
30. }

32. fn closest_points_aux(
33. points: &[Point],
34. mut start: usize,
35. mut end: usize,
36. ) -> Option<(Point, Point)> {
37. let n = end - start;

39. if n <= 1 {
40. return None;
41. }

43. if n <= 3 {
44. // bruteforce
45. let mut min = dist(&points[0], &points[1]);
46. let mut pair = (points[0], points[1]);

48. for i in 0..n {
49. for j in (i + 1)..n {
50. let new = dist(&points[i], &points[j]);
51. if new < min {
52. min = new;
53. pair = (points[i], points[j]);
54. }
55. }
56. }
57. return Some(pair);
58. }

60. let mid = (start + end) / 2;
61. let left = closest_points_aux(points, start, mid);
62. let right = closest_points_aux(points, mid, end);

64. let (mut min_dist, mut pair) = match (left, right) {
65. (Some((l1, l2)), Some((r1, r2))) => {
66. let dl = dist(&l1, &l2);
67. let dr = dist(&r1, &r2);
68. if dl < dr {
69. (dl, (l1, l2))
70. } else {
71. (dr, (r1, r2))
72. }
73. }
74. (Some((a, b)), None) => (dist(&a, &b), (a, b)),
75. (None, Some((a, b))) => (dist(&a, &b), (a, b)),
76. (None, None) => unreachable!(),
77. };

79. let mid_x = points[mid].0;
80. while points[start].0 < mid_x - min_dist {
81. start += 1;
82. }
83. while points[end - 1].0 > mid_x + min_dist {
84. end -= 1;
85. }

87. let mut mids: Vec<&Point> = points[start..end].iter().collect();
88. mids.sort_by(|a, b| f64_cmp(&a.1, &b.1));

90. for (i, e) in mids.iter().enumerate() {
91. for k in 1..8 {
92. if i + k >= mids.len() {
93. break;
94. }

96. let new = dist(e, mids[i + k]);
97. if new < min_dist {
98. min_dist = new;
99. pair = (**e, *mids[i + k]);
100. }
101. }
102. }

104. Some(pair)
105. }

107. #[cfg(test)]
108. mod tests {
109. use super::closest_points;
110. use super::Point;

112. fn eq(p1: Option<(Point, Point)>, p2: Option<(Point, Point)>) -> bool {
113. match (p1, p2) {
114. (None, None) => true,
115. (Some((p1, p2)), Some((p3, p4))) => (p1 == p3 && p2 == p4) || (p1 == p4 && p2 == p3),
116. _ => false,
117. }
118. }

120. macro_rules! assert_display {
121. ($left: expr, $right: expr) => {
122. assert!(
123. eq($left, $right),
124. "assertion failed: `(left == right)`\nleft: `{:?}`,\nright: `{:?}`",
125. $left,
126. $right
127. )
128. };
129. }

131. #[test]
132. fn zero_points() {
133. let vals: [Point; 0] = [];
134. assert_display!(closest_points(&vals), None::<(Point, Point)>);
135. }

137. #[test]
138. fn one_points() {
139. let vals = [(0., 0.)];
140. assert_display!(closest_points(&vals), None::<(Point, Point)>);
141. }

143. #[test]
144. fn two_points() {
145. let vals = [(0., 0.), (1., 1.)];
146. assert_display!(closest_points(&vals), Some(((0., 0.), (1., 1.))));
147. }

149. #[test]
150. fn three_points() {
151. let vals = [(0., 0.), (1., 1.), (3., 3.)];
152. assert_display!(closest_points(&vals), Some(((0., 0.), (1., 1.))));
153. }

155. #[test]
156. fn list_1() {
157. let vals = [
158. (0., 0.),
159. (2., 1.),
160. (5., 2.),
161. (2., 3.),
162. (4., 0.),
163. (0., 4.),
164. (5., 6.),
165. (4., 4.),
166. (7., 3.),
167. (-1., 2.),
168. (2., 6.),
169. ];
170. assert_display!(closest_points(&vals), Some(((2., 1.), (2., 3.))));
171. }

173. #[test]
174. fn list_2() {
175. let vals = [
176. (1., 3.),
177. (4., 6.),
178. (8., 8.),
179. (7., 5.),
180. (5., 3.),
181. (10., 3.),
182. (7., 1.),
183. (8., 3.),
184. (4., 9.),
185. (4., 12.),
186. (4., 15.),
187. (7., 14.),
188. (8., 12.),
189. (6., 10.),
190. (4., 14.),
191. (2., 7.),
192. (3., 8.),
193. (5., 8.),
194. (6., 7.),
195. (8., 10.),
196. (6., 12.),
197. ];
198. assert_display!(closest_points(&vals), Some(((4., 14.), (4., 15.))));
199. }

201. #[test]
202. fn vertical_points() {
203. let vals = [
204. (0., 0.),
205. (0., 50.),
206. (0., -25.),
207. (0., 40.),
208. (0., 42.),
209. (0., 100.),
210. (0., 17.),
211. (0., 29.),
212. (0., -50.),
213. (0., 37.),
214. (0., 34.),
215. (0., 8.),
216. (0., 3.),
217. (0., 46.),
218. ];
219. assert_display!(closest_points(&vals), Some(((0., 40.), (0., 42.))));
220. }
221. }