算法分析与设计——最近点对问题分治算法

---

一、问题

最近点对问题要求在一个包含n个点的集合中，找出距离最近的两个点，本题只考虑二维平面上的最近点对问题。

二、解析

假设点在二维平面上以笛卡尔直角坐标系给出坐标pi(xi,yi)，那么两点间的距离是标准欧几里得距离，即d(pi,pj)=sqrt((xi-xj)2+(yi-yj)2)；

最直接的思路是暴力求解所有点对之间的距离，找出最小的那一对，当然d(pi,pj)= d(pj,pi)，而且平方根函数是单调增的，所以我们可以只计算i<j的那部分点对的距离函数值不开根的那个值，即(xi-xj)2+(yi-yj)2，对程序进行优化；

但更好的思路是利用分治思想，此处引用《算法分析与设计基础》第三版Anany Levitin著，潘彦译，清华大学出版社P149-151

三、设计

输入：n个(n>=2)点的二维笛卡尔坐标
输出：最近点对以及他们之间的距离

1. 如果n=2，则直接计算两点间距离，返回该值，算法结束；
2. 如果n=3，则比较三点之间的最小距离，返回该值和该点对（n=3划分后必有n=1，会对递归造成困难）；
3. 按x坐标/y坐标进行区域划分：以各点x坐标的中位数为分割线划分成两个区域area1和area2；
4. 分别计算area1和area2的点集中的最近点对距离，对每个area都可以进行第3，4步的递归，例如area1可以根据中位数再分割为area1.1和area1.2，直到点集内点的个数<=3则进入1，2步求出结果；
5. 假设递归至最小的area1.1和area1.2中的最近对距离分别为d1和d2，合并两个区域，找出合并后的最近对距离，这个值一定在d1、d2或者在两个区域相分割线相邻两侧的点对之中；
6. 讨论合并两个区域时最近对是哪一个值，检查思路在解析中已经说明；
7. 合并成大区域之后，继续返回，直到整个点集都被合并完成。

四、分析

时间复杂度：排序算法为sort快排，复杂度为O(nlogn)，二分的分治策略是O(logn)，分治内的归并过程由于每个点都只扫描不超过6个点，为O(1)，遍历n个点，所以内部复杂度不超过O(n)，故整体不超过O(nlogn);

五、代码

#include<iostream>
#include<vector>
#include<cmath>
#include<time.h>
#include<cstring>
#include<algorithm>
 
using namespace std;
const int maxn = 5;//定义点的数量
const int mod = 10;//定义点值的范围

struct Point{//点集
	double x;
	double y;
};
struct PointPair{//记录最近对的可能值
    Point a, b;
    double d;
};

bool Compx(const Point &p1, const Point &p2)
{
	return p1.x < p2.x;
}
bool Compy(const Point &p1, const Point &p2)
{
	return p1.y < p2.y;
}
double Distance(const Point &a, const Point &b)
{
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
//归并左右两个点集
void Merge(vector<Point> &v, int left, int mid, int right)
{
	vector<Point> vl(v.begin()+left,v.begin()+mid);
	int i = left;
	int j = 0;
	int k = mid;
	while(i < right){
		if(j < vl.size() && (k == right || vl[j].y <= v[k].y))
			v[i++] = vl[j++];
		if(k < right && (j == vl.size() || vl[j].y > v[k].y))
			v[i++] = v[k++];
	}
}

//随机生成点
void CreatePoints(Point Points[], int pointNumber){
	srand(time(NULL));//随机化随机数种子
    for (int i = 0; i<pointNumber; i++){
		Points[i].x = rand()%mod;
		Points[i].y = rand()%mod;
        //生成带2位小数的点值
        //Points[i].x = rand()%(mod*100)/100.0;
		//Points[i].y = rand()%(mod*100)/100.0;
		//每1000个数据乘以一个特定的数，将数据尽量分散，减少重复
		Points[i].x *= ((i / 1000) + 1);
		Points[i].y *= ((i / 1000) + 1);
		//遍历已经生成的所有点，一旦发现重合则重新生成该点
		for (int j = 0; j < i; j++){
			if (Points[j].x == Points[i].x && Points[j].y == Points[i].y) {
				i--;
				break;
			}
		}
	}
}

PointPair Closest(vector<Point> &vx, int left, int right)
{
    PointPair ans;
    //当点集内只有2-3个点时特殊考虑，直接用暴力解答
    if(right - left == 2){
		if(vx[left].y > vx[right-1].y){
			swap(vx[left],vx[right-1]);
		}
        ans.a= vx[left];
        ans.b= vx[right-1];
        ans.d=Distance(vx[left],vx[right-1]);
        return ans;
    }
	if(right - left == 3){
		sort(vx.begin()+left,vx.begin()+right,Compy); 
		double d1 = Distance(vx[left],vx[left+1]);
		double d2 = Distance(vx[left],vx[right-1]);
		double d3 = Distance(vx[left+1],vx[right-1]);
        //记录最近对
        if(min({d1,d2,d3})==d1){
            ans.a = vx[left];
            ans.b = vx[left+1];
        }
        else if(min({d1,d2,d3})==d2){
            ans.a = vx[left];
            ans.b = vx[right-1]; 
        }
        else if(min({d1,d2,d3})==d3){
            ans.a = vx[left+1];
            ans.b = vx[right-1];
        }
        //返回最近距离
		ans.d=min({d1,d2,d3});
        return ans;
    }
    //若点集超过3个，要用递归进行分治
	int mid = (left + right) / 2; 
	double mx = vx[mid].x;
	PointPair pl = Closest(vx,left,mid);
	PointPair pr = Closest(vx,mid,right);
	if(pl.d<pr.d)
        ans = pl;
    else
        ans = pr;
    double d = min(pl.d, pr.d);
    //查找两个区间的中间部分是否存在最小值
    Merge(vx,left,mid,right);
	vector<Point> vp;
	for(int i = left; i < right; ++i){
		if(abs(vx[i].x - mx) < d)
			vp.push_back(vx[i]);
	}
	for(int i = 0; i < vp.size(); ++i){
		for(int j = i + 1; j < i+1+6 && vp.size(); ++j){//可以只遍历相邻六个点
			if(vp[j].y - vp[i].y >= d)
				break;
			else{
				double dm = Distance(vp[i],vp[j]);
				if(dm < d){
                    d = dm;
                    ans.a = vp[i];
                    ans.b = vp[j];
                }	
			}
		}
	}
    ans.d = d;
    return ans;
}

int main()
{
	vector<Point> v;
    Point points[maxn];
    CreatePoints(points, maxn);
    for (int i = 0; i < maxn;i++){
		v.push_back(points[i]);
        cout << "(" << points[i].x << "," << points[i].y << "); ";
    }
    cout << endl;
    sort(v.begin(),v.end(),Compx);
    PointPair ans = Closest(v, 0, v.size());
    cout << ans.d << endl;
    cout << "(" << ans.a.x << "," << ans.a.y << ") & (" << ans.b.x << "," << ans.b.y << ")" << endl;
    system("pause");
    return 0;
}