真正的平面上最近点对的n log(n)算法

平面上一堆二维坐标点，找其中两点间的最小距离，这个问题网上一大把n log(n) * log(n)的方法，这些算法的典型例子 http://yzmduncan.iteye.com/blog/1432880 ，原因是在求得中心线两边的最小距离s后，都要对处于中线左右两边s距离以内的点按y坐标排序，这种方法不是n log(n)复杂度（参见维斯在“数据结构与算法分析，c语言版”,中文版p281），要是n log(n)复杂度，要在递归之前就对点按y坐标排序，网上少见这种算法的源码，其中一个似乎有些问题的代码见tge7618291的博客 http://blog.csdn.net/tge7618291/article/details/4753126?reload#comments， 但他在递归中反复利用同样两个事先已经按x, y坐标排好序的p, q数组，这样就不用分配新内存，但是不太好理解，另外，更糟糕的是，他的源代码编译后的程序执行，我发现了一组和最普通的n平方算法标准答案不一样的一组数据点，我留言给作者，但作者似乎已经废弃了他的博客。。。

        我按照 算法导论和维斯在“数据结构与算法分析，c语言版”上介绍的思想，实现了n log(n)复杂度算法，即：递归之前按x, y坐标排好序X, Y数组, 作为参数，其后每次递归调用里把Y数组中每个点按在中心线的左右两边分成YL, YR子数组（每个子数组仍然按y坐标排序）， 把这两个数组分别传人对左右两边的子递归的函数中，而X数组保持不变，这样就比tge7618291的代码好理解，代码如下

// p280_nearestPtPair.cpp : 定义控制台应用程序的入口点。
<pre name="code" class="cpp">// nearPts_ZHENGLI.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"

#include "stdafx.h"
#include <stdlib.h>
//#include <iostream>
//#include <cstdio>
//#include <cstring>
#include <cmath>
#include <stack>
#include <list>
using namespace std;


// 似乎http://blog.csdn.net/lu1012123053/article/details/9825175方法也是n log(n)算法 ，但没细看
// 下面方法部分代码 参考了 http://blog.csdn.net/tge7618291/article/details/4753126?reload#comments
#define min(a,b)    (((a) < (b)) ? (a) : (b))
typedef struct _Point { 
double x;
double y;
//Point(double a, double b):x(a), y(b) {} 
} Point;

typedef struct _PointMap {
Point point;
int index; /* 这个成员的作用见文章后边的解释 */
#define poi_x   point.x
#define poi_y   point.y
} PointMap;

const double EPS=0.00000000001;

const double MAXDIS=9e300;

static doubleDis( 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;
}

//n * n的方法，这个方法算出来的应该是标准答案
double n2MinDis( const PointMap M[], int n )
{
//assert( n >= 1 );
if ( n == 1 )
return 0.0;
int i, j;
pair<Point, Point> pt2;
stack<pair<Point, Point>> stk;
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] ) ) ) {
pt2 = make_pair(M[i].point, M[j].point);
stk.push(pt2);
mindis = tmp;
}
    pt2 = stk.top();
double dist = sqrt( mindis );
printf("the nearest 2  points are: (%f, %f) (%f, %f), dist = %f\n", pt2.first.x, pt2.first.y, pt2.second.x, pt2.second.y, dist);
return dist;
}


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



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;
}

// 以下2个函数是tge7618291的方法，参见 http://blog.csdn.net/tge7618291/article/details/4753126?reload#comments，
//他在递归中反复利用同样两个事先已经按x, y坐标排好序的p, q数组，这样就不用分配新内存，但是不太好理解，
//另外，更糟糕的是，他的源代码编译后的程序执行，我发现了一组和最普通的n平方
//算法标准答案不一样的一组数据点， 见main函数中
static voidMerge( 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 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;
}//-------------------以上2个函数是tge7618291的方法


// -------------------------------------我的方法,每次递归都分配新的存放按Y座标排序的数组
typedef double (*findNearestPtsFunc)( PointMap X[], PointMap Y[],   int l, int r ) ;

