问题场景：在应用中，常用诸如点、圆等简单的几何对象代表现实世界中的实体。在涉及这些几何对象的问题中，常需要了解其邻域中其他几何对象的信息。例如，在空中交通控制问题中，若将飞机作为空间中移动的一个点来看待，则具有最大碰撞危险的2架飞机，就是这个空间中最接近的一对点。这类问题是计算几何学中研究的基本问题之一。

      问题描述：给定平面上n个点，找其中的一对点，使得在n个点的所有点对中，该点对的距离最小。严格地说，最接近点对可能多于1对。为了简单起见，这里只限于找其中的一对。

      1、一维最接近点对问题

       算法思路：

       这个问题很容易理解，似乎也不难解决。我们只要将每一点与其他n-1个点的距离算出，找出达到最小距离的两个点即可。然而，这样做效率太低，需要O(n^2)的计算时间。在问题的计算复杂性中我们可以看到，该问题的计算时间下界为Ω(nlogn)。这个下界引导我们去找问题的一个θ(nlogn)算法。采用分治法思想，考虑将所给的n个点的集合S分成2个子集S1和S2，每个子集中约有n/2个点，然后在每个子集中递归地求其最接近的点对。在这里，一个关键的问题是如何实现分治法中的合并步骤，即由S1和S2的最接近点对，如何求得原集合S中的最接近点对，因为S1和S2的最接近点对未必就是S的最接近点对。如果组成S的最接近点对的2个点都在S1中或都在S2中，则问题很容易解决。但是，如果这2个点分别在S1和S2中，则对于S1中任一点p，S2中最多只有n/2个点与它构成最接近点对的候选者，仍需做n^2/4次计算和比较才能确定S的最接近点对。因此，依此思路，合并步骤耗时为O(n^2)。整个算法所需计算时间T(n)应满足:　T(n)=2T(n/2)+O(n^2)。它的解为T(n)=O(n^2)，即与合并步骤的耗时同阶，这不比用穷举的方法好。从解递归方程的套用公式法，我们看到问题出在合并步骤耗时太多。这启发我们把注意力放在合并步骤上。      

      设S中的n个点为x轴上的n个实数x1,x2,..,xn。最接近点对即为这n个实数中相差最小的2个实数。我们显然可以先将x1,x2,..,xn排好序，然后，用一次线性扫描就可以找出最接近点对。这种方法主要计算时间花在排序上，在排序算法已经证明，时间复杂度为O(nlogn)。然而这种方法无法直接推广到二维的情形。因此，对这种一维的简单情形，我们还是尝试用分治法来求解，并希望能推广到二维的情形。假设我们用x轴上某个点m将S划分为2个子集S1和S2，使得S1={x∈S|x≤m}；S2={x∈S|x>m}。这样一来，对于所有p∈S1和q∈S2有p<q。递归地在S1和S2上找出其最接近点对{p1,p2}和{q1,q2}，并设d=min{|p1-p2|,|q1-q2|}，S中的最接近点对或者是{p1,p2}，或者是{q1,q2}，或者是某个{p3,q3}，其中p3∈S1且q3∈S2。如图所示。

     如果S的最接近点对是{p3,q3}，即|p3-q3|<d，则p3和q3两者与m的距离不超过d，即|p3-m|<d，|q3-m|<d，也就是说，p3∈(m-d,m]，q3∈(m,m+d]。由于在S1中，每个长度为d的半闭区间至多包含一个点（否则必有两点距离小于d），并且m是S1和S2的分割点，因此(m-d,m]中至多包含S中的一个点。同理，(m,m+d]中也至多包含S中的一个点。由图可以看出，如果(m-d,m]中有S中的点，则此点就是S1中最大点。同理，如果(m,m+d]中有S中的点，则此点就是S2中最小点。因此，我们用线性时间就能找到区间(m-d,m]和(m,m+d]中所有点，即p3和q3。从而我们用线性时间就可以将S1的解和S2的解合并成为S的解。也就是说，按这种分治策略，合并步可在O(n)时间内完成。这样是否就可以得到一个有效的算法了呢？还有一个问题需要认真考虑，即分割点m的选取，及S1和S2的划分。选取分割点m的一个基本要求是由此导出集合S的一个线性分割，即S=S1∪S2 ，S1∩S2=Φ，且S1={x|x≤m}；S2={x|x>m}。容易看出，如果选取m=[max(S)+min(S)]/2，可以满足线性分割的要求。选取分割点后，再用O(n)时间即可将S划分成S1={x∈S|x≤m}和S2={x∈S|x>m}。然而，这样选取分割点m，有可能造成划分出的子集S1和S2的不平衡。例如在最坏情况下，|S1|=1，|S2|=n-1，由此产生的分治法在最坏情况下所需的计算时间T(n)应满足递归方程:
   T(n)=T(n-1)+O(n)
    它的解是T(n)=O(n^2)。这种效率降低的现象可以通过分治法中“平衡子问题”的方法加以解决。即通过适当选择分割点m，使S1和S2中有大致相等个数的点。自然地，我们会想到用S的n个点的坐标的中位数来作分割点。在选择算法中介绍的选取中位数的线性时间算法使我们可以在O(n)时间内确定一个平衡的分割点m。

       本程序确定平衡点采用m=[max(S)+min(S)]/2方法。如果需要利用中位数作分割点，看结合笔者博文《0005算法笔记——线性时间选择》改写。

     一维最接近临近点对问题程序清单如下：

