分治法求最近点对

Dvide and Conquer

Implement the algorithm for the closest pair problem in your favourite language.

INPUT: n points in a plane.

OUTPUT: The pair with the least Euclidean distance.

算法的思想：

首先对所有的点按照x坐标进行从小到大排序，然后利用二分法进行分割，直到区间范围变为2 or 3（即只有2个点或只有3个点），然后可以求这2个点或是3个点的距离（很小的子问题）。

对于每一部分子问题，求其中间mid左右各min范围内的点的距离，如图：

对于这幅图中，只需将绿色的点与右边3个红色的点比较即可（格子按1/2&距离分割，可以保证每个格子内只有1个点）。

对于这幅图，只需将绿色的点与左边6个点比较即可

这里，当我们按照y坐标进行排序后（不考虑x坐标），因此对于每一个点，只需要将其与附近的11个点比较即可。

算法的合并时间：边界附近的点进行排序需要O(nlogn)。

如果每一次迭代保持2个序列，一个按照x进行排序，一个按照y进行排序，像归并排序一样，将预先排序好的2个方向的序列合并起来，实际上只需要0(n)时间。

因此算法的时间复杂度为T（n）=2T（n）+O(n)=0（nlogn）。

input:

-5,-5
-4,-4
-3,-6
-3,-3
-2,-3
-1,-5
0,-3
1,-4
0.5,-3
0.5,-4

output:

#include<iostream>  
#include<string.h>
#include<math.h>
#include <fstream>
#include<iomanip>
#include<algorithm> 
using namespace std;
int k=0; 


struct Record
{
	int c1,c2;
	double ds;
}record[10000];


struct Points
{
	double x,y;
}*points;


double dis(int a,int b)
{
	double dx = points[a].x-points[b].x,dy = points[a].y-points[b].y;
	return sqrt(dx*dx+dy*dy);
}


double min(double a,double b)
{
	return a>b?b:a;
}

   
int cmpx(const Points &a,const Points &b)    
{  
	if(a.x<b.x)   
		return 1;    
	else  
		return 0;   
}  

int cmpy(const Points &a,const Points &b)    
{  
	if(a.y<b.y)   
		return 1;    
	else  
		return 0;   
}  

void storage(double c,int l,int r)
{	//记录最近点对，用于输出
	record[k].c1 = l;
	record[k].c2 = r;
	record[k++].ds = c;
}


double ClosestPair(int l, int r)
{
	if(1==r-l)
	{
		double d = dis(l,r);
		storage(d,l,r);
		return d;
	}
	if(2==r-l)
	{
		double ll = dis(l,l+1);
		double rr = dis(l+1,r);
		storage(ll,l,l+1);
		storage(rr,l+1,r);
		return min(ll,rr);
	}
	int mid = l+((r-l)>>1);
	double dl = ClosestPair(l,mid);
	double dr = ClosestPair(mid+1,r);
	double Min=min(dl,dr);
	int h1=mid,h2=mid;
	while((points[mid].x-points[h1].x<=Min && h1>=l) || (points[h2].x-points[mid].x<=Min && h2<=r))
	{	//得到mid附近2&距离内的区间范围
		h1--;
		h2++;
	}
	sort(&points[l],&points[r],cmpy); //将mid附近2&距离内的点按照y坐标从小到大排序
	for(int i=h1+1;i<h2;i++)
	{
		int k=(i+11)>h2?h2:(i+11);	//只需要比较附近11个点即可
		for(int j=i+1;j<k;j++)
		{
			double d=dis(i,j);
			if(d>=Min)
				break;
			Min = d;
			storage(d,i,j);
		}
	}
	return Min;
}


int main()
{
	ifstream infile("./points.txt", ios::in); 
    char buf[30];
	double (*data)[2];	
	char *subarr=NULL;
	char *delims = ",";
	int n;
	cout<<"Please input the number of points(eg.input 10 int this test):"<<endl;
	cin>>n;
	points = new struct Points[n];
	data = new double[n][2];
    if(!infile.is_open())
	{
		cout<<"Error opening file!";
		exit(1);
	}
	int i=0,j=0,p=0;
	while(infile.getline(buf,30)!=NULL)
	{
		j = 0;
		subarr = strtok(buf,delims);
		while(subarr!=NULL)
		{
			data[i][j] = atof(subarr);	
			subarr = strtok(NULL,delims);
			j++;
		}
		
		points[p].x = data[i][0];
		points[p++].y = data[i][1];
		i++;
		
	}
	infile.close();
	sort(&points[0],&points[n],cmpx); //按横坐标从小到大排序
	double value = ClosestPair(0,n-1);
	int x,y;
	for(i=0;i<k;i++)
	{
		if(record[i].ds==value)
		{
			
			x = record[i].c1;
			y = record[i].c2;
			cout<<"points("<<points[x].x<<","<<points[x].y<<")and("<<points[y].x<<","<<points[y].y<<")"<<endl;
		}
	}
	cout<<"have the shortest distance :"<<value<<endl;
    return 0;
}