编程之美之寻找最近点对

转载于  http://sxnuwhui.blog.163.com/blog/static/137068373201264104915935/?COLLCC=3097019025&COLLCC=3097019024&

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

    一种简单的想法是暴力枚举每两个点，记录最小距离，显然，时间复杂度为O(n^2)。

    在这里介绍一种时间复杂度为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)。

    在二维情形下，类似的，利用分治法，但是难点在于如何实现线性的合并？

    由上图可见，形成的宽为2d的带状区间，最多可能有n个点，合并时间最坏情况下为n^2,。但是，P1和P2中的点具有以下稀疏的性质，对于P1中的任意一点，P2中的点必定落在一个d X 2d的矩形中，且最多只需检查六个点（鸽巢原理）。

    这样，先将带状区间的点按y坐标排序，然后线性扫描，这样合并的时间复杂度为O(nlogn)，几乎为线性了。

算法描述：已知集合S中有n个点，分治法的思想就是将S进行拆分，分为2部分求最近点对。算法每次选择一条垂线L，将S拆分左右两部分为SL和SR，L一般取点集S中所有点的中间点的x坐标来划分，这样可以保证SL和SR中的点数目各为n/2，

（否则以其他方式划分S，有可能导致SL和SR中点数目一个为1，一个为n-1，不利于算法效率，要尽量保持树的平衡性）

依次找出这两部分中的最小点对距离：δL和δR，记SL和SR中最小点对距离δ = min（δL，δR），如图1：

     以L为中心，δ为半径划分一个长带，最小点对还有可能存在于SL和SR的交界处，如下图2左图中的虚线带，p点和q点分别位于SL和SR的虚线范围内，在这个范围内，p点和q点之间的距离才会小于δ，最小点对计算才有意义。

Figure 2

      对于SL虚框范围内的p点，在SR虚框中与p点距离小于δ的顶多只有六个点，就是图二右图中的2个正方形的6的顶点。这个可以反推证明，如果右边这2个正方形内有7个点与p点距离小于δ，例如q点，则q点与下面正方形的四个顶点距离小于δ，则和δ为SL和SR中的最小点对距离相矛盾。因此对于SL虚框中的p点，不需求出p点和右边虚线框内所有点距离，只需计算SR中与p点y坐标距离最近的6个点，就可以求出最近点对，节省了比较次数。

（否则的话，最坏情形下，在SR虚框中有可能会有n/2个点，对于SL虚框中的p点，每次要比较n/2次，浪费了算法的效率）

HDOJ1007：http://acm.hdu.edu.cn/showproblem.php?pid=1007

代码如下：

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

struct Point
{
    double x;
    double y;
}point[N];
int n;
int tmpt[N];

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)
            tmpt[k++] = i;
    }
    sort(tmpt,tmpt+k,cmpy);
    //线性扫描
    for(i = 0; i < k; i++)
    {
        for(j = i+1; j < k && point[tmpt[j]].y-point[tmpt[i]].y<d; j++)//  理解！！
        {
            double d3 = dis(tmpt[i],tmpt[j]);
            if(d > d3)
                d = d3;
        }
    }
    return d;
}


int main()
{
    while(true)
    {
        scanf("%d",&n);
        if(n==0)
            break;
        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;
}