数字之魅：寻找二维平面上的最近的点对

在二维平面上的n个点中，如何快速的找出最近的一对点，就是最近点对问题。





初看这个题，可能感觉有点儿复杂。
方案一：蛮力法。数组中总共包含Ｎ个数，所以我们可以把平面内所有的点按Ｘ轴排序，然后依次算出后一个坐标与前面所有左边的距离，然后用Ｍin和position来记录最近的距离和两个坐标。该方案和在一维空间求两个最近点的距离有点儿类似，其时间复杂度为：O(N*N).
方案二：在一维空间里，我们知道如果数组有序，我们可以很快找出最近的两个点。我们可以用O(N*logN)的时间复杂度来对数据进行排序[快，堆，归]。排完序后只需要O(N)的时间复杂度就可以得到最小的差值。但是该方法不能用在二维，因为在二维空间中，X轴间距离最小的两个点距离不一定是最短的，如下图所示。


那还有什么方法吗？这里我们采用了分治法。我们利用数组的中位数[中位数的求解细节可以参考这里]，将数组分成Left和Right两个部分，要么来自Right部分，要么来自Left部分，或者来自Left和Right。这里我们就可以将时间复杂度降低到O(N*logN)。

这里主要复杂在合并的地方，其思想如下：

   第一步：首先找出点集数据中的中位数median，按X坐标来划分;用median对点集数据进行划分，左边的为data1，右边的为data2;

   第二步：对分成的两部分分别求出data1和data2中的最近点对，记为：MinDis1和MinDis2；

   第三步：求出data1和data2最近点对距离的较小值：MinDis = min{MinDis1,MinDis2};

   第四步：找出data2中y值前6大的点[利用鸽巢原理，因为其距离不可能小于Min，而我们只考虑在Min*2*Min的方框内的数据]，

对于data1中的点，与data2中的每一个点计算距离MinDisO, 如果MinDisO < MinDis，就改变MinDis的值，MinDis=MinDis0；说明在最小的距离一个来自左边一个来自右边，在合并的时候产生。

参考代码如下：

#include <iostream>
#include <vector>
#include <algorithm>
#include <cmath>
using namespace std;
// 顶点信息
struct Point
{
double m_x, m_y;
Point():m_x(0.0),m_y(0.0) {}
Point(double x, double y):m_x(x),m_y(y){}
bool operator==(const Point& p) const
{return m_x==p.m_x && m_y==p.m_y;}
};

ostream& operator<<(ostream& os, const Point& p)
{
return os << "(" << p.m_x << "," << p.m_y << ")";
}

// 插入排序
template<class T, class Pr>
void insert_sort(vector<T> &vec, int l, int r, Pr pred)
{
int i, j;
for (i=l+1; i<=r; i++)
{
T tmp = vec[i];
for (j=i-1; j>=l && pred(tmp,vec[j]); j--)
vec[j+1]=vec[j];
vec[j+1] = tmp;
}
}

// 找到key所在的位置
template<class T>
int get_position(vector<T> &vec, int l, int r, T key)
{
for (int i=l; i<=r; i++)
if (key == vec[i])
return i;
return -1;
}

// 按第一个元素对vec进行划分
template<class T, class Pr>
int partition(vector<T> &vec, int l, int r, Pr pred)
{
int i, j;
for (i=l+1,j=l; i<=r; i++)
{
if (pred(vec[i],vec[l]))
{
++j;
swap(vec[i],vec[j]);
}
}
swap(vec[j],vec[l]);
return j;
}

