用分治法求解三维空间中的最近点对

原创 2012年03月24日 20:23:28

题目如下:

用分治法求解下面的问题

 

输入:

P=(p(1),p(2),…,p(n))为三维空间中n个不同的点,即

P(i)=(x(i), y(i), z(i)) ,1≤i  ≤n

输出:

距离最近的两点。


 

所有的过程与寻找二维空间中的最近点对类似(见算法导论第二版591页),只是在找Y内的最短距离时,需要考虑的紧随其后的点的数目不同。

(1)Divide:我们按照y值的大小来分割三维空间。平面y=middleY代表将点集P分割的中间曲面,得到两个集合pL和pR;

(2)Conquer:分别对pL和pR集合求最近点对,和最近的距离,lMinDistance和rMinDistance,然后将当前的最近距离设置为lMinDistance和rMinDistance的最小值currentMinDistance=min{lMinDistance,rMinDistance};

(3)Merge:建立一个数组middleP和middleArrayY,分别代表平面y=middleY两侧y值的距离不超过currentMinDistance的范围内所有的点和所有的y值的集合。我们在middleP中找最近的点对。可以证明,对于middleP中的每一个点,仅需要检查紧随其后的15个点即可。

证明:我们可以在平面y=middleY的两侧做两个边长为currentMinDistance的正方体,分别记为L1和R1。假定在某一级递归调用中,最近的点对des1属于pL,des2属于pR,则des1和des2的距离Distance(des1,des2)一定严格小于min{lMinDistance,rMinDistance}。点des1必在平面y=middleY上,或者在平面的左侧;点des2必在平面y=middleY上,或者在平面的右侧。由此可知,des1和des2在以y=middleY为中心,两个在平面上投影重合的正方体L1和R1围城的区域θ内。

下面证明,至多有16个点位于θ内。因为L1中的所有点之间的距离至少为currentMinDistance,L1属于PL,PL之间的所有点之间的距离至少为currentMinDistance,所以至多有8个点可能位于该正方形内。有8个点的情况,恰恰是八个点都是顶点时。假如有9个点位于正方形内的话,由正方体的图形结构可知,第九个点的距离必定会和前8个点其中的一个的距离小于currentMinDistance,与PL中所有点之间的距离至少为currentMinDistance矛盾。类似的,R1中也至多有8个点位于该正方形内。所以,P中至多有16个点位于区域θ内,就很容易看出,仅需检查数组middleP中每个点之后的15个点,(这15个点有7个已经在求pL中的最小距离时检查过了),所以其实只需要检查后8个点即可。(实现时,还是按照15个实现的。)

接着对算法流程的论述,求出middleY中的最小距离后,需要和上面的两部递归求的的最小距离(lMinDistance和rMinDistance)比较,得到的最小值即为算法最后输出的最小距离。

 

实现代码如下:

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>

using namespace std;

class Point
{
public:
	double x;
	double y;
	double z;
};

// [low, high]
int randomInteger(int low, int high)
{
	return low+rand()%(high-low+1);
}

void swap(double & a, double & b)
{
	double tmp=a;
	a=b;
	b=tmp;
}

int quickSortPartition(Point * p, int low, int high)
{
	int random=randomInteger(low, high);
	swap(p[random], p[low]);
	int i, j;
	i=low, j=high;
	double pivot=p[low].y;
	for(j=low; j<=high; j++)
	{
		if(pivot>p[j].y)
		{
			i=i+1;
			swap(p[i], p[j]);
		}
	}
	swap(p[i], p[low]);
	return i;
}


int partition(double * p, int len, int low, int high)
{
	int i, j;
	int random=randomInteger(low, high);
	swap(p[low], p[random]);
	i=low, j=high;
	double pivot=p[low];
	for(j=low; j<=high; j++)
	{
		if(pivot>p[j])
		{
			i=i+1;
		}
	}
	swap(p[random], p[low]);
	return i;
}

void quickSort(Point * p, int low, int high)
{
	if(low<high)
	{
		int par=quickSortPartition(p, low, high);
		quickSort(p, low, par-1);
		quickSort(p, par+1, high);
	}
}

double findMiddle(double * p, int len, int ith)
{
	if(len==1)
		return p[0];
	int random=partition(p, len, 0, len-1);
	if(random == ith)
	{
		return p[random];
	}
	else if(random < ith)
	{
		return findMiddle(p+random+1, len-random-1, ith-random-1);
	}
	else
	{
		return findMiddle(p, random, ith);
	}
}



double Distance(const Point & p1, const Point & p2)
{
	return sqrt( (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y) + (p1.z-p2.z)*(p1.z-p2.z) );
}

