最近點對問題的分治演算法分析與實現

 為了方便起見，在這也先說明一下什麼是＂最近點對問題＂．設有N個點，則易知
從這N個點中任意選取2個點作為一對一共有N(N-1)/2種組合, "最近點對問題"就是從這
N(N-1)/2個點對中求出相距最短的一對, 並求出它們之間的距離.

為了方便演算法的描述, 也把源碼貼一下:
typedef struct _Point {
    double x;
    double y;
} Point;

typedef struct _PointMap {
    Point point;
    int index; /* 這個成員的作用見文章後邊的解釋 */
#define poi_x   point.x
#define poi_y   point.y
} PointMap;
 
跟據問題描述, 我們很容易就會得到一種O(N^2)的解法, 即是: 窮舉所有的點對, 從中
出距離最短的. 該演算法用c語言描述如下:

 

const double EPS=0.00000000001;

const double MAXDIS=9e300;

MinDis( const PointMap M[], int n )
{
    assert( n >= 1 );
    if ( n == 1 )
        return 0.0;
    int i, j;
    double mindis = MAXDIS, tmp;
    for ( i = 0; i < n; i++ )
        for ( j = i + 1; j < n; j++ )
            if ( mindis > ( tmp = Dis( M[i], M[j] ) ) )
                mindis = tmp;
    return sqrt( mindis );
}
夠簡單吧? 呵呵.............

開始進入正題, 我們知道分治的思想就是把一個問題遞迴分解成較小的問題, 然後從子問
題的解中構建出原問題的解. 根據這個思想我們分析"最近點對問題"可知, 我們很容易就
可以對該問題進行分治. 我們可以先取某一個x值作為分界線把所有的點分為兩部分
                          |     
          ... .... .......|.. ..... .... 
          .. ..... .......|...  ...... . 
          ...  ...... . . |..... ......  
          ..... ......  ..|.. .... ... . 
          .. .... ... . ........ . ... . 
          ..... . ... . ..|......... ... 
          ......... ......|............. 
          ...   ... ......| 
                          |                                    
                          x
, 我們可以看出, 最近的點對要麼在分界線的左邊; 要麼在右邊; 也有可能是一個
點在左邊, 另一個點有右邊. 對分界線的左右兩部分遞迴此過程, 此問題也就可以解決了
. 現在的關鍵問題是, 對於第三種情況, 我們應該怎樣有效地求解呢?
                           |     
          ... .... ...|....|.. .|.... .... 
          .. ..... ...|....|... | ...... . 
          ...  ...... |. . |....|. ......  
          ..... ......|  . |.. .|... ... . 
          .. .... ... |. . .....|. . ... . 
          ..... . ... |.. .|....|..... ... 
          ......... ..|....|....|......... 
          ...   ... ..|....|    |
                           |
                 /____/|/____/     
              mdisLRC| mdisLRC       
                           x
                      /_________/
                         雙道帶
假設mdisLRC是遞迴左右兩子部分遞迴返回的兩個最短距離中的較短者, 諾最短點對出現
在第三種情況中, 這兩個只可能出現在中的雙道帶中. 因此我們可以只檢驗雙道帶中
的點. 為了快速地從雙道帶中選出最短的點對(如果它們比midisLRC還要短), 我們可事先
對這些點相對y進行排序, 從而在選取其中一點後可以快速地判斷另一點是否可能成為選
擇, 使得這對點對最短, 即: 如果這兩個點的y軸距離大於mdisLRC, 則它和排在它後邊的
點都不可能成為最短的點(要記得這些點已經按y進行升序排序過了).
該演算法的時間界為: O(NlogN).

還是來看代碼吧:
double
MinDis( const PointMap M[], int n )
{
    assert( n >= 1 );
    if ( n == 1 )
        return 0.0;

    PointMap *P = malloc( n * sizeof(PointMap) );
    if ( P == NULL )
        Error( "Out Of Space!! " );
    PointMap *Q = malloc( n * sizeof(PointMap) );
    if ( Q == NULL )
        Error( "Out Of Space!! " );

    memcpy( P, M, n * sizeof(M[0]) );
    qsort( P, n, sizeof(P[0]), cmp_x ); /* 遞增快速排序 */
    int i;
    /* 對已按x排序的各個項, 分配一個唯一標記index, 因為有可能存在x座標相等的多
     * 個項 */
    for ( i = 0; i < n; i++ )
        P[i].index = i;
    memcpy( Q, P, n * sizeof(P[0]) );
    qsort( Q, n, sizeof(Q[0]), cmp_y );
   
    double mdis = doMinDis( P, Q, 0, n - 1 );

    free( P );
    free( Q );

    return sqrt( mdis );
}

