计算几何——距离平面上最近的点对
给定平面上 n n n 个点,找出其中的一对点的距离,使得在这 n n n 个点的所有点对中,该距离为所有点对中最小的。
分治法
算法每次输入的参数为 P ⊆ Q P \subseteq Q P⊆Q,数组 X X X和 Y Y Y,其中 P P P为待解决的点集合, Q Q Q为全体点集, X X X是 P P P按照 x x x坐标排序的结果, Y Y Y是按照 y y y排序的结果.
如果 ∣ P ∣ = 1 |P| =1 ∣P∣=1,那么过程返回 ∞ \infty ∞。如果 ∣ P ∣ = 2 |P| = 2 ∣P∣=2,那么返回这个点对的距离。
如果 ∣ P ∣ > 2 |P| > 2 ∣P∣>2,则进入分治过程。
分解
找出一条垂线段 L L L,使得 P P P在 L L L的左侧的 ∣ P l ∣ |P_l| ∣Pl∣等于右侧的 ∣ P r ∣ |P_r| ∣Pr∣。并且 P l P_l Pl都在L的左侧, P r P_r Pr都在 L L L的右侧。类似的,将数组 X X X和 Y Y Y分别划分成 X r X_r Xr X l X_l Xl Y r Y_r Yr Y l Y_l Yl。
解决
分别解决 P l P_l Pl X l X_l Xl Y l Y_l Yl和 P r P_r Pr X r X_r Xr Y r Y_r Yr的两个子问题,答案记为 A l A_l Al和 A r A_r Ar。并记 A = min ( A l , A r ) A = \min(A_l,A_r) A=min(Al,Ar)。
合并
主要问题是寻找最近的点对,分别位于 P l P_l Pl和 P r P_r Pr上。我们建立数组 Y ′ ⊆ Y Y' \subseteq Y Y′⊆Y,使得 Y ′ Y' Y′中的点的 x x x坐标到 L L L的距离都小于 A A A。并且保证 Y ′ Y' Y′中的元素还是按照 y y y坐标从小到大排序。紧接着,我们遍历 Y ′ Y' Y′中的元素,依次考虑之后的 8 8 8个点,即可。
代码
O ( n log n ) O(n \log n) O(nlogn)
#include <bits/stdc++.h>
#define FR freopen("in.txt", "r", stdin)
using namespace std;
typedef long long ll;
int n;
struct Point
{
double x;
double y;
bool operator<(const Point &p) const
{
return x < p.x;
}
} points[200005];
int X[200005];
int Y[200005];
int XT[200005];
int YT[200005];
bool cmpX(int a, int b)
{
return points[a].x < points[b].x;
}
bool cmpY(int a, int b)
{
return points[a].y < points[b].y;
}
const double INF = 10e19;
double dis(int a, int b)
{
return sqrt((points[a].x - points[b].x) * (points[a].x - points[b].x) + (points[a].y - points[b].y) * (points[a].y - points[b].y));
}
double merge(int l, int r)
{
if (l == r)
return INF;
if (l == r - 1)
return dis(l, r);
// split
int mid = (l + r) >> 1;
int lp = l;
int rp = mid + 1;
for (int i = l; i <= r; i++)
{
if (X[i] <= mid)
XT[lp++] = X[i];
else
XT[rp++] = X[i];
}
for (int i = l; i <= r; i++)
{
X[i] = XT[i];
}
lp = l;
rp = mid + 1;
for (int i = l; i <= r; i++)
{
if (Y[i] <= mid)
YT[lp++] = Y[i];
else
YT[rp++] = Y[i];
}
for (int i = l; i <= r; i++)
{
Y[i] = YT[i];
}
double al = merge(l, mid);
double ar = merge(mid + 1, r);
double as = min(al, ar);
lp = l;
rp = mid + 1;
int pt = l;
while (lp <= mid && rp <= r)
{
if (points[Y[lp]].y < points[Y[rp]].y)
YT[pt++] = Y[lp++];
else
YT[pt++] = Y[rp++];
}
while (lp <= mid)
{
YT[pt++] = Y[lp++];
}
while (rp <= r)
{
YT[pt++] = Y[rp++];
}
for (int i = l; i <= r; i++)
{
Y[i] = YT[i];
}
pt = l;
for (int i = l; i <= r; i++)
{
if (fabs((((double)points[Y[i]].x)) - points[mid].x) < as)
{
YT[pt++] = Y[i];
}
}
for (int i = l; i < pt; i++)
{
for (int j = i + 1; j <= i + 8 && j < pt; j++)
{
as = min(as, dis(YT[i], YT[j]));
}
}
return as;
}
int main()
{
scanf("%d", &n);
for (int i = 1; i <= n; i++)
{
scanf("%lf %lf", &points[i].x, &points[i].y);
X[i] = Y[i] = i;
}
sort(points + 1, points + 1 + n);
sort(X + 1, X + n + 1, cmpX);
sort(Y + 1, Y + n + 1, cmpY);
printf("%.4lf", merge(1, n));
return 0;
}
O ( n log 2 n ) O(n \log^2 n) O(nlog2n)
每次归并都内置排序
#include <bits/stdc++.h>
#define FR freopen("in.txt", "r", stdin)
using namespace std;
typedef long long ll;
int n;
struct Point
{
double x;
double y;
bool operator<(const Point &p) const
{
return x < p.x;
}
} points[200005];
int temp[200005];
const double INF = 10e19;
double dis(int a, int b)
{
return sqrt((points[a].x - points[b].x) * (points[a].x - points[b].x) + (points[a].y - points[b].y) * (points[a].y - points[b].y));
}
bool cmp(int a, int b)
{
return points[a].y < points[b].y;
}
double merge(int l, int r)
{
if (l == r)
return INF;
if (l == r - 1)
return dis(l, r);
// split
int mid = (l + r) >> 1;
double al = merge(l, mid);
double ar = merge(mid + 1, r);
double as = min(al, ar);
int pt = 0;
for (int i = l; i <= r; i++)
{
if (fabs(points[i].x - points[mid].x) < as)
{
temp[pt++] = i;
}
}
sort(temp, temp + pt, cmp);
for (int i = 0; i < pt; i++)
{
for (int j = i + 1; j <= i + 8 && j < pt; j++)
{
as = min(as, dis(temp[i], temp[j]));
}
}
return as;
}
int main()
{
scanf("%d", &n);
for (int i = 1; i <= n; i++)
{
scanf("%lf %lf", &points[i].x, &points[i].y);
}
sort(points + 1, points + 1 + n);
printf("%.4lf", merge(1, n));
return 0;
}