double myMinDis( PointMap X[], PointMap Y[], int l, int r) {
if ((r - l) == 1){
double ret = Dis(Y[0],Y[1]);return ret;
}
else if (2 == (r -l)) {
double tmp_dis1 = Dis(Y[0],Y[ 1]);   
double tmp_dis2 = Dis(Y[1],Y[2]);   
double tmp_dis3 = Dis(Y[0],Y[2]);     
double ret = min(min(tmp_dis1,tmp_dis2),tmp_dis3);return ret;
}
int i, j, k, mid;
double mdisleft, mdisright, mdisLRC, tmp;
mid = (l + r)/ 2; /* 前提: 以x的升序相吻合, 即index是递增的.这个mid值为界线，index <= mid的,为左半部分，反之右半部 */
int ttlen = r - l + 1;
PointMap *YL = (PointMap *)malloc( ttlen * sizeof(PointMap) );// 似乎可以以ttlen / 2为大小分配
if ( YL == NULL ){
printf( "yl Out Of Space!! " );
return -1;
}
PointMap *YR = (PointMap *)malloc(ttlen * sizeof(PointMap) );//(n - mid)
if ( YR == NULL ){
printf( "yr Out Of Space!! " );
return -1;
}

// 把按Y值递增排序的第二个参数数组Y[left -->rig]内所有元素，分成在mid左边和右边的两部分，左边
// 的放在数组YL中，中，右边的放在数组YR中 ，而且每部分还是按y值递增的，
for ( i = 0, j = 0, k = 0/*mid + 1*/; i < ttlen; i++ ){
//printf("Y[%d].index = %d\n", i, Y[i].index);  
if ( Y[i].index <= mid )   
YL[j++] = Y[i];   
else    
YR[k++] = Y[i];
}
mdisleft = myMinDis(X, YL, l, mid  );
mdisright = myMinDis( X, YR,  mid + 1, r );
mdisLRC = Min( mdisleft, mdisright );
// YinStrip 表示所有x坐标和中心点的x坐标相差小于mdisLRC(即左右两边递归求的的最小距离）的带状范围内的点，
///按y坐标递增的. */
PointMap *YinStrip = (PointMap *)malloc(ttlen  * sizeof(PointMap) );
if ( YinStrip == NULL ) {
printf( "YinStrip Out Of Space!! " );
return -1;
}
for ( i = 0, j = 0; i < ttlen; i++ ){
if ( abs(X[Y[i].index].poi_x - X[mid].poi_x) < mdisLRC ) 
YinStrip[j++] = Y[i]; 
}
// j是在strip带范围之内的点的个数,每个点 m 需要和其后最多7个点相比较，如果这个点其后点数少于7个，靠 s < j判断
for (int m = 0; m < j; m++ )
for( int s = (m + 1); (s <= (m + 7) && s < j); s++){
tmp = Dis(YinStrip[m], YinStrip[s]);
if ( tmp < mdisLRC)
mdisLRC = tmp;
}
// 注意把本次调用内的内存释放掉
if (YL){
free(YL);
YL = 0;
}

if (YR){
free(YR);
YR = 0;
}

if (YinStrip) {
free(YinStrip);
YinStrip = 0;
}
return mdisLRC;
}

double getMindis(const PointMap M[], int n, findNearestPtsFunc  findNearestPts) 
{
PointMap *X = (PointMap *)malloc( n * sizeof(PointMap) );  
if ( X == NULL )     
printf( "Out Of Space!! " );  
PointMap *Y = (PointMap *)malloc( n * sizeof(PointMap) );   
if ( Y == NULL )    
printf( "Out Of Space!! " ); 
memcpy( X, M, n * sizeof(M[0]) );  
qsort( X, n, sizeof(X[0]), cmp_x ); /* 递增快速排序 */ 
int i;    
/* 对已按x排序的各个项, 分配一个唯一标记index, 因为有可能存在x坐标相等的多个项 */   
for ( i = 0; i < n; i++ )     
X[i].index = i;
for (int j = 0; j < n; j++){
//printf("x[%d].index = %d, pt(x, y)= (%f, %f)\n", j, X[j].index, X[j].poi_x, X[j].poi_y);
} 

memcpy( Y, X, n * sizeof(X[0]) ); 
qsort( Y, n, sizeof(Y[0]), cmp_y );// 初始Y数组是按y值递增顺序排列的，但同时，每个点的index又指明了在X 中的序号 
for (int j = 0; j < n; j++){
//printf("Y[%d].index = %d, poi_y= %f\n", j, Y[j].index, Y[j].poi_y); 
}  
double mdis = findNearestPts( X, Y, 0, n - 1);  
free( X );  
free( Y );
return sqrt( mdis );
}