/*
 * ===  FUNCTION  ===========================================================
 *         Name:  doMinDis
 *  Description:  演算法思想: 分治演算法
 *                 已耗用時間: O(NlogN)
 * ==========================================================================
 */
static double
doMinDis( PointMap P[], PointMap Q[],
        int left, int right )
{
    assert( right - left >= 0 );
    if ( right - left == 0 )
        return MAXDIS;
    else if ( right - left == 1 )
        return Dis( Q[left], Q[right] );
    int i, j, k, mid;
    double mdisleft, mdisright, mdisLRC, tmp;
    mid = ( left + right ) / 2; /* 前提: 以x的升序相吻合, 即index是遞增的. */
    for ( i = left, j = left, k = mid + 1; i <= right; i++ )
        if ( Q[i].index <= mid ) /* 以index=mid作為分界,劃分為兩個部分,劃分 */
            P[j++] = Q[i]; /* + 出來的兩個部分, 各自的項的值按y都是遞增的. */
        else
            P[k++] = Q[i];

    mdisleft = doMinDis( Q, P, left, mid );
    mdisright = doMinDis( Q, P, mid + 1, right );
    mdisLRC = Min( mdisleft, mdisright );

    /* P左右部分分別是對y座標有序的, 將其合并到Q. */
    Merge( Q, P, left, mid + 1, right ); /* 目的: 使得整體相對y是有序的 */
    for ( i = left, k = left; i <= right; ++i ) /* 找出雙道帶中的點 */
        if ( fabs( Q[i].poi_x - Q[mid].poi_x ) < mdisLRC )
            P[k++] = Q[i];
    for ( i = left; i < k; i++ ) /* 各項對y是排序的,大大減少了這裡已耗用時間 */
        for ( j = i + 1; j < k && P[j].poi_y - P[i].poi_y < mdisLRC; j++ )
        /* 對poi_y的測試, 如果y座標已大於mdisLRC, 則點對距離肯定大於mdisLRC */
            if ( (tmp = Dis( P[i], P[j] )) < mdisLRC )
                    mdisLRC = tmp;
    return mdisLRC;
}

/*
 * ===  FUNCTION  ======================================================================
 *         Name:  Merge
 *  Description:  把TmpAs合并到P, TmpAs的左右部分各自相對y是升序的
 * =====================================================================================
 */
static void
Merge( PointMap P[], PointMap TmpAs[], int LeftPos,
        int RightPos, int RightEnd )
{
    int LeftEnd, TmpPos;
    LeftEnd = RightPos - 1;
    TmpPos = LeftPos;

    while ( LeftPos <= LeftEnd && RightPos <= RightEnd )
        if ( TmpAs[LeftPos].poi_y < TmpAs[RightPos].poi_y )
            P[TmpPos++] = TmpAs[LeftPos++];
        else
            P[TmpPos++] = TmpAs[RightPos++];

    while ( LeftPos <= LeftEnd )
        P[TmpPos++] = TmpAs[LeftPos++];
    while ( RightPos <= RightEnd )
        P[TmpPos++] = TmpAs[RightPos++];
}

static double
Min( double lh, double rh )
{
    return lh < rh ? lh : rh;
}

static double
Dis( PointMap po1, PointMap po2 ) /* 返回距離的平方和 */
{
    double xdis = po1.poi_x - po2.poi_x;
    double ydis = po1.poi_y - po2.poi_y;
    return xdis * xdis + ydis * ydis;
}

static int
cmp_x( const void *lh, const void *rh )
{
    double tmp = ((PointMap*)lh)->poi_x
        - ((PointMap*)rh)->poi_x;
    if ( tmp > 0 )
        return 1;
    else if ( fabs( tmp ) < EPS )
        return 0;
    else
        return -1;
}

static int
cmp_y( const void *lh, const void *rh )
{
    double tmp = ((PointMap*)lh)->poi_y
        - ((PointMap*)rh)->poi_y;
    if ( tmp > 0 )
        return 1;
    else if ( fabs( tmp ) < EPS )
        return 0;
    else
        return -1;
}