《算法分析与设计》之Closest Pair of Points Problem（最近点对问题）d

       今天的算法课上主要讲了最近点对的问题，在老师的讲解下对这个问题有了一个基础的认识和了解，这里就先对这个问题做一个简单的总结吧。

最近点对问题介绍：

       最近点对问题说来其实很简单，主要就是在二维平面内的n个点中，找出（欧式）距离最近的两个点来。

问题解决思路：

       解决这个问题最直接的方法就是“蛮力法”（Brute Force）了，其实现过程即为：

1. 计算出二维平面内所有点对（任意两点的组合）之间的距离
2. 逐一比较这些距离的大小，找出其中的最短的距离，则这两点即为所求的“最近点对”

       当然，蛮力法虽然很容易想到，理解和记忆起来也不困难，在平面内n个点的数目较少时不失为一种解决该问题的简单方法，但是当问题的规模扩大时，这个方法所需要的大量计算和比较，将耗去相当大的时间和资源。

况且，理想环境下的“小数据”情况在现实生活中几乎不存在，面对大量的数据时，找到一个更为简单和快速的方法尤为重要。

这里，就要引入解决这个问题所需要的一种算法设计思想——分治策略（Divide and Conquer）

分治策略简介：

       这里简单介绍一下分治策略的思想：

       分治策略求解问题的过程一般可以分为如下三个步骤：

• 分解问题（Divide）：

       这个步骤将一个问题划分成若干个子问题（注：这些子问题的规模最好是差不多的），这些子问题在形式上和原问题是一样的，只是规模相对来说更小，方便解决。

       例如在我国，要对全国总计约13亿的人口进行一次普查是一个相当复杂且麻烦的问题，但是如果将这个问题的规模由“全国”这一范围依次划分成对各省级单位、市级单位、县区级……进行人口普查，那么，同样是人口普查，但是难度上就小了许多了。所以对问题进行分解，往往是解决一个问题的开始。

• 解决问题（Conquer）：

       这个步骤需要做的工作是“递归地求解出子问题”。这个步骤是要和上一步骤中的分解问题共同进行的。

       还是上面的那个例子，当我们接到任务说要对全国人口进行一次普查时，我们首先将这一问题分解成对全国各个省份的人口进行普查，然后尝试着来解决；可是我们会发现，即使是按省份来普查人口，这个范围还是太大了，于是我们只好继续分解问题，按城市、按县区级来划分都还是没有将问题的规模缩小到足以让我们在短期内能够普查完这种划分下的所有人口，那么说明问题还是没有被分解到“足够小”，所以我们还需要将这种“递归思想”继续下去，直至我们将人口普查的范围缩小到了街道或者村级这样的范围，我们差不多就可以进行这一范围内的人口普查工作了。当然，问题规模也不要太小，例如小到一个家庭的这种程度，这样统计起来虽然快，但未免有些太浪费资源了……

       上面的例子说了这么多废话，其实主要想表达的就是“递归求解子问题，如果子问题的规模足够小，则停止递归，直接求解。”（出自《算法导论（原书第3版）》第四章 分治策略 P37）

• 合并问题（Combine）

       虽然这个步骤并没有在分治策略的名字中有所体现，但是这个步骤却也是不容忽视的。因为小问题的解决并不意味着大问题的解决，就像上面例子中所提到的，我们在第二个步骤完成后得到的还只是全国各个街道或者村的人口普查结果，而不是我们最终想要得到的全国人口普查结果。

       所以要想真正解决原问题，还得把得到的小问题（即子问题）的解依次合并成大问题（即原问题）的解才行。

       好了，到这里为止，分治策略思想的基本情况就介绍完了，下面就开始利用这个思想来解决今天的“最近点对问题”。

       根据上面对分治策略思想的介绍，我们要解决最近点对问题的第一步就要从“分解问题”开始。那么这个问题要怎么分解呢？

       对于平面内的n个点，要怎样将这个平面划分成更小的平面呢？如果是按照我们学习数学时划分象限的方法将平面划分成上下左右四个小平面的话就会是这样的：（注：以下图片均来自老师授课时所用的PPT课件）

       上面这个划分看起来好像没什么问题，我们的确是按照分治策略把原来的平面分成了规模差不多的四个小平面，但需要注意的是，划分平面时，我们需要考虑平面内的点的一些特殊分布情况，像是下面这样的：

       这个时候四个小问题的规模就完全不一样了，就像是工厂把生产任务分给四个员工，分配任务的方式这么不平均，必然会招致员工的不满了。

