#include <algorithm>
#include <vector>
#include <iostream>
#include <math.h>
#include <Windows.h>
using namespace std;
// point class representing point in 2D plane
struct Point
{
int x;
int y;
Point() = default;
Point(int xx, int yy) : x(xx), y(yy) {}
};
// compare function using x coordinate
bool comparX(const Point &p1, const Point &p2)
{
return p1.x < p2.x;
}
// compare function using y coordinate
bool comparY(const Point &p1, const Point &p2)
{
return p1.y < p2.y;
}
// euclidean distance function
float dist(const Point &p1, const Point &p2)
{
return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}
double bruteForce(vector<Point> &points, int start, int end)
{
double minDist = DBL_MAX;
for (int i = start; i < end; ++i)
{
for (int j = i + 1; j < end; ++j)
{
if (dist(points[i], points[j]) < minDist)
minDist = dist(points[i], points[j]);
}
}
return minDist;
}
// each recursive function takes two arrays, one sorted by x coordinate, one by y coordinate
// so the in each recursive function, no more need to sort the strip.
// time complexity is reduced to O(nlog(n))
double closestUtil(vector<Point> &points_x, vector<Point> &points_y, int start, int end)
{
int n = end - start;
// base case: if there are 2 or 3 points, use brute force
if (n <= 3)
return bruteForce(points_x, start, end);
int mid = n / 2; // the mid point
vector<Point> points_y_left;
vector<Point> points_y_right;
for (int i = start; i < end; ++i)
{
if (points_y[i].x < points_x[mid].x)
points_y_left.push_back(points_y[i]);
if (points_y[i].x > points_x[mid].x)
points_y_right.push_back(points_y[i]);
}
double delta1 = closestUtil(points_x, points_y_left, start, mid); // left part
double delta2 = closestUtil(points_x, points_y_right, mid + 1, end); // right part
double delta = min(delta1, delta2);
// strip is the points set that all points are closer than delta to the
// line passing through the mid line
vector<Point> strip;
for (int i = start; i < end; ++i)
if (abs(points_y[i].x - points_x[mid].x) < delta)
strip.push_back(points_y[i]);
// O(n) complexity to find the minimum in the strip
for (int i = 0; i < strip.size(); ++i)
for (int j = i + 1; j < strip.size() && strip[j].y - strip[i].y < delta; ++j)
if (dist(strip[i], strip[j]) < delta)
delta = dist(strip[i], strip[j]);
return delta;
}
double closestPair(vector<Point> &points)
{
vector<Point> points_y(points);
sort(points.begin(), points.end(), comparX); // points sorted by x coordinate
sort(points_y.begin(), points_y.end(), comparY); // points sorted by y coordinate
return closestUtil(points, points_y, 0, points.size());
}
int main()
{
vector<Point> p{ { 2, 3 },{ 12, 30 },{ 40, 50 },{ 5, 1 },{ 12, 10 },{ 3, 4 },{ 100, 3 },{ 10, 10 } };
cout << closestPair(p) << endl;
system("PAUSE");
return 0;
}
reference:
http://www.cs.princeton.edu/~wayne/kleinberg-tardos/pdf/05DivideAndConquerI.pdf
http://www.geeksforgeeks.org/closest-pair-of-points-onlogn-implementation/