背景：

计算机应用中经常采用点、圆等简单的几何对象表示物理实体，常需要了解其邻域中其他几何对象的信息

• 例如：在空中交通控制中，若将飞机作为空间中的一个点来处理，则具有最大碰撞危险的两架飞机所处的点对，就是这个空间中距离最接近的一对点。

这类问题是计算几何学中研究的基本问题之一。

我们着重考虑平面上的最接近点对问题。

最接近点对问题：

给定平面上的n个点，找出其中的一对点，使得在n个点组成的所有点对中，该点对的距离最小。

直观解法：

• 将每一个点与其他n-1个点的距离算出，找出最小距离。
• 时间复杂度：T(n)=n(n-1)/2+n=O(n2)，ps：求出n(n-1)/2个点对距离的时间+求出最短距离的时间

分治法：

• 分解：将n个点的集合分成大小近似相等的两个子集。
• 求解：递归地求解两个子集内部的最接近点对。
• 合并（关键问题）：从子空间内部最接近点对，和两个子空间之间的最接近点对中，选择最接近点对。

1、一维最接近点对问题

算法思路：

考虑将所给的n个点的集合S分成2个子集S1和S2，每个子集中约有n/2个点，然后在每个子集中递归地求其最接近的点对。

关键的问题是分治法中的合并步骤：

• 由S1和S2的最接近点对，还要求得原集合S中的最接近点对，因为S1和S2的最接近点对未必就是S的最接近点对。
• 如果这2个点分别在S1和S2中，则对于S1中任一点p，S2中最多只有n/2个点与它构成最接近点对的候选者。
• 仍需做n^2/4(n/2*n/2)次计算和比较才能确定S的最接近点对。
• 合并步骤耗时为O(n^2)，整个算法所需计算时间T(n)应满足:　T(n)=2T(n/2)+O(n^2)。
• 因此，问题出在合并步骤耗时太多。

用x轴上某个点m将S划分为2个子集S1和S2，使得S1={x∈S|x≤m}；S2={x∈S|x>m}。

因此，所有p∈S1和q∈S2有p<q，递归地在S1和S2上找出其最接近点对{p1,p2}和{q1,q2}，并设d=min{|p1-p2|,|q1-q2|}。

S中的最接近点对或者是{p1,p2}，或者是{q1,q2}，或者是某个{p3,q3}，其中p3∈S1且q3∈S2。

如果S的最接近点对是{p3,q3}，即|p3-q3|<d，则p3和q3两者与m的距离不超过d，即|p3-m|<d，|q3-m|<d。

可得，p3∈(m-d,m]，q3∈(m,m+d]。

由于在S1中，每个长度为d的半闭区间至多包含一个点（否则必有两点距离小于d）。

并且m是S1和S2的分割点，因此(m-d,m]中至多包含S中的一个点，且为S1中最大点。

同理，(m,m+d]中也至多包含S中的一个点，且为S2中最小点。

因此，我们用线性时间就能找到区间(m-d,m]和(m,m+d]中所有点，即p3和q3。

按这种分治策略，合并步可在O(n)时间内完成。

还需考虑的一个问题：分割点m的选取，及S1和S2的划分。

• 选取分割点m的一个基本要求是由此导出集合S的一个线性分割。
• 即S=S1∪S2 ，S1∩S2=Φ，且S1={x|x≤m}；S2={x|x>m}。
• 容易看出，如果选取m=[max(S)+min(S)]/2，可以满足线性分割的要求。
• 选取分割点后，再用O(n)时间即可将S划分成S1={x∈S|x≤m}和S2={x∈S|x>m}。

划分可能出现的问题：造成划分出的子集S1和S2的不平衡

• 例如在最坏情况下，|S1|=1，|S2|=n-1。
• 由此，产生的最坏情况下所需的计算时间T(n)=T(n-1)+O(n)。
• 该递归方程的解为：T(n)=O(n^2)
• 可用分治法中“平衡子问题”的方法来解决，如用S的n个点的坐标的中位数来作分割点。
• 用选取中位数的线性时间算法，可在O(n)时间内确定一个平衡的分割点m。

