/**
* Closest pair
* Input: N points (x,y)
* Output: A pair of points (x1,y1),(x2,y2) which are closest.
* Algorithm: Line sweep (see http://w...content-available-to-author-only...r.com/tc?module=Static&d1=tutorials&d2=lineSweep)
*
* Coder: Erel Segal http://t...content-available-to-author-only...s.fm/rentabrain
* Created: 2010-12-15
*/
#include <iostream>
#include <algorithm>
#include <set>
#include <map>
#include <vector>
#include <time.h>
#include <math.h>
#include <assert.h>
#include <limits>
using namespace std;
template <class Iterator> ostream& print (ostream& out, Iterator iFrom, Iterator iTo) {
for (; iFrom!=iTo; ++iFrom) out << *iFrom << " ";
return (out << endl);
}
template <class T> ostream& operator<< (ostream& out, const vector<T>& v) {
return print (out, v.begin(), v.end());
}
template <class T> ostream& operator<< (ostream& out, const set<T>& v) {
return print (out, v.begin(), v.end());
}
const int INF = 999999999;
struct Point {
int x, y;
Point() {x=y=INF;}
Point(int _x, int _y) { x=_x; y=_y; }
friend istream& operator>> (istream& in, Point& r) { in >> r.x >> r.y; return in; }
friend ostream& operator<< (ostream& out, const Point& r) { out << r.x << "," << r.y; return out; }
};
bool lefter(const Point& A, const Point& B) { return A.x < B.x; }
bool righter(const Point& A, const Point& B) { return A.x > B.x; }
bool lower(const Point& A, const Point& B) { return A.y < B.y; }
bool upper(const Point& A, const Point& B) { return A.y > B.y; }
float euclidean_distance(const Point& A, const Point& B) { return sqrt( (A.x-B.x)*(A.x-B.x) + (A.y-B.y)*(A.y-B.y) ); }
class PointCollection {
vector<Point> myPoints;
typedef vector<Point>::const_iterator MyPointsIterator;
public:
PointCollection(int N=0): myPoints(N) {}
void clear() { myPoints.clear(); }
void add(const Point& p) {
myPoints.push_back(p);
}
void closestPairBruteForce() const {
Point p1, p2;
float shortestDistance = INF;
for (MyPointsIterator i=myPoints.begin(); i!=myPoints.end(); ++i) {
for (MyPointsIterator j=myPoints.begin(); j!=i; ++j) {
float currentDistance = euclidean_distance(*i, *j);
if (currentDistance<shortestDistance) {
shortestDistance = currentDistance;
p1 = *i;
p2 = *j;
}
}
}
cout << "Brute Force Closest pair: " << p1 << " - " << p2 << " (distance: " << shortestDistance << ")" << endl;
}
/// @see http://w...content-available-to-author-only...r.com/tc?module=Static&d1=tutorials&d2=lineSweep
void closestPairLineSweep() const {
Point p1, p2;
float shortestDistance = INF;
//cout << "update shortest distance to " << shortestDistance << endl;
// A. Order myPoints by the x coordinate:
vector<Point> myPointsOrderedByX(myPoints);
::sort(myPointsOrderedByX.begin(), myPointsOrderedByX.end(), &lefter);
//cout << "myPointsOrderedByX: " << myPointsOrderedByX;
// B. Maintain an active set of points - candidates to be a closest-pair with the current point.
// The set is ordered by the y coordinate:
typedef set<Point, bool(*)(const Point&,const Point&)> ActiveSetType;
typedef ActiveSetType::iterator ActiveSetIterator;
ActiveSetType myActiveSet (&lower);
MyPointsIterator leftestPointInActiveSet = myPointsOrderedByX.begin();
// C. Line sweep the points by X order:
for (MyPointsIterator point=myPointsOrderedByX.begin(); point!=myPointsOrderedByX.end(); ++point) {
// C-1. Scan the active set for points that may be closer to the current point than the shortest distance:
Point lowestPoint (point->x, point->y - shortestDistance);
Point highestPoint (point->x, point->y + shortestDistance);
ActiveSetIterator lowest = myActiveSet.lower_bound(lowestPoint);
if (lowest!=myActiveSet.end()) {
ActiveSetIterator highest = myActiveSet.upper_bound(highestPoint);
//cout << "scan active set between " << lowestPoint << "<=" << *lowest << " and " << *highest << "<" << highestPoint << endl;
for (ActiveSetIterator other=lowest; other!=highest; ++other) {
//cout << " compare " << *point << " with " << *other << endl;
float currentDistance = euclidean_distance(*point,*other);
if (currentDistance<shortestDistance) {
shortestDistance = currentDistance;
//cout << " update shortest distance to " << shortestDistance << endl;
p2 = *other;
p1 = *point;
}
}
}
// C-2. Update the active set - keep only the points whose distance with the current point is at most the shortest distance
//cout << "insert " << *point << " to active set" << endl;
myActiveSet.insert(*point);
while (point->x - leftestPointInActiveSet->x > shortestDistance) {
//cout << " remove " << *leftestPointInActiveSet << " from active set" << endl;
myActiveSet.erase(*leftestPointInActiveSet);
++leftestPointInActiveSet;
}
}
cout << "Line Sweep Closest pair: " << p1 << " - " << p2 << " (distance: " << shortestDistance << ")" << endl;
}
friend ostream& operator<< (ostream& out, PointCollection points) {
return (out << points.myPoints);
}
void fillRandomX(int N) {
Point p;
for (int i=0; i<N; ++i) {
p.x = rand()%1000;
p.y = 500;
add(p);
}
}
void fillRandomY(int N) {
Point p;
for (int i=0; i<N; ++i) {
p.y = rand()%1000;
p.x = 500;
add(p);
}
}
void fillRandomXY(int N) {
Point p;
for (int i=0; i<N; ++i) {
p.y = rand()%10000;
p.x = rand()%10000;
add(p);
}
}
void fillRandomXYWithCheckpoints(int N) {
// In this testcase, the min distance should always be 0, as there
// are two identical points.
Point checkpoint(rand()%1000, rand()%1000);
int i1 = rand()%N;
int i2 = rand()%N;
Point p;
for (int i=0; i<N; ++i) {
if (i==i1 || i==i2) {
add(checkpoint);
} else {
p.y = rand()%1000;
p.x = rand()%1000;
add(p);
}
}
}
};
int main() {
srand(time(NULL));
PointCollection points;
points.clear();
points.fillRandomXY(1000);
points.closestPairLineSweep();
points.closestPairBruteForce();
}