分治——最近点对问题

利用分治方法的经典问题——最近点对问题（Closest pair of points problem）

问题描述
n个点在公共空间中，求出所有点对的欧几里得距离最小的点对。

问题分析
该题直观的解决方法便是Brute Force(暴力求解)。

minDist = infinity
for i = 1 to length(P) - 1
 for j = i + 1 to length(P)
  let p = P[i], q = P[j]
  if dist(p, q) < minDist:
   minDist = dist(p, q)
   closestPair = (p, q)
return closestPair

利用分治思想进行求解。首先分析题目，符合分治法的适用条件，规模越小容易求解，同时具有最优子结构。
分治法求解
     分解
          对所有的点按照x坐标（或者y）从小到大排序（排序方法时间复杂度O(nlogn)。
          根据下标进行分割，使得点集分为两个集合。
     解决
          递归的寻找两个集合中的最近点对。
          取两个集合最近点对中的最小值min(disleft，disright) min(dis_{left}， dis_{right})min(dis left，disright)。
     合并
          最近距离不一定存在于两个集合中，可能一个点在集合A，一个点在集合B，而这两点间距离小于dis。

1998年，由周玉林、熊鹏荣、朱洪教授提出了平面最近点对的一个改进算法，针对Preparata-Shamos算法提出的6个点，又证明其实只需要4个点就可以确定最近点对距离，该证明提出2个定理，利用更加准确的半径画圈，证明了只要对左半域上的每个点p，检验右半域y坐标与p最近的至多4个点即可（上下个两个）。具体证明可以参考《求平面点集最近点对的一个改进算法》。
根据以上的优化，可以在合并时，通过检测与左半域点p的y坐标相邻的2个或者3个，即使用4点或者6点来检测，一般为了省事，只求与p点y坐标上界或者下界右半域连续的6个、4个点即可。

时间复杂度

代码实现
struct point {
    double x;
    double y;
    point(double x, double y) :x(x), y(y) {}
    point() { return; }
};

bool cmp_x(const point & A, const point & B)  // 比较x坐标
{
    return A.x < B.x;
}

bool cmp_y(const point & A, const point & B)  // 比较y坐标
{
    return A.y < B.y;
}

double distance(const point & A, const point & B)
{
    return sqrt(pow(A.x - B.x, 2) + pow(A.y - B.y, 2));
}
/*
* function: 合并，同第三区域最近点距离比较
* param: points 点的集合
*        dis 左右两边集合的最近点距离
*        mid x坐标排序后，点集合中中间点的索引值
*/
double merge(vector<point> & points, double dis, int mid)
{
    vector<point> left, right;
    for (int i = 0; i < points.size(); ++i)  // 搜集左右两边符合条件的点
    {
        if (points[i].x - points[mid].x <= 0 && points[i].x - points[mid].x > -dis)
            left.push_back(points[i]);
        else if (points[i].x - points[mid].x > 0 && points[i].x - points[mid].x < dis)
            right.push_back(points[i]);
    }
    sort(right.begin(), right.end(), cmp_y);
    for (int i = 0, index; i < left.size(); ++i)  // 遍历左边的点集合，与右边符合条件的计算距离
    {
        for (index = 0; index < right.size() && left[i].y < right[index].y; ++index);
        for (int j = 0; j < 7 && index + j < right.size(); ++j)  // 遍历右边坐标上界的6个点
        {
            if (distance(left[i], right[j + index]) < dis)
                dis = distance(left[i], right[j + index]);
        }
    }
    return dis;
}

double closest(vector<point> & points)
{
    if (points.size() == 2) return distance(points[0], points[1]);  // 两个点
    if (points.size() == 3) return min(distance(points[0], points[1]), min(distance(points[0], points[2]), 
        distance(points[1], points[2])));  // 三个点
    int mid = (points.size() >> 1) - 1;
    double d1, d2, d;
    vector<point> left(mid + 1), right(points.size() - mid - 1);
    copy(points.begin(), points.begin() + mid + 1, left.begin());  // 左边区域点集合
    copy(points.begin() + mid + 1, points.end(), right.begin());  // 右边区域点集合
    d1 = closest(left);
    d2 = closest(right);
    d = min(d1, d2);
    return merge(points, d, mid);
}

int main()
{
    int count;
    printf("点个数：");
    scanf("%d", &count);
    vector<point> points;
    double x, y;
    for (int i = 0; i < count; ++i)
    {
        printf("第%d个点", i);
        scanf("%lf%lf", &x, &y);
        point p(x, y);
        points.push_back(p);
    }
    sort(points.begin(), points.end(), cmp_x);
    printf("最近点对值：%lf", closest(points));
    return 0;
}