

 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) 


 
 
type point = { x : float; y : float }
 
 
let cmpPointX (a : point) (b : point) = compare a.x b.x 
let cmpPointY (a : point) (b : point) = compare a.y b.y 
 
 
let distSqrd (seg : (point * point) option) =
  match seg with
  | None -> max_float
  | Some(line) ->
    let a = fst line in
    let b = snd line in
 
    let dx = a.x -. b.x in
    let dy = a.y -. b.y in
 
    dx*.dx +. dy*.dy
 
 
let dist seg =
  sqrt (distSqrd seg)
 
 
let shortest l1 l2 =
  if distSqrd l1 < distSqrd l2 then
    l1
  else
    l2
 
 
let halve l =
  let n = List.length l in
  BatList.split_at (n/2) l
 
 
let rec closestBoundY from maxY (ptsByY : point list) =
  match ptsByY with
  | [] -> None
  | hd :: tl ->
    if hd.y > maxY then
      None
    else
      let toHd = Some(from, hd) in
      let bestToRest = closestBoundY from maxY tl in
      shortest toHd bestToRest
 
 
let rec closestInRange ptsByY maxDy =
  match ptsByY with
  | [] -> None
  | hd :: tl ->
    let fromHd = closestBoundY hd (hd.y +. maxDy) tl in
    let fromRest = closestInRange tl maxDy in
    shortest fromHd fromRest
 
 
let rec closestPairByX (ptsByX : point list) =
   if List.length ptsByX < 2 then
       None
   else
       let (left, right) = halve ptsByX in
       let leftResult = closestPairByX left in
       let rightResult = closestPairByX right in
 
       let bestInHalf = shortest  leftResult rightResult in
       let bestLength = dist bestInHalf in
 
       let divideX = (List.hd right).x in
       let inBand = List.filter(fun(p) -> abs_float(p.x -. divideX) < bestLength) ptsByX in
 
       let byY = List.sort cmpPointY inBand in
       let bestCross = closestInRange byY bestLength in
       shortest bestInHalf bestCross      
 
 
let closestPair pts =
  let ptsByX = List.sort cmpPointX pts in
  closestPairByX ptsByX
 
 
let parsePoint str =
  let sep = Str.regexp_string "," in
  let tokens = Str.split sep str in
  let xStr = List.nth tokens 0 in
  let yStr = List.nth tokens 1 in
 
  let xVal = (float_of_string xStr) in
  let yVal = (float_of_string yStr) in
 
  { x = xVal; y = yVal }
 
 
let loadPoints filename =
  let ic = open_in filename in
  let result = ref [] in
  try
    while true do
      let s = input_line ic in
      if s <> "" then
        let p = parsePoint s in
        result := p :: !result;
    done;
    !result
  with End_of_file ->
    close_in ic;
    !result
;;
 
let loaded = (loadPoints "Points.txt") in
let start = Sys.time() in
let c = closestPair loaded in
let taken = Sys.time() -. start in
Printf.printf "Took %f [s]\n" taken;
 
match c with
| None -> Printf.printf "No closest pair\n"
| Some(seg) -> 
  let a = fst seg in
  let b = snd seg in
 
  Printf.printf "(%f, %f) (%f, %f) Dist %f\n" a.x a.y b.x b.y (dist c)
 
 