求平面最近點對的核心思想乃是二分，用遞歸實現。具體操作如下：
若點的個數很少（比如小于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);
}