int _tmain(int argc, _TCHAR* argv[])
{
/*PointMap M[] = {{{1.0, 3.0}, 0}, {{1.5, 3}, 0}, {{5,7}, 0}, {{5.25,7}, 0}, {{3,3}, 0}, {{3,6}, 0},
/*{{1.0, 3.0}, 0}, {{1.5, 3}, 0}, {{5,7}, 0}, {{5.25,7}, 0}, *//*{{3,3}, 0},
{{1.1, 3.0}, 0}, {{1.9, 3}, 0}, {{5.7,7}, 0}, {{5.25,7.3}, 0}, {{39,63}, 0}, {{39.1,62.9}, 0},
{{1.36, 113.0}, 0}, {{11.9, 3}, 0}, {{1.7,17}, 0}, {{9.25,17.3}, 0}, {{9,6.3}, 0}, {{1,6}, 0}
};
*/
// tge7618291的算法算出来有问题的点对
PointMap M[] = {{{1.1, 3.0}, 0}, {{1.9, 3}, 0}, {{5.7,7}, 0}, {{5.25,7.3}, 0}, {{39,63}, 0}, {{31,6}, 0},
{{1.36, 113.0}, 0}, {{11.9, 3}, 0}, {{1.7,17}, 0}, {{9.25,17.3}, 0}, {{9,6.3}, 0}, {{1,6}, 0}
};

/*PointMap M[] = {{{91.0, 3.0}, 0}, {{11.5, 3}, 0}, {{59,37}, 0}, {{75.25,7}, 0}, {{113,3}, 0}, {{0.7,86}, 0},
{{719.0, 0.3}, 0}, {{15, 73}, 0}, {{8.5,7}, 0}, {{525,0.17}, 0}, {{3.11,73}, 0},{{1,8}, 0}, {{1.091,8.05}, 0},
{{1.1, 3.0}, 0}, {{1.9, 3}, 0}, {{5.7,7}, 0}, {{5.25,7.3}, 0}, {{39,63}, 0}, {{31,6}, 0},
{{1.36, 113.0}, 0}, {{11.9, 3}, 0}, {{1.7,17}, 0}, {{9.25,17.3}, 0}, {{9,6.3}, 0}, {{1,6}, 0},

{{5, 73}, 0}, {{85,7}, 0}, {{5,0.17}, 0}, {{9.11,73}, 0},{{1,98}, 0}, {{0.91,18.05}, 0},
{{1.91, 36}, 0}, {{99, 3133}, 0}, {{5, 9970}, 0}, {{25,7.3}, 0}, {{9,63}, 0}, {{31,1}, 0},
{{1.36, 3.0}, 0}, {{119, 3}, 0}, {{17,17}, 0}, {{925,17.3}, 0}, {{90,6.3}, 0}, {{109, 6}, 0}
};*/

int n = sizeof(M) / sizeof(PointMap);
for (int j = 0; j < n; j++){
printf("M[%d].index = %d, pt(x, y)= (%f, %f)\n", j, M[j].index, M[j].poi_x, M[j].poi_y);
}
double disn2= n2MinDis( M,  n );
double dis = getMindis(M,  n , doMinDis );//MinDis(  M,  n );

double my = getMindis(M,  n , myMinDis );
if (disn2 != my)
printf("\n\n---------------------my method 答案%f != 标准答案 %f\n", my, disn2);
else if (disn2 != dis)
printf("\n\n---------------------tge7618291's method 答案 %f != 标准答案 %f\n", dis, disn2);
else
printf("\n\n---------------------3 method get SAME result: %f\n", my);
system("pause");
return 0;
}