算法分析与设计——分治法最近点对

分治法最近点对

分治法

分治法将一个难以直接解决的大问题划分成一些规模较小的子问题，分别求解各个子问题，再合并子问题的解得到原问题的解。

一般来说，分治法的求解过程由以下三个阶段组成：

1. 划分：把规模为n的原问题划分为k个规模较小的子问题。
2. 求解子问题：各个子问题的解法与原问题的解法通常是相同的，可以用递归的方法求解各个子问题，有时递归也可以用循环来实现。
3. 合并：把各个子问题的解合并起来，合并的代价因情况不同有很大的差异，分治算法的效率很大程度上依赖于合并的实现。

最近点对问题

问题描述

设p₁=(x₁,y₁)，p₂=(x₂,y₂)，···，p_n=(x_n,y_n)是平面上n个点构成的集合S，最近点对问题就是找出集合S中距离最近的点对。严格地说，最接近的点对可能多于一个，简单起见，找出一个即可。

应用实例

假设在一个金属片上钻n个大小一样的洞，如果洞的距离太近，金属就可能会断裂。可以通过任意两个洞的最小距离来估算金属断裂的概率。这种最小距离问题实际上也就是最近点对问题。

想法

最近点对的分治策略如下：

（1）划分：将集合S分成两个子集S₁和S₂根据平衡子问题原则，每个子集中大约有n/2个点，设集合S的最近对是p_i和p_j，则会出现以下三种情况：

1. p_i∈S₁，p_j∈S₁，即最近对均在集合S₁中；
2. p_i∈S₂，p_j∈S₂，即最近对均在集合S₂中；
3. p_i∈S₁，p_j∈S₂，即最近对分别在集合S₁和S₂中。

（2）求解子问题：对于划分阶段的情况1和情况2可递归求解，如果最近对分别在两个子集中，在后边做讨论。

（3）合并：比较在划分阶段三种情况下的最近点对，取三者之中距离较小者为原问题的解。

下面讨论划分阶段的情况3：

为了将平面上的点集S分割为点的个数大致相同的两个集合，选取垂直线x=m来作分割线，其中m为S中各点x坐标的中位数。由此将S分割为S₁和S₂。递归地在S₁和S₂中求解最近对问题，分别得到S₁中的最近距离d₁和S₂中的最近距离d₂，令d=min(d₁，d₂。若S的最近对(p,q)之间的距离小于d，则p和q必定分属于集合S₁和S₂。那么，不妨设p∈S₁，q∈S₂，则p和q距离直线x=m的距离均小于d，所以，可以将求解限制在以x=m为中心，宽为2d的垂直带P1和P2中，垂直带之外的任何点对之间的距离都一定大于d。

假设点p(x,y)是集合P₁和P₂中y坐标最小的点，则点p可能在P₁中也可能在P₂中，现在需要找出来和点p之间小于d的点，显然这样的点的y坐标一定位于区间[y,y+d]之间，而且这样的点不会超过8个，因为P1和P2中点的相互之间的距离至少为d，如下图所示。

所以，可以将P₁和P₂中的点按照y坐标升序排列，顺序地处理P₁和P₂中的点p(x,y)，在y坐标区间[y,y+d]内最多取出8个候选点，计算它们和点p之间的距离。

算法

简单起见，假设点集S已按照x坐标升序排列，伪代码描述如下：

输入：按x坐标升序排列的n个点的集合S
输出：最近点对的距离

1. 如果n等于2，则返回(x₁,y₁)和(x₂,y₂)之间的距离，算法结束；
2. 划分：m=S中各点x坐标的中位数；
3. d₁=计算{(x₁,y₁)，···，(x_m,y_m)}的最近对距离；
4. d₂=计算{(x_m,y_m)，···，(x_n,y_n)}的最近对距离；
5. d=min{d₁,d₂}；
6. 依次考察集合S中的点p(x,y)，如果x<=x_m并且x>=x_m-d，则将点p放入集合P₁中；如果x>x_m并且x<=x_m+d，则将点p放入集合P₂中；
7. 将集合P₁和P₂按y坐标升序排列；
8. 对集合P₁和P₂中的每个点p(x,y)，在y坐标区间[y,y+d]内最多取出8个候选点，计算与点p的最近距离d₃；
9. 返回min{d，d₃}；

源代码

#include<iostream>
#include<math.h>
#include<stdlib.h>

using namespace std;

int M[2], N[2], P[2], Q[2];

typedef struct point
{
   
   
    double x,y;
}point;

void Qsortx(point *A, int low, int high)
  {
   
   
    double tempx, tempy;
    int i = low, j = high;
    if(low < high)
      {
   
   
        tempx = A[i].x;
        tempy = A[i].y;
        while(i != j)
         {
   
   
            while(j > i && A[j].x > tempx)
               j--;
            if(i < j)
              {
   
   
                A[i].x = A[j].x;
                A[i].y = A[j].y;
                i++;
              }
            while(i < j && A[i].x < tempx)
              i++;
            if(i<j)
              {
   
   
                A[j].x = A[i].x;
                A[j].y = A[i].y;
                j--;
              }
            A[i].x = tempx;
            A[i].y = tempy