Closest Pair of Points Problem（最近点对问题）

最近点对问题介绍

在二维平面内的n个点中，找出（欧式）距离最近的两个点

问题解决的初步思路

蛮力法（Brute Force），实现过程大致为：

1. 遍历二维平面内的所有点，并计算任意两点之间的距离
2. 注意比较这些距离的大小，找出其中的最短距离，则这两点即为所求的"最近点对"

总共有 C(n, 2) 种组合，即 n×(n-1)/2，时间复杂度为 O(n²)

暴力枚举的伪代码:

蛮力法很容易想到，但算法的时间复杂度很高，这时就要引入一种算法设计思想 —— 分治策略(Divide and Conquer)

最近点对问题分治法初步分析

步骤: 设集合S用数组 S[1, …, n] 表示

1. 分解: 根据x坐标对S中的点升序排序，以经过点 S[n/2] 的垂直线L将S均分为两个子集:
   S_l = S[1, …, m] S_r = S[m+1, …, n]，其中，m = n / 2
2. 递归: 用同样的分治算法分别求 S_l 和 S_r 的最近点对的距离d_l 和d_r
3. 组合: 求 S_l 中的点和 S_r 中的点之间的最近距离 d’ ，则S中点对的最近距离为 min{d_l, d_r, d’}

分析:

如何求 d’ ?
一种最简单的方法就是把 S_l 中的所有点与 S_r 中的所有点计算距离，但时间复杂度太高

如何筛选出尽量少的中间点对?
我们知道组合步骤是在递归步骤之后，所以此时 d_l 和 d_r 已经被求出来了，所以我们在求 d’ 是可以利用 d_l 和 d_r
       令 d’ = min (d_l, d_r)，T_l 和 T_r 分别为 S_l 和 S_r 中与垂线的距离 ≤ d‘ 的点组成的子集，T = T_l ∪ T_r

这样子不就简化了算法吗，但隐隐之中感觉似乎漏掉了些什么…

接下来探讨最坏情况，如果非常多的点都是集中在分割线的附近，这时候仅仅只是排除掉了几个点而已，可以看到最坏情况下依然是 O(n²)

最近点对问题分治法进一步分析

首先我们要知道几个性质，看看能不能利用这些性质来设计高效的算法

假设如果存在 P_l ∈ S_l，P_r ∈ S_r，使得 dist(P_l, P_r) = D 且 D < d’ [ dist() 是求两点之间距离的函数]，于是 (P_l, P_r) 必定满足:

• P_l ∈ T_l，P_r ∈ T_r，即这两点必定分布在垂线两侧
       如果 (P_l, P_r) 同属于T_l或T_r 同一侧，即在T_l或T_r 一侧能够找到两个点，使得这两点的距离为D(<d’)，这与 S_l 和 S_r 的最近点对距离为 d_l 和 d_r 矛盾 (别忘了d’ = min{d_l, d_r})

• P_l，P_r位于以垂线为中轴线的 2d’ × d’ 的矩形区域内
       看我用红线连接的那两个点，没有框在矩形区域内，可以看见它们之间的距离一定会大于 d’，即 dist(P_l, P_r) > d’ ，与之前的假设 dist(P_l, P_r) = D 且 D < d’ 矛盾
  

• 以垂线为中轴线的 2d’ × d’ 的矩形区域中至多含T中的6个点，且 T_l 和 T_r 中各至多有4个
       我们取 2d’ × d’ 矩形的右边一半出来，即为 d’ × d’ 的小矩形，假设已经有四个点分布在这个矩形的四个角，那么矩形内部无法再有任何点存在，因为这个小矩形是 ∈ T_r，假设如果小矩形内还有一个红点，令该点距离其余四点距离均最远，那么该红点必定在小矩形的中心位置，距离其余四个点的距离为 (√2d’) / 2 < d’ ，我们要始终记得d’本身的含义，它是T_l，T_r 中的最小值，不可能有比d’还要小的距离出现，矛盾，所以小矩形至多只能放4个点，左边的小矩形也至多能放4个点，两个小矩形会有两个点重合，所以 2d’ × d’ 的矩形中至多只能有6个点
  

