背景：

計算機應用中經常採用點、圓等簡單的幾何物件表示物理實體，常需要了解其鄰域中其他幾何物件的信息

• 例如：在空中交通控制中，若將飛機作為空間中的一個點來處理，則具有最大碰撞危險的兩架飛機所處的點對，就是這個空間中距離最接近的一對點。

這類問題是計算幾何學中研究的基本問題之一。

我們著重考慮平面上的最接近點對問題。

最接近點對問題：

給定平面上的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];
}