确定平衡点采用m=[max(S)+min(S)]/2的方法，程序如下：

#include <ctime>
#include <iostream>

using namespace std;

//点对结构体
struct Pair {
    float d;//点对距离
    float d1, d2;//点对坐标
};

float Max(float s[], int p, int q);

float Min(float s[], int p, int q);

template<class Type>
void Swap(Type &x, Type &y);

template<class Type>
int Partition(Type s[], Type x, int l, int r);

Pair Cpair(float s[], int l, int r);

int main() {
    srand((unsigned) time(NULL));
    float s[] = {20.14, 3.04, 8.85, 31.72, 40.97, 81.27, 41.15, 45.13, 25.5, 81.84, 3.96, 5.18, 30.82, 72.23, 21.13,
                 57.59, 76.39, 60.28, 87.88, 13.67, 1.22, 7.82, 9.23, 29.09, 30.15};
    Pair d;
    d = Cpair(s, 0, 24);
    cout << endl << "最近点对坐标为：(d1:" << d.d1 << ",d2:" << d.d2 << ")";
    cout << endl << "这两点距离为：" << d.d << endl;
    return 0;
}


float Max(float s[], int l, int r)//返回s[]中的最大值
{
    float s_max = s[l];
    for (int i = l + 1; i <= r; i++)
        if (s_max < s[i])
            s_max = s[i];
    return s_max;
}

float Min(float s[], int l, int r)//返回s[]中的最小值
{
    float s_min = s[l];
    for (int i = l + 1; i <= r; i++)
        if (s_min > s[i])
            s_min = s[i];
    return s_min;
}

template<class Type>
void Swap(Type &x, Type &y) {
    Type temp = x;
    x = y;
    y = temp;
}

template<class Type>
int Partition(Type s[], Type x, int l, int r) {
    int i = l - 1, j = r + 1;

    while (true) {
        while (s[++i] < x && i < r);
        while (s[--j] > x);
        if (i >= j) {
            break;
        }
        Swap(s[i], s[j]);
    }
    return j;
}

//返回s[]中的具有最近距离的点对及其距离
Pair Cpair(float s[], int l, int r) {
    Pair min_d = {99999, 0, 0};//最短距离
    if (r - l < 1) return min_d;
    float m1 = Max(s, l, r), m2 = Min(s, l, r);
    float m = (m1 + m2) / 2;//找出点集中的中位数

    //将点集中的各元素按与m的大小关系分组
    int j = Partition(s, m, l, r);

    Pair d1 = Cpair(s, l, j), d2 = Cpair(s, j + 1, r);//递归
    float p = Max(s, l, j), q = Min(s, j + 1, r);

    //返回s[]中的具有最近距离的点对及其距离
    if (d1.d < d2.d) {
        if ((q - p) < d1.d) {
            min_d.d = (q - p);
            min_d.d1 = q;
            min_d.d2 = p;
            return min_d;
        } else return d1;
    } else {
        if ((q - p) < d2.d) {
            min_d.d = (q - p);
            min_d.d1 = q;
            min_d.d2 = p;
            return min_d;
        } else return d2;
    }
}

该算法的分割步骤和合并步骤总共耗时O(n)。因此，算法耗费的计算时间T(n)满足递归方程：                         

解得：T(n)=O(nlogn)

2、二维最接近点对问题

①选取二维平面的一条垂直平分线L：x=m作为分割线。

━其中m为S中各点x坐标的中位数，由此将S分割为S1={p∈S|px≤m}和S2={p∈S|px>m}。

②递归地在S1和S2上找出其最小距离d1和d2。

━现设d=min(d1,d2)。若S的最接近点对(p,q)之间的距离d(p,q)<d，则p必∈S1和q必∈S2。

③如果用符号P1和P2分别表示直线 L 的左右两边宽为d的区域，则p∈P1，q∈P2。

问题分析：

• 在一维的情形，距分割点距离为d的2个区间(m-d,m](m,m+d]中最多各有S中一个点。
• 因而这2个点成为唯一的末检查过的最接近点对候选者。
• 在二维的情形，在最坏情况下有n2/4对这样的候选者。
• 但是，P1和P2中的点具有稀疏性质，它使我们不必检查所有这n^2/4对候选者。