//2d10-1 一维最邻近点对问题
#include "stdafx.h"
#include <ctime>
#include <iostream> 
using namespace std; 

const int L=100;
//点对结构体
struct Pair
{
	float d;//点对距离
	float d1,d2;//点对坐标
};
float Random();
int input(float s[]);//构造S
float Max(float s[],int p,int q);
float Min(float s[],int p,int q);
template <class Type>
void Swap(Type &x,Type &y);
template <class Type>
int Partition(Type s[],Type x,int l,int r);
Pair Cpair(float s[],int l,int r);

int main()
{
	srand((unsigned)time(NULL));
	int m;
	float s[L];
	Pair d;
	m=input(s);
	d=Cpair(s,0,m-1);
	cout<<endl<<"最近点对坐标为： (d1:"<<d.d1<<",d2:"<<d.d2<<")";
	cout<<endl<<"这两点距离为： "<<d.d<<endl;
	return 0;
}


float Random()
{
	float result=rand()%10000;
	 return result*0.01;
}

int input(float s[])
{
	int length;
	cout<<"输入点的数目： ";
	cin>>length;
	cout<<"点集在X轴上坐标为：";
	for(int i=0;i<length;i++)
	{
		s[i]=Random();
		cout<<s[i]<<" ";
	}
	
	return length;
}


float Max(float s[],int l,int r)//返回s[]中的最大值
{
	float s_max=s[l];
	for(int i=l+1;i<=r;i++)
		if(s_max<s[i])
			s_max=s[i];
	return s_max;
}

float Min(float s[],int l,int r)//返回s[]中的最小值
{
	float s_min=s[l];
	for(int i=l+1;i<=r;i++) 
		if(s_min>s[i])
			s_min=s[i];
	return s_min;
}

template <class Type>
void Swap(Type &x,Type &y)
{
	Type temp = x;
	x = y;
	y = temp;
}

template <class Type>
int Partition(Type s[],Type x,int l,int r)
{
	int i = l - 1,j = r + 1;

	while(true)
	{
		while(s[++i]<x && i<r);
		while(s[--j]>x);
		if(i>=j)
		{
			break;
		}
		Swap(s[i],s[j]);
	}
	return j;
}

//返回s[]中的具有最近距离的点对及其距离
Pair Cpair(float s[],int l,int r)
{
	Pair min_d={99999,0,0};//最短距离

	if(r-l<1) return min_d;
	float m1=Max(s,l,r),m2=Min(s,l,r);

	float m=(m1+m2)/2;//找出点集中的中位数

	//将点集中的各元素按与m的大小关系分组
	int j = Partition(s,m,l,r);

	Pair d1=Cpair(s,l,j),d2=Cpair(s,j+1,r);//递归
	float p=Max(s,l,j),q=Min(s,j+1,r);

	//返回s[]中的具有最近距离的点对及其距离
	if(d1.d<d2.d)
	{
		if((q-p)<d1.d)
		{
			min_d.d=(q-p);
			min_d.d1=q;
            min_d.d2=p;
			return min_d;
		}
		else return d1;
	}
	else
	{
		if((q-p)<d2.d)
		{
			min_d.d=(q-p);
			min_d.d1=q;
            min_d.d2=p;
			return min_d;
		}
		else return d2;
	}
}

