求平面最近点对的核心思想乃是二分,用递归实现。具体操作如下:
若点的个数很少(比如小于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->x - arg2->x;
double y_diff = arg1->y - 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->x - 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->y - 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);
}