问题的优化方案：

①考虑P1中的任意一点，它若与P2中的点q构成最接近点对的候选者，则必有：distance(p，q)＜d。

②P2中满足条件的点一定落在矩形R中，矩形R的大小为：d×2d。

③由d的定义可知：P2中任何2个点（qi∈S）的距离都不小于d，由此可以推出矩形R中最多只有6个S中的点。

④因此，在分治法的合并步骤中最多只需要检查6×n/2=3n个候选者。

数学证明：

（一）合并时最多只需检查6×n/2=3n个候选者（矩形R中最多只有6个S中的点）

补充知识：

抽屉原理：把m个元素任意放入n（n<m）个集合，则一定有一个集合呈至少要有k个元素。

其中 k= [m/n]([]表示向上取整）。

• 将矩形R长为2d的边3等分，将长为d的边2等分。
• 由此导出6个(d/2)×(2d/3)的小矩形，如图。
• 若矩形R中有多于6个S中的点，由抽屉原理易知，至少有一个小矩形中有2个以上S中的点
• 设u，v是位于同一小矩形中的2个点，则，
• 即：distance(u,v)<d，这与d的意义相矛盾（P2内不可能再出现距离<d的两个点），命题得证。

⑤如何确定需要检查的6个点？

• 可以将p和P2中所有S2的点投影到垂直线L上。
• 由于能与p点一起构成最接近点对候选者的S2中的点一定在矩形R中。
• 所以它们在直线 L 上的投影点与 p 在 L 上投影点的距离小于d。
• 根据上述分析，这种投影点最多只有6个。因此，若将区域P1和P2中所有S中的点按其y坐标排好序。
• 则：对P1中的所有点，只需一次扫描就可以找出所有候选者。
• 对排好序的点作一次扫描，可以找出所有最接近点对的候选者。
• 对P1中每个点，最多只需检查P2中排好序的相继6个点。

#include<time.h>
#include<iostream>
#include<cmath>

using namespace std;
const int M = 50;

//用类PointX和PointY表示依x坐标和y坐标排好序的点  
class PointX {
public:
    int operator<=(PointX a) const { return (x <= a.x); }

    int ID; //点编号
    float x, y; //点坐标
};

class PointY {
public:
    int operator<=(PointY a) const { return (y <= a.y); }

    int p; //同一点在数组x中的坐标
    float x, y; //点坐标
};

float Random();

template<class Type>
float dis(const Type &u, const Type &v);

bool Cpair2(PointX X[], int n, PointX &a, PointX &b, float &d);

void closest(PointX X[], PointY Y[], PointY Z[], int l, int r, PointX &a, PointX &b, float &d);

template<typename Type>
void Copy(Type a[], Type b[], int left, int right);

template<class Type>
void Merge(Type c[], Type d[], int l, int m, int r);

template<class Type>
void MergeSort(Type a[], Type b[], int left, int right);

int main() {
    srand((unsigned) time(NULL));
    int length;

    cout << "请输入点对数：";
    cin >> length;

    PointX X[M];
    cout << "随机生成的二维点对为：" << endl;

    for (int i = 0; i < length; i++) {
        X[i].ID = i;
        X[i].x = Random();
        X[i].y = Random();
        cout << "(" << X[i].x << "," << X[i].y << ") ";
    }

    PointX a;
    PointX b;
    float d;

    Cpair2(X, length, a, b, d);

    cout << endl;
    cout << "最邻近点对为：(" << a.x << "," << a.y << ")和(" << b.x << "," << b.y << ") " << endl;
    cout << "最邻近距离为： " << d << endl;

    return 0;
}

float Random() {
    float result = rand() % 10000;
    return result * 0.01;
}

//平面上任意两点u和v之间的距离可计算如下  
template<class Type>
inline float dis(const Type &u, const Type &v) {
    float dx = u.x - v.x;
    float dy = u.y - v.y;
    return sqrt(dx * dx + dy * dy);
}

bool Cpair2(PointX X[], int n, PointX &a, PointX &b, float &d) {
    if (n < 2) return false;

    PointX *tmpX = new PointX[n];
    MergeSort(X, tmpX, 0, n - 1);

    PointY *Y = new PointY[n];
    for (int i = 0; i < n; i++) //将数组X中的点复制到数组Y中
    {
        Y[i].p = i;
        Y[i].x = X[i].x;
        Y[i].y = X[i].y;
    }

    PointY *tmpY = new PointY[n];
    MergeSort(Y, tmpY, 0, n - 1);

    PointY *Z = new PointY[n];
    closest(X, Y, Z, 0, n - 1, a, b, d);

    delete[]Y;
    delete[]Z;
    delete[]tmpX;
    delete[]tmpY;
    return true;
}

void closest(PointX X[], PointY Y[], PointY Z[], int l, int r, PointX &a, PointX &b, float &d) {
    if (r - l == 1) //两点的情形
    {
        a = X[l];
        b = X[r];
        d = dis(X[l], X[r]);
        return;
    }

    if (r - l == 2) //3点的情形
    {
        float d1 = dis(X[l], X[l + 1]);
        float d2 = dis(X[l + 1], X[r]);
        float d3 = dis(X[l], X[r]);

        if (d1 <= d2 && d1 <= d3) {
            a = X[l];
            b = X[l + 1];
            d = d1;
            return;
        }

        if (d2 <= d3) {
            a = X[l + 1];
            b = X[r];
            d = d2;
        } else {
            a = X[l];
            b = X[r];
            d = d3;
        }
        return;
    }

    //多于3点的情形，用分治法   
    int m = (l + r) / 2;
    int f = l, g = m + 1;

    //在算法预处理阶段，将数组X中的点依x坐标排序，将数组Y中的点依y坐标排序  
    //算法分割阶段，将子数组X[l:r]均匀划分成两个不想交的子集，取m=(l+r)/2  
    //X[l:m]和X[m+1:r]就是满足要求的分割。  
    for (int i = l; i <= r; i++) {
        if (Y[i].p > m) Z[g++] = Y[i];
        else Z[f++] = Y[i];
    }

    closest(X, Z, Y, l, m, a, b, d);
    float dr;

    PointX ar, br;
    closest(X, Z, Y, m + 1, r, ar, br, dr);

    if (dr < d) {
        a = ar;
        b = br;
        d = dr;
    }

    Merge(Z, Y, l, m, r);//重构数组Y

    //d矩形条内的点置于Z中  
    int k = l;
    for (int i = l; i <= r; i++) {
        if (fabs(X[m].x - Y[i].x) < d) {
            Z[k++] = Y[i];
        }
    }

    //搜索Z[l:k-1]  
    for (int i = l; i < k; i++) {
        for (int j = i + 1; j < k && Z[j].y - Z[i].y < d; j++) {
            float dp = dis(Z[i], Z[j]);
            if (dp < d) {
                d = dp;
                a = X[Z[i].p];
                b = X[Z[j].p];
            }
        }
    }
}

template<class Type>
void Merge(Type c[], Type d[], int l, int m, int r) {
    int i = l, j = m + 1, k = l;
    while ((i <= m) && (j <= r)) {
        if (c[i] <= c[j]) {
            d[k++] = c[i++];
        } else {
            d[k++] = c[j++];
        }
    }

    if (i > m) {
        for (int q = j; q <= r; q++) {
            d[k++] = c[q];
        }
    } else {
        for (int q = i; q <= m; q++) {
            d[k++] = c[q];
        }
    }
}

template<class Type>
void MergeSort(Type a[], Type b[], int left, int right) {
    if (left < right) {
        int i = (left + right) / 2;
        MergeSort(a, b, left, i);
        MergeSort(a, b, i + 1, right);
        Merge(a, b, left, i, right);//合并到数组b
        Copy(a, b, left, right);//复制回数组a
    }
}

template<typename Type>
void Copy(Type a[], Type b[], int left, int right) {
    for (int i = left; i <= right; i++)
        a[i] = b[i];
}