所以这个时候我们就要考虑，是不是太心急了？一下子就把问题划分成了四个小问题是不是跨度有些太大了？不如把解决问题的急切心情舒缓一下，先把原平面划分成像是下面这样的两个平面看看吧：

       这个时候我们就只要分别求出上图中左右两个子平面中的距离最近的两个点之间的距离就可以了，如下图所示：

       当然，这是在子平面内点数较小的情况下才可以直接求解最近距离的点对的。如果子平面内点数仍然较多，那么还需要递归地对左右两个子平面继续划分下去，直至可以较快求得最近距离点对为止。

       如上图中所计算得到的，左平面中计算得到最近距离为12（单位就暂且忽略吧……），右平面中计算所得的最近距离为21，这样一比较，12明显小于21啊，那是不是说明这个平面内的最近点距就是12了呢？

       我的回答是：不一定。细心的网友一定会发现，在左右平面的分割线附近，还有两个相距非常近的点的存在，因为分割线而使得这两个点分别位于左右两个平面中了，但是在原平面中，它们则有可能才是相距最近的两个点，验证一下吧：

       计算出来了，果然，这两个点之间的距离只有8，比之前求得的最近距离12还要小。

       那么问题来了，要怎么样做才能避免此类情况的发生？先别急着看下面的答案，试着自己想一想。解决的方法其实很简单，我们只要想一下问题是怎么产生的，就能明白如何解决了。

解决思路：

       由于正是两个平面的分割线反而将平面内最近的两点给分割开来了，那么我们只要想办法将这两个点再次合并到一个平面内就可以了。可是如果把左右两个平面再次合并成一个平面，然后计算平面内两点间最近的距离的话，这不又是变成了“蛮力法”的求解方法了吗？

       是的，如果直接再次合并两个平面并计算平面内所有点对距离，这对解决这个问题在速度上提升上将会没有任何效果，反而让前面的步骤变得多余。

       这时，细心的网友（ =(￣ω￣)=细心的网友好厉害……）会再次发现，其实只要计算分割线附近的点的距离就可以了。

       好的，问题又来了，这个“附近”到底是怎么定义的呢？多近才能算是附近？

       再次陷入困惑的网友请不要忽略了我们前面辛辛苦苦计算出来的左右两个平面中的点的最近距离！

       是的，既然要在分割线附近找比左右两个平面内计算出的最近距离还要近的距离，那么范围自然是不能超过目前已知的最近距离了。

       所以，分割线“附近”的定义如下图所示：

       分割线附近的范围定义好了以后，我们就可以通过计算落在这个范围内的点之间的距离并找出它们的最近距离了。

       最后，再次与之前已经求得的最近距离进行比较就可以得到整个平面内的最近点对了！

       至此，今天关于最近点对问题的讲解就结束了，后面有时间会把相应的算法再补充上来。（2015/04/01）

PS:这个日期说了这样的话，总觉得以后不补充上来也能有很好的借口呐(～￣▽￣)～。

（2015/04/15 补充）

时隔半个月（时间过得好快……）来补充一下当初说好的算法……

特别声明：这个算法来自网络（根据搜索及相关链接，该程序作者为CHEN, Yue 如信息有误，欢迎指正）