程序运行结果如下：

      该算法的分割步骤和合并步骤总共耗时O(n)。因此，算法耗费的计算时间T(n)满足递归方程：

       解此递归方程可得T(n)=O(nlogn)。

         2、二维最接近点对问题

     将以上过程推广到二维最接近点对问题，设S中的点为平面上的点，它们都有2个坐标值x和y。为了将平面上点集S线性分割为大小大致相等的2个子集S1和S2，我们选取一垂直线l:x=m来作为分割直线。其中m为S中各点x坐标的中位数。由此将S分割为S1={p∈S|px≤m}和S2={p∈S|px>m}。从而使S1和S2分别位于直线l的左侧和右侧，且S=S1∪S2 。由于m是S中各点x坐标值的中位数，因此S1和S2中的点数大致相等。递归地在S1和S2上解最接近点对问题，我们分别得到S1和S2中的最小距离d1和d2。现设d=min(d1,d2)。若S的最接近点对(p,q)之间的距离d(p,q)<d则p和q必分属于S1和S2。不妨设p∈S1，q∈S2。那么p和q距直线l的距离均小于d。因此，我们若用P1和P2分别表示直线l的左边和右边的宽为d的2个垂直长条，则p∈S1，q∈S2，如图所示:

     距直线l的距离小于d的所有点

       在一维的情形，距分割点距离为d的2个区间(m-d,m](m,m+d]中最多各有S中一个点。因而这2点成为唯一的末检查过的最接近点对候选者。二维的情形则要复杂些，此时，P1中所有点与P2中所有点构成的点对均为最接近点对的候选者。在最坏情况下有n2/4对这样的候选者。但是P1和P2中的点具有以下的稀疏性质，它使我们不必检查所有这n^2/4对候选者。考虑P1中任意一点p,它若与P2中的点q构成最接近点对的候选者，则必有d(p,q)<d。满足这个条件的P2中的点有多少个呢？容易看出这样的点一定落在一个d×2d的矩形R中，如下图所示:

包含点q的dX2d矩形R

     由d的意义可知P2中任何2个S中的点的距离都不小于d。由此可以推出矩形R中最多只有6个S中的点。事实上，我们可以将矩形R的长为2d的边3等分，将它的长为d的边2等分，由此导出6个（d/2）×（2d/3）的矩形。如左图所示:

矩阵R中点的稀疏性

     若矩形R中有多于6个S中的点，则由鸽舍原理易知至少有一个δ×2δ的小矩形中有2个以上S中的点。设u,v是这样2个点，它们位于同一小矩形中，则：

因此d(u,v)≤5d/6<d 。这与d的意义相矛盾。也就是说矩形R中最多只有6个S中的点。图4(b)是矩形R中含有S中的6个点的极端情形。由于这种稀疏性质，对于P1中任一点p，P2中最多只有6个点与它构成最接近点对的候选者。因此，在分治法的合并步骤中，我们最多只需要检查6×n/2=3n对候选者，而不是n^2/4对候选者。这是否就意味着我们可以在O(n)时间内完成分治法的合并步骤呢？现在还不能作出这个结论，因为我们只知道对于P1中每个S1中的点p最多只需要检查P2中的6个点，但是我们并不确切地知道要检查哪6个点。为了解决这个问题，我们可以将p和P2中所有S2的点投影到垂直线l上。由于能与p点一起构成最接近点对候选者的S2中点一定在矩形R中，所以它们在直线l上的投影点距p在l上投影点的距离小于d。由上面的分析可知，这种投影点最多只有6个。因此，若将P1和P2中所有S的点按其y坐标排好序，则对P1中所有点p，对排好序的点列作一次扫描，就可以找出所有最接近点对的候选者，对P1中每一点最多只要检查P2中排好序的相继6个点。

程序清单如下：

//2d10-2 二维最邻近点对问题
#include "stdafx.h"
#include<time.h>
#include<iostream> 
#include<cmath>

using namespace std;
const int M=50;

//用类PointX和PointY表示依x坐标和y坐标排好序的点
class PointX {
	public: 
		int operator<=(PointX a)const
		{ return (x<=a.x); }
		int ID; //点编号
		float x,y; //点坐标 
};

class PointY { 
	public: 
		int operator<=(PointY a)const
		{ return(y<=a.y); }
		int p; //同一点在数组x中的坐标 
		float x,y; //点坐标
};

float Random();
template <class Type>
float dis(const Type&u,const Type&v); 

bool Cpair2(PointX X[], int n,PointX& a,PointX& b, float& d);
void closest(PointX X[],PointY Y[],PointY Z[], int l, int r,PointX& a,PointX& b,float& d);

template <typename Type> 
void Copy(Type a[],Type b[], int left,int right);

template <class Type>
void Merge(Type c[],Type d[],int l,int m,int r);

template <class Type>
void MergeSort(Type a[],Type b[],int left,int right);

int main()
{ 
	srand((unsigned)time(NULL));
	int length;

	cout<<"请输入点对数：";
	cin>>length;

	PointX X[M];
	cout<<"随机生成的二维点对为："<<endl;

	for(int i=0;i<length;i++)
	{
		X[i].ID=i;
		X[i].x=Random();
		X[i].y=Random();
		cout<<"("<<X[i].x<<","<<X[i].y<<") ";
	}

	PointX a; 
	PointX b; 
	float d;

	Cpair2(X,length,a,b,d); 

	cout<<endl;
	cout<<"最邻近点对为：("<<a.x<<","<<a.y<<")和("<<b.x<<","<<b.y<<") "<<endl;
	cout<<"最邻近距离为： "<<d<<endl;

	return 0;
}

