A Za, A Za, Fighting...

坚信:勤能补拙

2011分治-平面最近点对(附C++源代码)

求平面最近点对的核心思想乃是二分,用递归实现。具体操作如下:

     若点的个数很少(比如小于3或者小于5),就直接枚举。

     如点的个数很多,按现将所有点按X排序,并按X坐标平均的分成左右两个部分(假设分割线为X=nx),分别求出两边的最短距离minl与minr并令ans=min(minl,minr)。

     求出左右两边的最小值之后,剩下的工作就是合并。易见若该点集存在点对(a,b)的最近距离小于ans,则a,b一定分别在x=nx的两边,切nx-a.x与nx-b.x的绝对值肯定小于ans。

     据此我们可以将点集中所有X值在(nx-ans,nx+ans)的点都选出来,那么满足条件的(a,b)肯定都在其中。

     易见若存在(a,b)两点他们之间的距离小于ans,那么a.y-b.y的绝对值也肯定小于ans。

     综上存在(a,b)两点他们之间的距离小于ans那,(a,b)一定在一个长为2*ans宽为ans的矩形之中。而 且这个矩形被X=nx平分成两个ans*ans的矩形,由于无论是在左边还是在右边,任意两点的之间的距离总是小于等于ans的,所以两个ans*ans 的矩形中最多只有4个点(分别在四个顶点上),长为2*ans宽为ans的矩形最多有8个点。

     据此我们将所有X值在(nx-ans,nx+ans)的点按他们的Y值进行排序。依次看每个点与它之后的7个点的距离是否小于ans,若小于则更新ans,最后求出来的结果就是平面最近点对的距离。保留产生该距离的两个点即可得到最近点对。

     练手题目:Pku2107,Vijos1012

附C++代码(Pku2107):

#include <iostream>
#include <cmath>

const long maxsize = 100000;

typedef struct 

double x, y; 
} PointType;

long list[maxsize], listlen,n;
PointType point[maxsize];

int sortcmp(const void *,const void *); 
double dis(PointType,PointType);
double getmin(double,double);
int listcmp(const void *,const void *); 
double shortest(long,long);
int init(void);

int main() 

while (init())
   printf("%.2lf\n",shortest(0, n - 1)/2);    
return 0;
}

int sortcmp(const void *a, const void *b) 

if (((PointType*)a)->x < ((PointType*)b)->x)    
   return -1;   
else 
   return 1; 
}

double dis(PointType a, PointType b) 

return sqrt((a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y)); 
}

double getmin(double a, double b) 

return a<b?a:b;
}

int listcmp(const void *a, const void *b) 

if (point[*(int*)a].y < point[*(int*)b].y)    
   return -1;   
else 
   return 1; 
}

double shortest(long l, long r) 

if (r - l == 1) 
   return dis(point[l], point[r]);   
if (r - l == 2)    
   return getmin(getmin(dis(point[l], point[l+1]), dis(point[l], point[r])), dis(point[l+1], point[r]));   
long i, j, mid = (l + r) >> 1;   
double curmin = getmin(shortest(l, mid), shortest(mid + 1, r));   
listlen = 0;   
for (i = mid; i >= l && point[mid+1].x - point[i].x <= curmin; i --)    
   list[listlen++] = i;   
for (i = mid + 1; i <= r && point[i].x - point[mid].x <= curmin; i ++)    
   list[listlen++] = i;   
qsort(list, listlen, sizeof(list[0]), listcmp);   
for (i = 0; i < listlen; i ++)    
   for (j = i + 1; j < listlen && point[list[j]].y - point[list[i]].y <= curmin; j ++)     
    curmin = getmin(curmin, dis(point[list[i]], point[list[j]]));   
return curmin; 
}

int init(void)
{
int i;
scanf("%d", &n);      
for (i=0;i<n;i++) 
   scanf("%lf%lf",&point[i].x,&point[i].y);      
qsort(point,n,sizeof(point[0]),sortcmp);
return n;
}


