1. 程式人生 > >《演算法分析與設計》之Closest Pair of Points Problem(最近點對問題)d

《演算法分析與設計》之Closest Pair of Points Problem(最近點對問題)d

       今天的演算法課上主要講了最近點對的問題,在老師的講解下對這個問題有了一個基礎的認識和了解,這裡就先對這個問題做一個簡單的總結吧。

最近點對問題介紹:

       最近點對問題說來其實很簡單,主要就是在二維平面內的n個點中,找出(歐式)距離最近的兩個點來。

問題解決思路:

       解決這個問題最直接的方法就是“蠻力法”(Brute Force)了,其實現過程即為:
  1. 計算出二維平面內所有點對(任意兩點的組合)之間的距離
  2. 逐一比較這些距離的大小,找出其中的最短的距離,則這兩點即為所求的“最近點對”
       當然,蠻力法雖然很容易想到,理解和記憶起來也不困難,在平面內n個點的數目較少時不失為一種解決該問題的簡單方法,但是當問題的規模擴大時,這個方法所需要的大量計算和比較,將耗去相當大的時間和資源。
況且,理想環境下的“小資料”情況在現實生活中幾乎不存在,面對大量的資料時,找到一個更為簡單和快速的方法尤為重要。 這裡,就要引入解決這個問題所需要的一種演算法設計思想——分治策略(Divide and Conquer)

