C++: Closest-Pair Problem

Posted on September 1, 2017April 1, 2018 by TFE Times

Bjarne-stroustrup

The aim of this task is to 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(n²) algorithm (which we can call brute-force algorithm); the pseudocode (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, which is O(n log n); a pseudocode 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 : p_x ≤ xm }
  yR ← { p ∈ yP : p_x > 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 - p_x| < 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

/*
	Author: Kevin Bacon
	Date: 04/03/2014
	Task: Closest-pair problem
*/

#include <iostream>
#include <vector>
#include <utility>
#include <cmath>
#include <random>
#include <chrono>
#include <algorithm>
#include <iterator>

typedef std::pair<double, double> point_t;
typedef std::pair<point_t, point_t> points_t;

double distance_between(const point_t& a, const point_t& b) {
	return std::sqrt(std::pow(b.first - a.first, 2)
	+ std::pow(b.second - a.second, 2));
}

std::pair<double, points_t> find_closest_brute(const std::vector<point_t>& points) {
	if (points.size() < 2) {
		return { -1, { { 0, 0 }, { 0, 0 } } };
	}
	auto minDistance = std::abs(distance_between(points.at(0), points.at(1)));
	points_t minPoints = { points.at(0), points.at(1) };
	for (auto i = std::begin(points); i != (std::end(points) - 1); ++i) {
		for (auto j = i + 1; j < std::end(points); ++j) {
			auto newDistance = std::abs(distance_between(*i, *j));
			if (newDistance < minDistance) {
				minDistance = newDistance;
				minPoints.first = *i;
				minPoints.second = *j;
			}
		}
	}
	return { minDistance, minPoints };
}

std::pair<double, points_t> find_closest_optimized(const std::vector<point_t>& xP,
const std::vector<point_t>& yP) {
	if (xP.size() <= 3) {
		return find_closest_brute(xP);
	}
	auto N = xP.size();
	auto xL = std::vector<point_t>();
	auto xR = std::vector<point_t>();
	std::copy(std::begin(xP), std::begin(xP) + (N / 2), std::back_inserter(xL));
	std::copy(std::begin(xP) + (N / 2), std::end(xP), std::back_inserter(xR));
	auto xM = xP.at(N / 2).first;
	auto yL = std::vector<point_t>();
	auto yR = std::vector<point_t>();
	std::copy_if(std::begin(yP), std::end(yP), std::back_inserter(yL), [&xM](const point_t& p) {
				return p.first <= xM;
			});
	std::copy_if(std::begin(yP), std::end(yP), std::back_inserter(yR), [&xM](const point_t& p) {
				return p.first > xM;
			});
	auto p1 = find_closest_optimized(xL, yL);
	auto p2 = find_closest_optimized(xR, yR);
	auto minPair = (p1.first <= p2.first) ? p1 : p2;
	auto yS = std::vector<point_t>();
	std::copy_if(std::begin(yP), std::end(yP), std::back_inserter(yS), [&minPair, &xM](const point_t& p) {
				return std::abs(xM - p.first) < minPair.first;
			});
	auto result = minPair;
	for (auto i = std::begin(yS); i != (std::end(yS) - 1); ++i) {
		for (auto k = i + 1; k != std::end(yS) &&
		((k->second - i->second) < minPair.first); ++k) {
			auto newDistance = std::abs(distance_between(*k, *i));
			if (newDistance < result.first) {
				result = { newDistance, { *k, *i } };
			}
		}
	}
	return result;
}

void print_point(const point_t& point) {
	std::cout << "(" << point.first
	<< ", " << point.second
	<< ")";
}

int main(int argc, char * argv[]) {
	std::default_random_engine re(std::chrono::system_clock::to_time_t(
	std::chrono::system_clock::now()));
	std::uniform_real_distribution<double> urd(-500.0, 500.0);
	std::vector<point_t> points(100);
	std::generate(std::begin(points), std::end(points), [&urd, &re]() {
		return point_t { 1000 + urd(re), 1000 + urd(re) };
	});
	auto answer = find_closest_brute(points);
	std::sort(std::begin(points), std::end(points), [](const point_t& a, const point_t& b) {
		return a.first < b.first;
	});
	auto xP = points;
	std::sort(std::begin(points), std::end(points), [](const point_t& a, const point_t& b) {
		return a.second < b.second;
	});
	auto yP = points;
	std::cout << "Min distance (brute): " << answer.first << " ";
	print_point(answer.second.first);
	std::cout << ", ";
	print_point(answer.second.second);
	answer = find_closest_optimized(xP, yP);
	std::cout << "\nMin distance (optimized): " << answer.first << " ";
	print_point(answer.second.first);
	std::cout << ", ";
	print_point(answer.second.second);
	return 0;
}

Output:

Min distance (brute): 6.95886 (932.735, 1002.7), (939.216, 1000.17)
Min distance (optimized): 6.95886 (932.735, 1002.7), (939.216, 1000.17)

Content is available under GNU Free Documentation License 1.2.

Related