119 lines
3.9 KiB
Rust
119 lines
3.9 KiB
Rust
//! We interpret complex numbers as points in the Cartesian plane, here. We also use the
|
|
//! [sweepline/plane sweep closest pairs algorithm][algorithm] instead of the divide-and-conquer
|
|
//! algorithm, since it's (arguably) easier to implement, and an efficient implementation does not
|
|
//! require use of unsafe.
|
|
//!
|
|
//! [algorithm]: http://www.cs.mcgill.ca/~cs251/ClosestPair/ClosestPairPS.html
|
|
extern crate num;
|
|
|
|
use num::complex::Complex;
|
|
use std::cmp::{Ordering, PartialOrd};
|
|
use std::collections::BTreeSet;
|
|
type Point = Complex<f32>;
|
|
|
|
/// Wrapper around `Point` (i.e. `Complex<f32>`) so that we can use a `TreeSet`
|
|
#[derive(PartialEq)]
|
|
struct YSortedPoint {
|
|
point: Point,
|
|
}
|
|
|
|
impl PartialOrd for YSortedPoint {
|
|
fn partial_cmp(&self, other: &YSortedPoint) -> Option<Ordering> {
|
|
(self.point.im, self.point.re).partial_cmp(&(other.point.im, other.point.re))
|
|
}
|
|
}
|
|
|
|
impl Ord for YSortedPoint {
|
|
fn cmp(&self, other: &YSortedPoint) -> Ordering {
|
|
self.partial_cmp(other).unwrap()
|
|
}
|
|
}
|
|
|
|
impl Eq for YSortedPoint {}
|
|
|
|
fn closest_pair(points: &mut [Point]) -> Option<(Point, Point)> {
|
|
if points.len() < 2 {
|
|
return None;
|
|
}
|
|
|
|
points.sort_by(|a, b| (a.re, a.im).partial_cmp(&(b.re, b.im)).unwrap());
|
|
|
|
let mut closest_pair = (points[0], points[1]);
|
|
let mut closest_distance_sqr = (points[0] - points[1]).norm_sqr();
|
|
let mut closest_distance = closest_distance_sqr.sqrt();
|
|
|
|
// the strip that we inspect for closest pairs as we sweep right
|
|
let mut strip: BTreeSet<YSortedPoint> = BTreeSet::new();
|
|
strip.insert(YSortedPoint { point: points[0] });
|
|
strip.insert(YSortedPoint { point: points[1] });
|
|
|
|
// index of the leftmost point on the strip (on points)
|
|
let mut leftmost_idx = 0;
|
|
|
|
// Start the sweep!
|
|
for (idx, point) in points.iter().enumerate().skip(2) {
|
|
// Remove all points farther than `closest_distance` away from `point`
|
|
// along the x-axis
|
|
while leftmost_idx < idx {
|
|
let leftmost_point = &points[leftmost_idx];
|
|
if (leftmost_point.re - point.re).powi(2) < closest_distance_sqr {
|
|
break;
|
|
}
|
|
strip.remove(&YSortedPoint {
|
|
point: *leftmost_point,
|
|
});
|
|
leftmost_idx += 1;
|
|
}
|
|
|
|
// Compare to points in bounding box
|
|
{
|
|
let low_bound = YSortedPoint {
|
|
point: Point {
|
|
re: ::std::f32::INFINITY,
|
|
im: point.im - closest_distance,
|
|
},
|
|
};
|
|
let mut strip_iter = strip.iter().skip_while(|&p| p < &low_bound);
|
|
loop {
|
|
let point2 = match strip_iter.next() {
|
|
None => break,
|
|
Some(p) => p.point,
|
|
};
|
|
if point2.im - point.im >= closest_distance {
|
|
// we've reached the end of the box
|
|
break;
|
|
}
|
|
let dist_sqr = (*point - point2).norm_sqr();
|
|
if dist_sqr < closest_distance_sqr {
|
|
closest_pair = (point2, *point);
|
|
closest_distance_sqr = dist_sqr;
|
|
closest_distance = dist_sqr.sqrt();
|
|
}
|
|
}
|
|
}
|
|
|
|
// Insert point into strip
|
|
strip.insert(YSortedPoint { point: *point });
|
|
}
|
|
|
|
Some(closest_pair)
|
|
}
|
|
|
|
pub fn main() {
|
|
let mut test_data = [
|
|
Complex::new(0.654682, 0.925557),
|
|
Complex::new(0.409382, 0.619391),
|
|
Complex::new(0.891663, 0.888594),
|
|
Complex::new(0.716629, 0.996200),
|
|
Complex::new(0.477721, 0.946355),
|
|
Complex::new(0.925092, 0.818220),
|
|
Complex::new(0.624291, 0.142924),
|
|
Complex::new(0.211332, 0.221507),
|
|
Complex::new(0.293786, 0.691701),
|
|
Complex::new(0.839186, 0.728260),
|
|
];
|
|
let (p1, p2) = closest_pair(&mut test_data[..]).unwrap();
|
|
println!("Closest pair: {} and {}", p1, p2);
|
|
println!("Distance: {}", (p1 - p2).norm_sqr().sqrt());
|
|
}
|