float Random()
{
	float result=rand()%10000;
	return result*0.01;
}

//平面上任意两点u和v之间的距离可计算如下
template <class Type>
inline float dis(const Type& u,const Type& v)
{
	float dx=u.x-v.x;
	float dy=u.y-v.y; 
	return sqrt(dx*dx+dy*dy); 
}

bool Cpair2(PointX X[], int n,PointX& a,PointX& b,float& d)
{
	if(n<2) return false;

	PointX* tmpX = new PointX[n];
	MergeSort(X,tmpX,0,n-1);

	PointY* Y=new PointY[n]; 
	for(int i=0;i<n;i++) //将数组X中的点复制到数组Y中
	{ 
		Y[i].p=i;
		Y[i].x=X[i].x;
		Y[i].y=X[i].y;
	} 

	PointY* tmpY = new PointY[n];
	MergeSort(Y,tmpY,0,n-1);

	PointY* Z=new PointY[n];
	closest(X,Y,Z,0,n-1,a,b,d); 

	delete []Y; 
	delete []Z;
	delete []tmpX;
	delete []tmpY;
	return true; 
}
void closest(PointX X[],PointY Y[],PointY Z[], int l, int r,PointX& a,PointX& b,float& d) 
{ 
	if(r-l==1) //两点的情形 
	{
		a=X[l];
		b=X[r];
		d=dis(X[l],X[r]);
		return; 
	} 

	if(r-l==2) //3点的情形 
	{
		float d1=dis(X[l],X[l+1]);
		float d2=dis(X[l+1],X[r]); 
		float d3=dis(X[l],X[r]); 

		if(d1<=d2 && d1<=d3) 
		{ 
			a=X[l];
			b=X[l+1];
			d=d1;
			return;
		} 

		if(d2<=d3)
		{ 
			a=X[l+1];
			b=X[r];
			d=d2;
		} 
		else { 
			a=X[l]; 
			b=X[r]; 
			d=d3; 
		}
		return; 
	} 

	//多于3点的情形，用分治法 
	int m=(l+r)/2; 
	int f=l,g=m+1; 

	//在算法预处理阶段，将数组X中的点依x坐标排序，将数组Y中的点依y坐标排序
	//算法分割阶段，将子数组X[l:r]均匀划分成两个不想交的子集，取m=(l+r)/2
	//X[l:m]和X[m+1:r]就是满足要求的分割。
	for(int i=l;i<=r;i++)
	{
		if(Y[i].p>m) Z[g++]=Y[i]; 
		else Z[f++]=Y[i];
	}

	closest(X,Z,Y,l,m,a,b,d);
	float dr;

	PointX ar,br;
	closest(X,Z,Y,m+1,r,ar,br,dr); 

	if(dr<d)
	{
		a=ar; 
		b=br; 
		d=dr; 
	} 

	Merge(Z,Y,l,m,r);//重构数组Y

	//d矩形条内的点置于Z中
	int k=l; 
	for(int i=l;i<=r;i++)
	{
		if(fabs(X[m].x-Y[i].x)<d)
		{ 
			Z[k++]=Y[i]; 
		}
	}

	//搜索Z[l:k-1]
	for(int i=l;i<k;i++) 
	{ 
		for(int j=i+1;j<k && Z[j].y-Z[i].y<d;j++) 
		{ 
			float dp=dis(Z[i],Z[j]);
			if(dp<d) 
			{ 
				d=dp; 
				a=X[Z[i].p];
				b=X[Z[j].p]; 
			}
		}
	} 
}

template <class Type>
void Merge(Type c[],Type d[],int l,int m,int r)
{
	int i = l,j = m + 1,k = l;
	while((i<=m)&&(j<=r))
	{
		if(c[i]<=c[j])
		{
			d[k++] = c[i++];
		}
		else
		{
			d[k++] = c[j++];
		}
	}

	if(i>m)
	{
		for(int q=j; q<=r; q++)
		{
			d[k++] = c[q];
		}	
	}
	else
	{
		for(int q=i; q<=m; q++)
		{
			d[k++] = c[q];
		}
	}
}

template <class Type>
void MergeSort(Type a[],Type b[],int left,int right)
{
	if(left<right)
	{
		int i = (left + right)/2;
		MergeSort(a,b,left,i);
		MergeSort(a,b,i+1,right);
		Merge(a,b,left,i,right);//合并到数组b
		Copy(a,b,left,right);//复制回数组a		
	}
}

template <typename Type> 
void Copy(Type a[],Type b[], int left,int right)
{
	for(int i=left;i<=right;i++) 
		a[i]=b[i]; 
}

程序结果：