《数据结构与算法分析：C语言描述》复习——第十章“算法设计技巧”——平面最近点对

2014.07.06 17:56

简介：

　　给定二维平面上n个点，求出距离最近的两个点的距离。

描述：

　　个人觉得这本书的第十章是最精华的，因为每种算法思想用一个典型问题来讲解。平面最近点对问题，自然是分治法的例子了。

　　暴力的O(n^2)枚举法自然不用多说，两层循环能得出结果，但效率上只能承担至多几千个点的计算。

　　分治法的基本思路是：如果我将n个点分为两批，那么距离最短的点要么同属一批（这样可以用子问题得到答案），要么各落在一批中（这样就必须合并子问题才能得到答案）。

　　分解为子问题的过程自然很简单，选出一个分界点mid，然后将n个点分为两部分，对每一部分递归求最近距离点对。

　　合并过程复杂一些，并且发生在递归调用结束之后。也就是：

　　　　递归左半部分的子问题 -> 递归右半部分的子问题 -> 合并结果得出最终答案

　　我们定义δ（delta）为左右两部分计算得出的最近距离的较小值。如果最终结果比δ还要小的话，那么它应该会出现在分界线mid处（否则的话，最终结果肯定早就在子问题中得到了，也就不可能比δ还小），而且是由离mid很近的一左一右两个点得到。因此以mid位置的点p[mid]为中心，δ为半径的正方形区域（不是圆形，因为不好计算）都是可能得到最小值的区域。

　　那么，我们从下标mid开始向两边扫描，找出所有落在这个可疑区域的点，于是可以得到左右两个点集P_左和P_右。然后将两个集合中的点一一比对就能得到结果了。

　　如果这样算法就完了，那么实际的程序运行甚至会比暴力算法还慢，因为数据如果一直是无序的，那就无法在该结束循环的时候提前结束循环。

　　排序1：在算法的一开始，就应该选取一个坐标轴（比如X）进行排序，这样我们在“那么，我们从下标mid开始向两边扫描”这句话的时候，就有了break的条件了。如果X坐标是有序的， 那么当某个p[i]的X坐标与p[mid]相差超过了δ，就可以停止扫描了。

　　排序2：“于是可以得到左右两个点集P_左和P_右”，尽管我们可以直接对两个候选集合中的点一一计算距离，得到最短距离。但这么做的话，效率仍然是暴力级别的。这些点全都落在了2δ为边长的正方形区域内，而且这个δ给人的印象很小很小，于是我们可能会认为P_左与P_右加起来点的数量很小，即使枚举也没关系的。在实际的递归调用里，一个看似小的浪费就会被放大到可观察的程度，然后你的代码就超时了。在排序1中X坐标已经有序了，那么可能变化的比较离谱的，只有Y坐标了。于是我们对候选集合中点按Y坐标进行排序。有序的数据，总能给break语句制造机会。这个排序让算法完整了。

　　为了确认算法的效率，我把代码修改了一句话之后提交到了ZOJ的题目“Quoit Design”运行，结果的确能AC掉。

实现：

// Divide-and-conquer solution for Closest Pair of Points problem.
// This piece of code is also accepted on Zhejiang University Online Judge
// Problem ID 2107.
// http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=1107
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <vector>
using namespace std;

struct Point {
    double x;
    double y;
    Point(double _x = 0.0, double _y = 0.0): x(_x), y(_y) {}
};

bool comparatorX(const Point &a, const Point &b)
{
    if (a.x == b.x) {
        return a.y < b.y;
    } else {
        return a.x < b.x;
    }
}

bool comparatorY(const Point &a, const Point &b)
{
    if (a.y == b.y) {
        return a.x < b.x;
    } else {
        return a.y < b.y;
    }
}

double dist(const Point &a, const Point &b)
{
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}

double closestPairOfPointsRecursive(const vector<Point> &points, int left, 
    int right)
{
    if (right - left + 1 == 2) {
        return dist(points[left], points[left + 1]);
    } else if (right - left + 1 == 3) {
        return min(dist(points[left], points[left + 1]), 
            min(dist(points[left + 1], points[left + 2]), 
            dist(points[left], points[left + 2])));
    }
    
    int mid = left + (right - left) / 2;
    double delta = min(closestPairOfPointsRecursive(points, left, mid), 
        closestPairOfPointsRecursive(points, mid + 1, right));
    vector<Point> p;

    int i, j;
    i = mid;
    while (i >= left && points[i].x > points[mid].x - delta) {
        p.push_back(points[i]);
        --i;
    }
    
    i = mid + 1;
    while (i <= right && points[i].x < points[mid].x + delta) {
        p.push_back(points[i]);
        ++i;
    }

    int np;
    double result = delta;

    sort(p.begin(), p.end(), comparatorY);
    np = (int)p.size();
    for (i = 0; i < np; ++i) {
        for (j = i + 1; j < np; ++j) {
            if (p[j].y - p[i].y > result) {
                break;
            }
            result = min(result, dist(p[i], p[j]));
        }
    }
    p.clear();
    
    return result;
}

double closestPairOfPoints(vector<Point> &points)
{
    int n = points.size();
    sort(points.begin(), points.end(), comparatorX);
    return closestPairOfPointsRecursive(points, 0, n - 1);
}

int main()
{
    vector<Point> points;
    int n;
    int i;
    
    while (scanf("%d", &n) == 1 && n > 0) {
        points.resize(n);
        for (i = 0; i < n; ++i) {
            scanf("%lf%lf", &points[i].x, &points[i].y);
        }
        // printf("%.2f
", closestPairOfPoints(points) / 2.0);
        printf("%.2f
", closestPairOfPoints(points));
    }
    
    return 0;
}