自己写的代码:

/*
 * Problem(classic):
 *    there're many points in a plane surface, find the nearest two points
 *    see: <CLRS> 33.4 section
 
*/
#include
<stdio.h>
#include
<stdlib.h>
#include
<string.h>
#include
<math.h>
#define INF 0x7FFFFFFF
#define NUM_MAX 100000
#define THRESHOLD 3

struct Point {
    
double x, y;
};
struct Point points[NUM_MAX];
int total, yindex[NUM_MAX];

double
min(
double arg1, double arg2)
{
    
return (arg1 <= arg2 ? arg1 : arg2);
}

double
distance(
struct Point *arg1, struct Point *arg2)
{
    
double x_diff = arg1->- arg2->x;
    
double y_diff = arg1->- arg2->y;
    
return sqrt(x_diff*x_diff + y_diff*y_diff);
}

int
compare_x(
const void *arg1, const void *arg2)
{
    
struct Point *p1 = (struct Point *)arg1;
    
struct Point *p2 = (struct Point *)arg2;
    
return (p1->- p2->x);
}

int
compare_y(
const void *arg1, const void *arg2)
{
    
struct Point *p1 = points + *(int *)arg1;
    
struct Point *p2 = points + *(int *)arg2;
    
return (p1->- p2->y);
}

void
init_preprocess()
{
    
int i;
    scanf(
"%d"&total);
    
for(i=0; i<total; ++i)
        scanf(
"%lf %lf"&points[i].x, &points[i].y);
    qsort(points, total, 
sizeof(struct Point), compare_x);
}

double
find_nearest(
int begin, int end)
{
    
int i, j;
    
double ret = INF;
    
if(end-begin+1 <= THRESHOLD) {
        
for(i=begin; i<end; ++i) {
            
for(j=i+1; j<=end; ++j)
                ret 
= min(ret, distance(points+i, points+j));
        }
        
return ret;
    }
    
int mid = begin + ((end - begin) >> 1);
    
double dis = min(find_nearest(begin, mid), find_nearest(mid+1, end));
    
int len = 0;
    
for(j=mid; j>=begin && points[mid+1].x-points[j].x<=dis; --j)
        yindex[len
++= j;
    
for(j=mid+1; j<=end && points[j].x-points[mid].x<=dis; ++j)
        yindex[len
++= j;
    qsort(yindex, len, 
sizeof(int), compare_y);
    ret 
= dis;
    
for(i=0; i<=len-7++i) {
        
for(j=i+1; j<=i+7++j)
            ret 
= min(ret, distance(points+yindex[i], points+yindex[j]));
    }
    
return ret;
}

double
brute_force(
int begin, int end)
{
    
double ret = INF;
    
int i, j;
    
for(i=begin; i<end; ++i) {
        
for(j=i+1; j<=end; ++j)
            ret 
= min(ret, distance(points+i, points+j));
    }
    
return ret;
}

int
main(
int argc, char **argv)
{
    init_preprocess();
    
double ret = find_nearest(0, total-1);
    printf(
"\nNearest Distance[Brute Force]: %f\n", brute_force(0, total-1));
    printf(
"\nNearest Distance: %f\n", ret);
}


posted on 2011-09-03 23:17 simplyzhao 阅读(1104) 评论(1)  编辑 收藏 引用 所属分类: R_找工复习2011

评论

# re: 2011分治-平面最近点对(附C++源代码) 2011-11-05 11:45 co

好东西啊,非常感谢!写的很简洁  回复  更多评论   


只有注册用户登录后才能发表评论。
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理


导航

<2011年11月>
303112345
6789101112
13141516171819
20212223242526
27282930123
45678910

统计

常用链接

留言簿(1)

随笔分类

随笔档案

搜索

最新评论

阅读排行榜

评论排行榜