double findMinDistance(Point * p, double * arrayX, double * arrayY, double * arrayZ, int n, Point & des1, Point & des2)
{
	double minDistance=INT_MAX;
	if(n<=3)
	{
		for(int i=0; i<n; i++)
		{
			for(int j=0; j<n; j++)
			{
				if(j != i)
				{
					if(Distance(p[i], p[j])<minDistance)
					{
						minDistance=Distance(p[i], p[j]);
						des1=p[i], des2=p[j];
					}
				}
			}
			return minDistance;
		}
	}
	else
	{
		double middleY=findMiddle(arrayY, n, (int)floor((double)(n/2)));
		// partition P set into PL and PR according to the middle value of y
		Point * pL=new Point[n];
		Point * pR=new Point[n];
		int lIndex=0, rIndex=0;
		double * leftArrayX=new double[n];
		double * leftArrayY=new double[n];
		double * leftArrayZ=new double[n];
		double * rightArrayX=new double[n];
		double * rightArrayY=new double[n];
		double * rightArrayZ=new double[n];
		for(int i=0; i<n; i++)
		{
			if(p[i].y<middleY)
			{
				pL[lIndex]=p[i];
				leftArrayX[lIndex]=p[i].x;
				leftArrayY[lIndex]=p[i].y;
				leftArrayZ[lIndex]=p[i].z;
				lIndex++;
			}
			else
			{
				pR[rIndex]=p[i];
				rightArrayX[rIndex]=p[i].x;
				rightArrayY[rIndex]=p[i].y;
				rightArrayZ[rIndex]=p[i].z;
				rIndex++;
			}
		}
		// divide and conquer
		double lMinDistance=findMinDistance(pL, leftArrayX, leftArrayY, leftArrayZ, lIndex, des1, des2);
		double rMinDistance=findMinDistance(pR, leftArrayX, leftArrayY, leftArrayZ, rIndex, des1, des2);
		double currentMinDistance=lMinDistance<rMinDistance ? lMinDistance : rMinDistance;
		Point * middleP=new Point[n];
		double * middleArrayY=new double[n];
		int middleIndex=0;
		for(int i=0; i<n; i++)
		{
			if(fabs(p[i].y-middleY)<=currentMinDistance)
			{
				middleP[middleIndex]=p[i];
				middleArrayY[middleIndex]=p[i].y;
			}
		}
		for(int i=0; i<middleIndex; i++)
		{
			for(int j=1; j<=15; j++)
			{
				if(Distance(p[i], p[i+j])<currentMinDistance)
				{
					currentMinDistance=Distance(p[i], p[i+j]);
					des1=p[i];
					des2=p[i+j];
				}
			}
		}
		minDistance=currentMinDistance;
		delete [] pL;
		delete [] pR;
		delete [] leftArrayX;
		delete [] leftArrayY;
		delete [] leftArrayZ;
		delete [] rightArrayX;
		delete [] rightArrayY;
		delete [] rightArrayZ;
		delete [] middleP;
		delete [] middleArrayY;
		return minDistance;
	}
}

int main()
{
	Point * p=new Point[4];
	for(int i=0; i<4; i++)
	{
		p[i].x=4-i;
		p[i].y=4-i;
		p[i].z=4-i;
	}
	quickSort(p, 0, 3);
	double * arrayX, * arrayY, * arrayZ;
	arrayX=new double[4];
	arrayY=new double[4];
	arrayZ=new double[4];
	for(int i=0; i<4; i++)
	{
		arrayX[i]=p[i].x;
		arrayY[i]=p[i].y;
		arrayZ[i]=p[i].z;
	}
	for(int i=0; i<4; i++)
	{
		cout<<p[i].x<<" "<<p[i].y<<" "<<p[i].z<<endl;
	}
	double minDistance;
	Point des1, des2;
	minDistance=findMinDistance(p, arrayX, arrayY, arrayZ, 4, des1, des2);
	cout << "The minimal distance is:" << minDistance <<endl;
	cout << "The two points that have the minimal distance is:" << endl;
	cout << "("<<des1.x<<","<<des1.y<<","<<des1.z<<")"<<endl;
	cout << "("<<des2.x<<","<<des2.y<<","<<des2.z<<")"<<endl;
	getchar();
	return 0;
}


版权声明:本文为博主原创文章,未经博主允许不得转载。

相关文章推荐

0007算法笔记——【分治法】最接近点对问题

问题场景:在应用中,常用诸如点、圆等简单的几何对象代表现实世界中的实体。在涉及这些几何对象的问题中,常需要了解其邻域中其他几何对象的信息。例如,在空中交通控制问题中,若将飞机作为空间中移动的一个点来看...

分治法求解最近点距

1.问题描述:最近点对问题很简单,就是给定一堆点(这里是二位坐标下),求解最近的点的距离,该问题可以用穷举法求解,双重循环就够了,就不说了,主要看一下分治法的代码,代码也很简单,有注释。2.分治法求解...

分治法最近点对问题

在二维平面上的n个点中,如何快速的找出最近的一对点,就是最近点对问题。     一种简单的想法是暴力枚举每两个点,记录最小距离,显然,时间复杂度为O(n^2)。     在这里介绍一种时间...

算法 最近点对 分治法

  • 2014-12-12 17:00
  • 4.92MB
  • 下载

POJ 3714 最近点对问题 分治法

题意:station和agent分别有n(1 题解:http://blog.csdn.net/w397090770/article/details/7295797 这里讲的很好。 #i...
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:深度学习:神经网络中的前向传播和反向传播算法推导
举报原因:
原因补充:

(最多只允许输入30个字)