分治策略簡介:

       這裡簡單介紹一下分治策略的思想:

       分治策略求解問題的過程一般可以分為如下三個步驟:

  • 分解問題(Divide):

       這個步驟將一個問題劃分成若干個子問題(注:這些子問題的規模最好是差不多的,這些子問題在形式上和原問題是一樣的,只是規模相對來說更小,方便解決。

       例如在我國,要對全國總計約13億的人口進行一次普查是一個相當複雜且麻煩的問題,但是如果將這個問題的規模由“全國”這一範圍依次劃分成對各省級單位、市級單位、縣區級……進行人口普查,那麼,同樣是人口普查,但是難度上就小了許多了。所以對問題進行分解,往往是解決一個問題的開始。

  • 解決問題(Conquer):

這個步驟需要做的工作是“遞迴地求解出子問題”。這個步驟是要和上一步驟中的分解問題共同進行的。

還是上面的那個例子,當我們接到任務說要對全國人口進行一次普查時,我們首先將這一問題分解成對全國各個省份的人口進行普查,然後嘗試著來解決;可是我們會發現,即使是按省份來普查人口,這個範圍還是太大了,於是我們只好繼續分解問題,按城市、按縣區級來劃分都還是沒有將問題的規模縮小到足以讓我們在短期內能夠普查完這種劃分下的所有人口,那麼說明問題還是沒有被分解到“足夠小”,所以我們還需要將這種“遞迴思想”繼續下去,直至我們將人口普查的範圍縮小到了街道或者村級這樣的範圍,我們差不多就可以進行這一範圍內的人口普查工作了。當然,問題規模也不要太小,例如小到一個家庭的這種程度,這樣統計起來雖然快,但未免有些太浪費資源了……

上面的例子說了這麼多廢話,其實主要想表達的就是“遞迴求解子問題,如果子問題的規模足夠小,則停止遞迴,直接求解。”(出自《演算法導論(原書第3版)》第四章 分治策略 P37)

  • 合併問題(Combine)

雖然這個步驟並沒有在分治策略的名字中有所體現,但是這個步驟卻也是不容忽視的。因為小問題的解決並不意味著大問題的解決,就像上面例子中所提到的,我們在第二個步驟完成後得到的還只是全國各個街道或者村的人口普查結果,而不是我們最終想要得到的全國人口普查結果。

所以要想真正解決原問題,還得把得到的小問題(即子問題)的解依次合併成大問題(即原問題)的解才行。

好了,到這裡為止,分治策略思想的基本情況就介紹完了,下面就開始利用這個思想來解決今天的“最近點對問題”。

根據上面對分治策略思想的介紹,我們要解決最近點對問題的第一步就要從“分解問題”開始。那麼這個問題要怎麼分解呢?

對於平面內的n個點,要怎樣將這個平面劃分成更小的平面呢?如果是按照我們學習數學時劃分象限的方法將平面劃分成上下左右四個小平面的話就會是這樣的:(注:以下圖片均來自老師授課時所用的PPT課件


上面這個劃分看起來好像沒什麼問題,我們的確是按照分治策略把原來的平面分成了規模差不多的四個小平面,但需要注意的是,劃分平面時,我們需要考慮平面內的點的一些特殊分佈情況,像是下面這樣的:


這個時候四個小問題的規模就完全不一樣了,就像是工廠把生產任務分給四個員工,分配任務的方式這麼不平均,必然會招致員工的不滿了。

所以這個時候我們就要考慮,是不是太心急了?一下子就把問題劃分成了四個小問題是不是跨度有些太大了?不如把解決問題的急切心情舒緩一下,先把原平面劃分成像是下面這樣的兩個平面看看吧:

這個時候我們就只要分別求出上圖中左右兩個子平面中的距離最近的兩個點之間的距離就可以了,如下圖所示:


當然,這是在子平面內點數較小的情況下才可以直接求解最近距離的點對的。如果子平面內點數仍然較多,那麼還需要遞迴地對左右兩個子平面繼續劃分下去,直至可以較快求得最近距離點對為止。

如上圖中所計算得到的,左平面中計算得到最近距離為12(單位就暫且忽略吧……),右平面中計算所得的最近距離為21,這樣一比較,12明顯小於21啊,那是不是說明這個平面內的最近點距就是12了呢?

我的回答是:不一定。細心的網友一定會發現,在左右平面的分割線附近,還有兩個相距非常近的點的存在,因為分割線而使得這兩個點分別位於左右兩個平面中了,但是在原平面中,它們則有可能才是相距最近的兩個點,驗證一下吧:


計算出來了,果然,這兩個點之間的距離只有8,比之前求得的最近距離12還要小。

那麼問題來了,要怎麼樣做才能避免此類情況的發生?先別急著看下面的答案,試著自己想一想。解決的方法其實很簡單,我們只要想一下問題是怎麼產生的,就能明白如何解決了。

解決思路:

由於正是兩個平面的分割線反而將平面內最近的兩點給分割開來了,那麼我們只要想辦法將這兩個點再次合併到一個平面內就可以了。可是如果把左右兩個平面再次合併成一個平面,然後計算平面內兩點間最近的距離的話,這不又是變成了“蠻力法”的求解方法了嗎?

是的,如果直接再次合併兩個平面並計算平面內所有點對距離,這對解決這個問題在速度上提升上將會沒有任何效果,反而讓前面的步驟變得多餘。

這時,細心的網友( =( ̄ω ̄)=細心的網友好厲害……)會再次發現,其實只要計算分割線附近的點的距離就可以了。

好的,問題又來了,這個“附近”到底是怎麼定義的呢?多近才能算是附近?

再次陷入困惑的網友請不要忽略了我們前面辛辛苦苦計算出來的左右兩個平面中的點的最近距離!

是的,既然要在分割線附近找比左右兩個平面內計算出的最近距離還要近的距離,那麼範圍自然是不能超過目前已知的最近距離了。

所以,分割線“附近”的定義如下圖所示:


分割線附近的範圍定義好了以後,我們就可以通過計算落在這個範圍內的點之間的距離並找出它們的最近距離了。

最後,再次與之前已經求得的最近距離進行比較就可以得到整個平面內的最近點對了!

至此,今天關於最近點對問題的講解就結束了,後面有時間會把相應的演算法再補充上來。(2015/04/01)

PS:這個日期說了這樣的話,總覺得以後不補充上來也能有很好的藉口吶(~ ̄▽ ̄)~。

(2015/04/15 補充)

時隔半個月(時間過得好快……)來補充一下當初說好的演算法……

特別宣告:這個演算法來自網路(根據搜尋及相關連結,該程式作者為

/**************求最近點對的分治演算法實現,輸入點對,輸出最近點對的距離**************/
#include <iostream>
#include<cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
using namespace std;

const int N = 100005;
const double MAX = 10e100;
const double eps = 0.00001;
typedef struct TYPE
{
    double x, y;
    int index;
} Point;
Point a[N], b[N], c[N];
double closest(Point *, Point *, Point *, int, int);
double dis(Point, Point);
int cmp_x(const void *, const void*);
int cmp_y(const void *, const void*);
int merge(Point *, Point *, int, int, int);
inline double min(double, double);

int main()
{
    int n, i;
    double d;
    scanf("%d", &n);
    while (n)
    {
        for (i = 0; i < n; i++)
            scanf("%lf%lf", &(a[i].x), &(a[i].y));
        qsort(a, n, sizeof(a[0]), cmp_x);
        for (i = 0; i < n; i++)
            a[i].index = i;
        memcpy(b, a, n *sizeof(a[0]));
        qsort(b, n, sizeof(b[0]), cmp_y);
        d = closest(a, b, c, 0, n - 1);
        printf("%.2lf\n", d);
        scanf("%d", &n);
    }
    return 0;
}

double closest(Point a[], Point b[], Point c[], int p, int q)
{
    if (q - p == 1)
        return dis(a[p], a[q]);
    if (q - p == 2)
    {
        double x1 = dis(a[p], a[q]);
        double x2 = dis(a[p + 1], a[q]);
        double x3 = dis(a[p], a[p + 1]);
        if (x1 < x2 && x1 < x3)
            return x1;
        else if (x2 < x3)
            return x2;
        else
            return x3;
    }
    int m = (p + q) / 2;
    int i, j, k;
    double d1, d2;
    for (i = p, j = p, k = m + 1; i <= q; i++)
        if (b[i].index <= m)
            c[j++] = b[i];
    //陣列c左半部儲存劃分後左部的點, 且對y是有序的.
    else
        c[k++] = b[i];
    d1 = closest(a, c, b, p, m);
    d2 = closest(a, c, b, m + 1, q);
    double dm = min(d1, d2);
    merge(b, c, p, m, q); //陣列c左右部分分別是對y座標有序的, 將其合併到b.
    for (i = p, k = p; i <= q; i++)
        if (fabs(b[i].x - b[m].x) < dm)
            c[k++] = b[i];
    //找出離劃分基準左右不超過dm的部分, 且仍然對y座標有序.
    for (i = p; i < k; i++)
    for (j = i + 1; j < k && c[j].y - c[i].y < dm; j++)
    {
        double temp = dis(c[i], c[j]);
        if (temp < dm)
            dm = temp;
    }
    return dm;
}

double dis(Point p, Point q)
{
    double x1 = p.x - q.x, y1 = p.y - q.y;
    return sqrt(x1 *x1 + y1 * y1);
}

int merge(Point p[], Point q[], int s, int m, int t)
{
    int i, j, k;
    for (i = s, j = m + 1, k = s; i <= m && j <= t;)
    {
        if (q[i].y > q[j].y)
            p[k++] = q[j], j++;
        else
            p[k++] = q[i], i++;
    }
    while (i <= m)
        p[k++] = q[i++];
    while (j <= t)
        p[k++] = q[j++];
    memcpy(q + s, p + s, (t - s + 1) *sizeof(p[0]));
    return 0;
}

int cmp_x(const void *p, const void *q)
{
    double temp = ((Point*)p)->x - ((Point*)q)->x;
    if (temp > 0)
        return 1;
    else if (fabs(temp) < eps)
        return 0;
    else
        return - 1;
}

int cmp_y(const void *p, const void *q)
{
    double temp = ((Point*)p)->y - ((Point*)q)->y;
    if (temp > 0)
        return 1;
    else if (fabs(temp) < eps)
        return 0;
    else
        return - 1;
}

inline double min(double p, double q)
{
    return (p > q) ? (q): (p);
}

非常感謝中南大學鄭瑾老師的授課,把這個問題講解的非常通透,以上講解都是按照老師上課時候的思路加上個人的理解整理出來的,如果我對這個問題的介紹存在任何疏漏,歡迎批評指正!
本文完全原創,所引用到的文字、程式碼及圖片均已標明出處,轉載請註明出處,謝謝合作!