根据以上的性质，我们可以得出:

     1)如果把T中的点按照Y坐标升序排序，任取一点T[i]，那么 T[i + 7] 与 T[i] 的距离肯定超过 d’

     2)如果把T_l中的点按照Y坐标升序排列 [T_l也就是取矩形的左边部分]，任取一点T_l[i]，那么 T_r[i + 5] 与 T_l[i] 的距离肯定超过 d’。T_r 有同样的性质

分治策略的伪代码:

算法的时间复杂度不好说，如果你有去看其他博客，你会发现有两种版本，原因在于对y轴坐标排序是放函数外，还是放函数内的不同，使得时间复杂度有两种写法，两者造成的区别其实是在于n的数据集大小

---

题目集

题目链接	难度	算法简介
P1429 平面最近点对（加强版）	★☆☆☆☆	模板题
ZOJ2107 Quoit Design	★☆☆☆☆	模板题
POJ3714 Raid	★☆☆☆☆	模板题

---

P1429 平面最近点对（加强版）

代码

• 语言支持：c++

c++ Code:

#include<iostream>
#include<algorithm>
#include<cmath>
#include<cstdio>
#include<vector>
#include<climits>
using namespace std;
typedef struct node{
	int x, y;
}Point;
vector<Point> T(300005); // 放函数体外，不然每次都要动态分配如此巨大的数组
/*
求两点间距离
param  : 两个点坐标
return : 距离(double型)
*/
double dist(Point p1, Point p2)
{
	return sqrt((double)(p1.x - p2.x) * (p1.x - p2.x) + (double)(p1.y - p2.y) * (p1.y - p2.y));
}

/*
sort()里的比较函数,对x进行升序排序，若x相同则依照y进行升序排序
param  : 两个点
return : 两点横纵坐标比较大小后的结果
*/
bool cmpx(Point a, Point b)
{
	if (a.x == b.x)
		return a.y < b.y;
	return a.x < b.x;
}

/*
sort()里的比较函数,对y进行排序
param  : 两个点
return : 两点纵坐标比较大小后的结果
*/
bool cmpy(Point a, Point b)
{
	return a.y < b.y;
}

/*
分治法求最近点对
param  : 点的数组，low, high 两个指针/标记
return : 最近点对的距离
*/
double Merge(vector<Point>& P, int low, int mid, int high, double d)
{
	int k = 0;
	for (int i = low; i <= high; i++){
		if (fabs(1.0 * P[mid].x - 1.0 * P[i].x) < d)
			T[k++] = P[i];
	}
	// k是T的大小
	sort(T.begin(), T.begin() + k, cmpy);
	double min_d = d;
	for (int i = 0; i < k - 1; i++){ // 边界是k - 1
		for (int j = i + 1; j < min(i + 7, k); j++){
			if (dist(T[i], T[j]) < min_d)
				min_d = dist(T[i], T[j]);
		}
	}
	return min_d;
}

double cp_merge_solution(vector<Point>& P, int low, int high)
{
	// 如果自己与自己的距离最近，那直接cout << 0.0000 << endl 就结束了
	if (low >= high) return (double)INT_MAX;
	if (low + 1 == high) // 这两点中间没有任何点，所有直接求距离即可
		return dist(P[low], P[high]);

	int mid = (high + low) / 2; // 分解

	// 分别获得两分区中最近点对距离，并取这两个距离之中的最小值
	double d = min(cp_merge_solution(P, low, mid), cp_merge_solution(P, mid + 1, high));
	return Merge(P, low, mid, high, d);
}

int main(){
	int N;
	scanf("%d", &N);
	vector<Point> P(N);
	for (int i = 0; i < N; i++)
		scanf("%d%d", &P[i].x, &P[i].y);
	// Divide and Merge
	sort(P.begin(), P.begin() + N, cmpx);
	printf("%.4lf\n", cp_merge_solution(P, 0, N - 1));
	return 0;
}

---

ZOJ上也有一道 最近点对问题 的模板题

ZOJ2107 Quoit Design

代码

• 语言支持：c++

c++ Code:

#include<iostream>
#include<algorithm>
#include<cmath>
#include<cstdio>
#include<vector>
#include<climits>
using namespace std;
typedef struct node{
	double x, y;
}Point;
vector<Point> T(100005);
double dist(Point p1, Point p2)
{
	return sqrt((p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y));
}

