平面最近点对
纠错
04 Oct 2015
描述
给了平面上若干个点的坐标，求最近的两个点之间的距离。
可供练习此问题的地方有:
思路
方法一：Brute Force
$O(n^2)$，不多说了
方法一的优化
事先对所有点按x坐标值排好序，然后通过x坐标差值缩小搜索范围。这个方法的实际效果是非常好的，但是最坏时间复杂度依然是$O(n^2)$，比如，当所有点的横坐标相等的时候。
该方法过了UVa 10245，但是ZOJ和HDU都TLE。下面是在UVa提交并AC的参考代码(C++):
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
const int MAXN = 100010;
struct Point
{
double x;
double y;
}P[MAXN];
int N;
bool cmpX(Point p1, Point p2)
{
if (p1.x != p2.x) {
return p1.x < p2.x;
}
return p1.y < p2.y;
}
double dist(Point p1, Point p2)
{
return sqrt( (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y) );
}
double getClosestDis(Point Px[], int n)
{
if (n == 1) {
return INFINITY;
}
if (n == 2) {
return dist(Px[0], Px[1]);
}
double d = INFINITY;
for (int i = 0; i < n; ++i) {
for (int j = i + 1; j < n; ++j) {
double temp = dist(Px[i], Px[j]);
if (Px[j].x - Px[i].x >= d) {
break;
}
if (temp < d) {
d = temp;
}
}
}
return d;
}
int main()
{
double rst;
while (scanf("%d", &N) && N) {
for (int i = 0; i < N; ++i) {
scanf("%lf%lf", &P[i].x, &P[i].y);
}
sort(P, P + N, cmpX);
rst = getClosestDis(P, N);
if (rst >= 10000.00) {
printf("INFINITY\n");
} else {
printf("%.4lf\n", rst);
}
//printf("%.2lf\n", rst / 2);
}
return 0;
}
方法二：Divide & Conquer
- 把所有点按x坐标值从小到大排序
- 按x坐标把点集划分为左、右两个子集
- 分别计算左、右两个子集的最近点对距离
- 合并，取左右子集最近点对距离的较小者（设为$d$）；同时还要考虑一个点在左子集，另一个点在右子集的情况——只需考虑出现在距离分割线$d$以内范围中的点，在二维平面上是一个长条形区域。筛选出该长条形区域内的所有点，按y坐标值对它们排序（排序的目的是剪枝），然后找这个点集中的最近点对，如果比$d$小，就更新$d$
以上2~4步递归地进行。
划分比较容易，关键是合并的复杂度。筛选出位于长条形区域内的点，须要扫一遍整个点集，$O(n)$；按y坐标值对它们排序，最快$O(n \log{n})$；然后找出最小距离，扫一遍，$O(n)$。所以，合并的复杂度是$O(n \log{n})$。整个算法的复杂度：$O(\log{n} * (n \log{n})) = O(n (\log{n})^2)$
分治的递归方程： \(T(n) = 2 T\left( \frac{n}{2} \right) + O(n \log{n})\)
该方法已经可以通过所有的OJ了。以HDU 1007为例，参考的AC代码(C++)如下:
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
const int MAXN = 100010;
struct Point
{
double x;
double y;
}P[MAXN], strip[MAXN];;
int N;
bool cmpX(Point p1, Point p2)
{
if (p1.x != p2.x) {
return p1.x < p2.x;
}
return p1.y < p2.y;
}
bool cmpY(Point p1, Point p2)
{
if (p1.y != p2.y) {
return p1.y < p2.y;
}
return p1.x < p2.x;
}
double dist(Point p1, Point p2)
{
return sqrt( (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y) );
}
double getStripDis(Point strip[], int n, double d)
{
double rst = d;
sort(strip, strip + n, cmpY);
for (int i = 0; i < n; ++i) {
for (int j = i + 1; j < n && (strip[j].y - strip[i].y) < rst; ++j) {
if (dist(strip[i], strip[j]) < rst) {
rst = dist(strip[i], strip[j]);
}
}
}
return rst;
}
double getClosestDis(Point Px[], int n)
{
if (n == 1) {
return INFINITY;
}
if (n == 2) {
return dist(Px[0], Px[1]);
}
int mid = n / 2;
double x_mid = Px[mid].x;
double d_left = getClosestDis(Px, mid);
double d_right = getClosestDis(Px + mid, n - mid);
double d = min(d_left, d_right);
int k = 0;
for (int i = 0; i < n; ++i) {
if (abs(Px[i].x - x_mid) < d) {
strip[k++] = Px[i];
}
}
double d_strip = getStripDis(strip, k, d);
double d_final = d < d_strip? d:d_strip;
return d_final;
}
int main()
{
double rst;
while (scanf("%d", &N) && N) {
for (int i = 0; i < N; ++i) {
scanf("%lf%lf", &P[i].x, &P[i].y);
}
sort(P, P + N, cmpX);
rst = getClosestDis(P, N);
printf("%.2lf\n", rst / 2);
}
return 0;
}
分治法的优化
注意到在合并操作中，每一次递归都要依y坐标值对strip数组内的点排序，使得合并的复杂度最坏能到$O(n\log{n})$。如果没有排序这一步，那么扫描一遍即可，复杂度即是构建strip数组的复杂度—$O(n)$。于是，就能将求解整个问题的复杂度降至$O(n\log{n})$。
这是可以做到的——通过预处理。为此，需要另开一个数组(设为Py)，在一开始就将所有的点按’Y’坐标值从小到大排好序。然后在递归过程中，按x坐标值做divide操作的同时，也对Py数组进行划分。
以HDU 1007为例，优化后的AC代码(C++)如下：
#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
using namespace std;
const int MAXN = 100010;
struct Point
{
double x;
double y;
}P[MAXN], P2[MAXN], strip[MAXN];
int N;
bool cmpX(Point p1, Point p2)
{
if (p1.x != p2.x) {
return p1.x < p2.x;
}
return p1.y < p2.y;
}
bool cmpY(Point p1, Point p2)
{
if (p1.y != p2.y) {
return p1.y < p2.y;
}
return p1.x < p2.x;
}
double dist(Point p1, Point p2)
{
return sqrt( (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y) );
}
double getStripDis(Point strip[], int n, double d)
{
double rst = d;
for (int i = 0; i < n; ++i) {
for (int j = i + 1; j < n && (strip[j].y - strip[i].y) < rst; ++j) {
if (dist(strip[i], strip[j]) < rst) {
rst = dist(strip[i], strip[j]);
}
}
}
return rst;
}
double getClosestDis(Point Px[], Point Py[], int n)
{
if (n == 1) {
return INFINITY;
}
if (n == 2) {
return dist(Px[0], Px[1]);
}
int mid = n / 2;
int x_mid = Px[mid].x;
Point Py_left[mid+5];
Point Py_right[n-mid+5];
int nl = 0, nr = 0;
// 对Py数组进行划分
for (int i = 0; i < n; ++i) {
if (Py[i].x <= Px[mid-1].x && nl < mid) {
Py_left[nl++] = Py[i];
} else {
Py_right[nr++] = Py[i];
}
}
// conquer
double d_left = getClosestDis(Px, Py_left, mid);
double d_right = getClosestDis(Px + mid, Py_right, n - mid);
double d = min(d_left, d_right);
// 构建strip数组
int k = 0;
for (int i = 0; i < n; ++i) {
if (abs(Py[i].x - x_mid) < d) {
strip[k++] = Py[i];
}
}
// 求strip内的最近点对
double d_strip = getStripDis(strip, k, d);
// 合并子问题的解
double d_final = d < d_strip? d:d_strip;
return d_final;
}
int main()
{
double rst;
while (scanf("%d", &N) && N) {
for (int i = 0; i < N; ++i) {
scanf("%lf%lf", &P[i].x, &P[i].y);
P2[i].x = P[i].x;
P2[i].y = P[i].y;
}
sort(P, P + N, cmpX);
sort(P2, P2 + N, cmpY);
rst = getClosestDis(P, P2, N);
printf("%.2lf\n", rst / 2);
}
return 0;
}
参考
[1] CMSC 451: Closest Pair of Points [2] Divide and Conquer | Set 2 (Closest Pair of Points) [3] Closest Pair of Points | O(nlogn) Implementation [4] 国立台湾师范大学: Closest Pair
-----EOF-----