/**************求最近点对的分治算法实现,输入点对，输出最近点对的距离**************/#include <iostream>#include<cmath>#include <cstdio>#include <cstdlib>#include <cstring>using namespace std;const int N = 100005;const double MAX = 10e100;const double eps = 0.00001;typedef struct TYPE{    double x, y;    int index;} Point;Point a[N], b[N], c[N];double closest(Point *, Point *, Point *, int, int);double dis(Point, Point);int cmp_x(const void *, const void*);int cmp_y(const void *, const void*);int merge(Point *, Point *, int, int, int);inline double min(double, double);int main(){    int n, i;    double d;    scanf("%d", &n);    while (n)    {        for (i = 0; i < n; i++)            scanf("%lf%lf", &(a[i].x), &(a[i].y));        qsort(a, n, sizeof(a[0]), cmp_x);        for (i = 0; i < n; i++)            a[i].index = i;        memcpy(b, a, n *sizeof(a[0]));        qsort(b, n, sizeof(b[0]), cmp_y);        d = closest(a, b, c, 0, n - 1);        printf("%.2lf\n", d);        scanf("%d", &n);    }    return 0;}double closest(Point a[], Point b[], Point c[], int p, int q){    if (q - p == 1)        return dis(a[p], a[q]);    if (q - p == 2)    {        double x1 = dis(a[p], a[q]);        double x2 = dis(a[p + 1], a[q]);        double x3 = dis(a[p], a[p + 1]);        if (x1 < x2 && x1 < x3)            return x1;        else if (x2 < x3)            return x2;        else            return x3;    }    int m = (p + q) / 2;    int i, j, k;    double d1, d2;    for (i = p, j = p, k = m + 1; i <= q; i++)        if (b[i].index <= m)            c[j++] = b[i];    //数组c左半部保存划分后左部的点, 且对y是有序的.    else        c[k++] = b[i];    d1 = closest(a, c, b, p, m);    d2 = closest(a, c, b, m + 1, q);    double dm = min(d1, d2);    merge(b, c, p, m, q); //数组c左右部分分别是对y坐标有序的, 将其合并到b.    for (i = p, k = p; i <= q; i++)        if (fabs(b[i].x - b[m].x) < dm)            c[k++] = b[i];    //找出离划分基准左右不超过dm的部分, 且仍然对y坐标有序.    for (i = p; i < k; i++)    for (j = i + 1; j < k && c[j].y - c[i].y < dm; j++)    {        double temp = dis(c[i], c[j]);        if (temp < dm)            dm = temp;    }    return dm;}double dis(Point p, Point q){    double x1 = p.x - q.x, y1 = p.y - q.y;    return sqrt(x1 *x1 + y1 * y1);}int merge(Point p[], Point q[], int s, int m, int t){    int i, j, k;    for (i = s, j = m + 1, k = s; i <= m && j <= t;)    {        if (q[i].y > q[j].y)            p[k++] = q[j], j++;        else            p[k++] = q[i], i++;    }    while (i <= m)        p[k++] = q[i++];    while (j <= t)        p[k++] = q[j++];    memcpy(q + s, p + s, (t - s + 1) *sizeof(p[0]));    return 0;}int cmp_x(const void *p, const void *q){    double temp = ((Point*)p)->x - ((Point*)q)->x;    if (temp > 0)        return 1;    else if (fabs(temp) < eps)        return 0;    else        return - 1;}int cmp_y(const void *p, const void *q){    double temp = ((Point*)p)->y - ((Point*)q)->y;    if (temp > 0)        return 1;    else if (fabs(temp) < eps)        return 0;    else        return - 1;}inline double min(double p, double q){    return (p > q) ? (q): (p);}

       非常感谢中南大学郑瑾老师的授课，把这个问题讲解的非常通透，以上讲解都是按照老师上课时候的思路加上个人的理解整理出来的，如果我对这个问题的介绍存在任何疏漏，欢迎批评指正！

版权声明：本文完全原创，所引用到的文字、代码及图片均已标明出处，转载请注明出处，谢谢合作！