// 顺序统计得到第k个元素的值
template<class T, class Pr>
T select(vector<T> &vec, int l, int r, int k, Pr pred)
{
int n = r-l+1;
if (n==1)
{
if (k!=0)
printf("Out of Boundary!\n");
return vec[l];
}
// 找中位数的中位数作为分割点
int cnt = n/5;
int tcnt = (n+4)/5;
int rem = n%5;
vector<T> group(tcnt);
int i, j;
for (i=0,j=l; i<cnt; i++,j+=5)
{
insert_sort(vec, j, j+4, pred);
group[i] = vec[j+2];
}
if (rem)
{
insert_sort(vec, j, j+rem-1, pred);
group[i] = vec[j+(rem-1)/2];
}
T key = select(group, 0, tcnt-1, (tcnt-1)/2, pred);
// 找到分割点的位置
int key_pos = get_position(vec, l, r, key);
swap(vec[key_pos], vec[l]);
// 用分割点对数组进行划分，小的在左边，大的在右边
int pos = partition(vec, l, r, pred);
int x = pos - l;
if (x == k) return key;
else if (x < k)
return select(vec, pos+1, r, k-x-1, pred);
else
return select(vec, l, pos-1, k, pred);
}

// 计算点a和b的距离
double dist(const Point& a, const Point& b)
{
double x = a.m_x-b.m_x;
double y = a.m_y-b.m_y;
return sqrt(x*x+y*y);
}

bool cmpX(const Point& a, const Point& b)
{
return a.m_x < b.m_x;
}

bool cmpY(const Point& a, const Point& b)
{
return a.m_y < b.m_y;
}

double minDifferent(vector<Point> p, int l, int r, vector<Point> &result)
{
// 按中位数进行划分后的子区域的元素个数都会减小到2或3，不会再到1
if ((r-l+1)==2)
{
result[0] = p[l];
result[1] = p[r];
if (cmpX(p[r],p[l])) swap(p[l], p[r]);
return dist(p[l], p[r]);
}
if ((r-l+1)==3)
{
insert_sort(p, l, r, cmpX);
double tmp1 = dist(p[l], p[l+1]);
double tmp2 = dist(p[l+1], p[l+2]);
double ret = min(tmp1, tmp2);
if (tmp1 == ret)
{
result[0] = p[l];
result[1] = p[l+1];
}
else
{
result[0] = p[l+1];
result[1] = p[l+2];
}
return ret;
}
// 大于3个点的情况
int mid = (r+l)>>1;
Point median = select(p, l, r, mid-l, cmpX);
vector<Point> res1(2), res2(2);
double min_l = minDifferent(p, l, mid, res1);
double min_r = minDifferent(p, mid+1, r, res2);
double minum = min(min_l, min_r);
if (minum == min_l)
{
result[0] = res1[0];
result[1] = res1[1];
}
else
{
result[0] = res2[0];
result[1] = res2[1];
}
// 对[p[mid+1]-minum, p[mid]+minum]的带状区域按y排序
vector<Point> yvec;
int i, j;
for (i=mid+1; i<=r; i++)
if (p[i].m_x - p[mid].m_x < minum)
yvec.push_back(Point(p[i]));
for (i=mid; i>=l; i--)
if (p[mid+1].m_x - p[i].m_x < minum)
yvec.push_back(Point(p[i]));
sort(yvec.begin(), yvec.end(), cmpY);
for (i=0; i<yvec.size(); i++)
{
// 至多只有与其后最多7个点的距离会小于minum
for (j=i+1; j<yvec.size() && yvec[j].m_y-yvec[i].m_y<minum &&
j<=i+7; j++)
{
double delta = dist(yvec[i],yvec[j]);
if (delta < minum)
{
minum = delta;
result[0] = yvec[i];
result[1] = yvec[j];
}
}
}
return minum;
}

int main()
{
int n, i, j, x, y;
vector<Point> result(2);
vector<Point> input;
cout<<"please input the number of your data:"<<endl;
cin >> n;
cout<<"please input your data:"<<endl;
for (i=0; i<n; i++)
{
cin >> x;
cin >> y;
input.push_back(Point(x,y));
}
double minum = minDifferent(input, 0, input.size()-1, result);
cout << "nearest point: " << result[0] << " and "
<< result[1] << endl;
cout << "distance: " << minum << endl;
system("pause");
return 0;
}

运行结果：