最近点对问题的分治算法分析与实现

为了方便起见，在这也先说明一下什么是＂最近点对问题＂．设有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;
}