bool cmpx(Point a, Point b)
{
	if (a.x == b.x)
		return a.y < b.y;
	return a.x < b.x;
}

bool cmpy(Point a, Point b)
{
	return a.y < b.y;
}

double Merge(vector<Point>& P, int low, int mid, int high, double d)
{
	int k = 0;
	for (int i = low; i <= high; i++){
		if (fabs(P[mid].x - P[i].x) < d)
			T[k++] = P[i];
	}
	// k是T的大小
	sort(T.begin(), T.begin() + k, cmpy);
	double min_d = d;
	for (int i = 0; i < k - 1; i++){ // 边界是k - 1
		for (int j = i + 1; j < min(i + 7, k); j++){
			if (dist(T[i], T[j]) < min_d)
				min_d = dist(T[i], T[j]);
		}
	}
	return min_d;
}

double cp_merge_solution(vector<Point>& P, int low, int high)
{
	if (low >= high) return (double)INT_MAX;
	if (low + 1 == high)
		return dist(P[low], P[high]);

	int mid = (high + low) / 2; // 分解

	double d = min(cp_merge_solution(P, low, mid), cp_merge_solution(P, mid + 1, high));
	return Merge(P, low, mid, high, d);
}

int main(){
	int N;
	while (scanf("%d", &N) && N){
		vector<Point> P(N);
		for (int i = 0; i < N; i++)
			scanf("%lf%lf", &P[i].x, &P[i].y);

		sort(P.begin(), P.begin() + N, cmpx);
		printf("%.2lf\n", cp_merge_solution(P, 0, N - 1) / 2.0);
	}
	return 0;
}

---

POJ3714 Raid

这题同样也是模板题，只有当两个点，一个是特工坐标，另一个是据点坐标才算两点距离，否则直接返回INT_MAX，很简单，但当时死活TLE…，有一些小小的细节，需要注意的点在代码中已经写了出来

代码

• 语言支持：c++

c++ Code:

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cmath>
#include<vector>
#include<climits>
using namespace std;
typedef long long ll;
struct Point{
	ll x, y;
	bool flag;
};
vector<Point> T(200004); // 放函数体外，否则会TLE
double dist(Point p, Point q)
{
	return sqrt(1.0 * (p.x - q.x) * (p.x - q.x) + 1.0 * (p.y - q.y) * (p.y - q.y));
}

bool cmpx(Point p, Point q)
{
	if (p.x == q.x)
		return p.y < q.y;
	return p.x < q.x;
}

bool cmpy(Point p, Point q)
{
	return p.y < q.y;
}

double Merge(vector<Point>& P, int low, int high)
{
	if (low >= high) return (double)INT_MAX;
	if (low + 1 == high){
		if (P[low].flag == P[high].flag) return (double)INT_MAX; // 判断两点是否不属于同一类 
		else return dist(P[low], P[high]);
	}
	
	int mid = (low + high) / 2;
	double d = min(Merge(P, low, mid), Merge(P, mid + 1, high));
	
	int cnt = 0;
	for (int i = low; i <= high; i++){
		if (fabs(1.0 * P[mid].x - 1.0 * P[i].x) < d){
			T[cnt++] = P[i];
		}
	}

	double min_d = d;
	sort(T.begin(), T.begin() + cnt, cmpy);
	for (int i = 0; i < cnt - 1; i++){ // 边界是 cnt - 1 
		for (int j = i + 1; j < min(i + 7, cnt); j++){
			if(T[i].flag != T[j].flag && dist(T[i], T[j]) < min_d) // 判断两点是否不属于同一类 
				min_d = dist(T[i], T[j]);
		}
	}
	return min_d;
}

int main(){
	int T, N;
	scanf("%d", &T);
	while (T--){
		scanf("%d", &N);
		vector<Point> P(2 * N + 3);
		for (int i = 0; i < N; i++){
			scanf("%lld%lld", &P[i].x, &P[i].y);	
		}
		for (int i = N; i < 2 * N; i++){
			scanf("%lld%lld", &P[i].x, &P[i].y);
			P[i].flag = 1;
		}
		sort(P.begin(), P.begin() + 2 * N, cmpx); // 这里记得是2 * N 
		printf("%.3lf\n", Merge(P, 0, 2 * N - 1));
	}
	return 0;
}