

 This page uses content from Wikipedia. The original article was at Closest pair of points problem. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)



Task

Provide a function to find the closest two points among a set of given points in two dimensions,   i.e. to solve the   Closest pair of points problem   in the   planar   case.

The straightforward solution is a   O(n2)   algorithm   (which we can call brute-force algorithm);   the pseudo-code (using indexes) could be simply:

bruteForceClosestPair of P(1), P(2), ... P(N)
if N < 2 then
  return ∞
else
  minDistance ← |P(1) - P(2)|
  minPoints ← { P(1), P(2) }
  foreach i ∈ [1, N-1]
    foreach j ∈ [i+1, N]
      if |P(i) - P(j)| < minDistance then
        minDistance ← |P(i) - P(j)|
        minPoints ← { P(i), P(j) } 
      endif
    endfor
  endfor
  return minDistance, minPoints
 endif

A better algorithm is based on the recursive divide&conquer approach,   as explained also at   Wikipedia's Closest pair of points problem,   which is   O(n log n);   a pseudo-code could be:

closestPair of (xP, yP)
               where xP is P(1) .. P(N) sorted by x coordinate, and
                     yP is P(1) .. P(N) sorted by y coordinate (ascending order)
if N ≤ 3 then
  return closest points of xP using brute-force algorithm
else
  xL ← points of xP from 1 to ⌈N/2⌉
  xR ← points of xP from ⌈N/2⌉+1 to N
  xm ← xP(⌈N/2⌉)x
  yL ← { p ∈ yP : px ≤ xm }
  yR ← { p ∈ yP : px > xm }
  (dL, pairL) ← closestPair of (xL, yL)
  (dR, pairR) ← closestPair of (xR, yR)
  (dmin, pairMin) ← (dR, pairR)
  if dL < dR then
    (dmin, pairMin) ← (dL, pairL)
  endif
  yS ← { p ∈ yP : |xm - px| < dmin }
  nS ← number of points in yS
  (closest, closestPair) ← (dmin, pairMin)
  for i from 1 to nS - 1
    k ← i + 1
    while k ≤ nS and yS(k)y - yS(i)y < dmin
      if |yS(k) - yS(i)| < closest then
        (closest, closestPair) ← (|yS(k) - yS(i)|, {yS(k), yS(i)})
      endif
      k ← k + 1
    endwhile
  endfor
  return closest, closestPair
endif



References and further readings

   Closest pair of points problem
   Closest Pair (McGill)
   Closest Pair (UCSB)
   Closest pair (WUStL)
   Closest pair (IUPUI) 


 
//! 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());
}
 