平面最近点对问题(分治)

 

平面最近点对问题是指：在给出的同一个平面内的所有点的坐标，然后找出这些点中最近的两个点的距离.

 

方法1：穷举

 

1）算法描述：已知集合S中有n个点，一共可以组成n(n-1)/2对点对，蛮力法就是对这n(n-1)/2对点对逐对进行距离计算，通过循环求得点集中的最近点对
2）算法时间复杂度：算法一共要执行 n(n-1)/2次循环，因此算法复杂度为O(n²)

代码实现：

利用两个for循环可实现所有点的配对,每次配对算出距离然后更新最短距离.

 

for (i=0 ; i < n ;i ++)｛
     for(j= i+1 ; j<n ;j ++)｛
              点i与点j的配对
       ｝
｝

　　 

方法2：分治

 

1) 把它分成两个或多个更小的问题；

 2) 分别解决每个小问题； 

3) 把各小问题的解答组合起来，即可得到原问题的解答。小问题通常与原问题相似，可以递归地使用分而治之策略来解决。

 

在这里介绍一种时间复杂度为O(nlognlogn)的算法。其实，这里用到了分治的思想。将所给平面上n个点的集合S分成两个子集S1和S2，每个子集中约有n/2个点。然后在每个子集中递归地求最接近的点对。在这里，一个关键的问题是如何实现分治法中的合并步骤，即由S1和S2的最接近点对，如何求得原集合S中的最接近点对。如果这两个点分别在S1和S2中，问题就变得复杂了。

为了使问题变得简单，首先考虑一维的情形。此时，S中的n个点退化为x轴上的n个实数x1，x2，...，xn。最接近点对即为这n个实数中相差最小的两个实数。显然可以先将点排好序，然后线性扫描就可以了。但我们为了便于推广到二维的情形，尝试用分治法解决这个问题。

假设我们用m点将S分为S1和S2两个集合，这样一来，对于所有的p(S1中的点)和q(S2中的点)，有p<q。

递归地在S1和S2上找出其最接近点对{p1,p2}和{q1,q2}，并设

　　　　　　　　　　　　　　　　　　d = min{ |p1-p2| , |q1-q2| } 　　　　　　　　　　　　　　　　　　

由此易知，S中最接近点对或者是{p1,p2}，或者是{q1,q2}，或者是某个{q3,p3}，如下图所示。

    如果最接近点对是{q3,p3}，即|p3-q3|<d，则p3和q3两者与m的距离都不超过d，且在区间(m-d,d]和(d,m+d]各有且仅有一个点。这样，就可以在线性时间内实现合并。

    此时，一维情形下的最近点对时间复杂度为O(nlogn)。

 

在二维的情况下:

我们仿照一维的情况先把所有点按照x(横坐标)从左到右升序排列.

以X横坐标中间的点作为分界线.将平面的点分成左边和右边,以上图为例,分为左边5个右边5个.

然后找到左边的点中最近点对的距离d1,和右边最近点对的距离d2。

令d=min{d1,d2}.那么d是当前整个平面的最近点对的距离吗?

答案不是!!

 

当前的d1,d2都是两边各自的点的最近距离,但是没有考虑到两边的点相互配对的情况,即是点对一个在左边区域,一个在右边区域.

如果我们这个时候把两边相互配对的所有情况全部罗列出来那么时间复杂度就变为第一种方法的O(n²).

这个时候我们便可以用上面找到的d来限制配对.即是明显超过d的两点不需要配对.如下：

那么在范围内的只有三个点,也就是说只有这个三个点相互配对的距离才可能比d小.那么再按照方法一一样在这些点中找到最短距离再跟d去比较.小的就是当前整个平面中所考虑的所有点中最近点对的距离。

然而在以上问题处理上,有一个问题尚未解决就是如何找到两边区域中的最近点对的距离d1,d2 ?

我们可以发现在处理上面这个问题的时候，和最开始处理所有平面的点最近点距离的问题类似，只是点的数目减半了.

那么我们可以递归以上问题.直到划分的区域中只有一个或者两个点.这样,该区域的最近点对距离就是无穷大或者该两点的距离。这也是递归终止的条件。

 

代码实现：

#include <iostream>  
#include <cstdio>  
#include <cstring>  
#include <cmath>  
#include <algorithm>  
using namespace std;
const double inf = 1e20;
const int maxn = 100005;

struct Point{
    double x, y;
}point[maxn];

int n, mpt[maxn];


//以x为基准排序
bool cmpxy(const Point& a, const Point& b){
    if (a.x != b.x)
        return a.x < b.x;
    return a.y < b.y;
}

bool cmpy(const int& a, const int& b){
    return point[a].y < point[b].y;
}

double min(double a, double b){
    return a < b ? a : b;
}

double dis(int i, int j){
    return sqrt((point[i].x - point[j].x)*(point[i].x - point[j].x) + (point[i].y - point[j].y)*(point[i].y - point[j].y));
}

double Closest_Pair(int left, int right){
    double d = inf;
    if (left == right)
        return d;
    if (left + 1 == right)
        return dis(left, right);
    int mid = (left + right) >> 1;
    double d1 = Closest_Pair(left, mid);
    double d2 = Closest_Pair(mid + 1, right);
    d = min(d1, d2);
    int i, j, k = 0;
    //分离出宽度为d的区间  
    for (i = left; i <= right; i++){
        if (fabs(point[mid].x - point[i].x) <= d)
            mpt[k++] = i;
    }
    sort(mpt, mpt + k, cmpy);
    //线性扫描  
    for (i = 0; i < k; i++){
        for (j = i + 1; j < k && point[mpt[j]].y - point[mpt[i]].y<d; j++){
            double d3 = dis(mpt[i], mpt[j]);
            if (d > d3)    d = d3;
        }
    }
    return d;
}

int main(){
    while (~scanf("%d", &n) && n){
        for (int i = 0; i < n; i++)
            scanf("%lf %lf", &point[i].x, &point[i].y);
        sort(point, point + n, cmpxy);
        printf("%.2lf\n", Closest_Pair(0, n - 1) / 2);
    }
    return 0;
}