G_Weber

分治法实现最近对算法

最近对问题：这个问题非常简单，在一个平面中有N个点，查找距离最近的两个点。用蛮力法解决最简单，两两匹配，找出最近的一对，不过用蛮力法只局限于点比较少的情况，因为蛮力法要达到O(N*N)的时间复杂度。如果用分治的算法来做可以达到O(nlogn) 的时间复杂度，实际上需要多一个O(nlogn)的时间来就行预排序。

这两天用分治的思想实现平面最近对算法，并同时写了一个蛮力查找的算法作比较，果然用分治的算法快了很多。写的过程主要使用STL的容器存储数据，而算法中使用到的排序算法也没有自己实现，就直接使用STL中的sort排序。写的过程都把自己的思路写在注释上面了。

一、代码实现：
算法声明文件：CloestCouple.h

//-----------------------------------------
// CloestCouple.h
// 最近对算法
// 使用分治法实现
// 洪庚伟 2010.01.03
// http://www.cnblogs.com/G_Weber
//-----------------------------------------
#ifndef _CLOESTCOUPLE_HGW_
#define _CLOESTCOUPLE_HGW_
#include<vector>
#include<algorithm>
#include<cmath>

//分治法计算最近对
class CCloestCouple
{
public:
	struct Point
	{
		int x;
		int y;

		Point():x(0),y(0){}
		Point(int ix,int iy)
			:x(ix),y(iy)
		{}

	};

	//最近对算法还需要对点按照x轴、y轴进行排序
	//直接使用stl算法排序
	//因此需要编写比较点的函数对象

	//按x轴非递减排序
	struct PointCmpX
	{
		bool operator()(const Point & lhr,const Point & rhr)
		{
			if(lhr.x < rhr.x)
				return true;
			else if(lhr.x > rhr.x)
				return false;
			else
				return lhr.y < rhr.y;
		}
	};

	//按y轴非递减排序
	struct PointCmpY
	{
		bool operator()(const Point & lhr,const Point & rhr)
		{
			if(lhr.y < rhr.y)
				return true;
			else if(lhr.y > rhr.y)
				return false;
			else 
				return lhr.x < rhr.x;
		}
	};

	typedef std::vector<Point> PointVector;

public:
	//给出平面的一个点的集合，分治找出最近的两个点
	//算法基本是在位运行
	void find(PointVector& Points,Point& from,Point& to);

	//计算两点距离的平方
	float distance(const Point &from,const Point &to);

	//另外编写了一个蛮力法找最近对，用来检查分治法的正确性
	//以及比较运行效率
	void ForceFind(PointVector& Points,Point& from,Point& to);

private:

	//对有序的点集合Points[begin,end)，递归返回最近对两个点from 和 to
	void _find(PointVector& Points,int begin,int end,Point& from,Point& to);

	//在left、right 两个集合中查找是否有距离小于
	//ldis，的两个分别分布于left和right的点对，
	//用他们的值改写 from to

	//由于分治找出最近对后，还有可能存在这样的最近对：
	//两个点位于不同的分治组合，因此要有这一步检查
	void compare(PointVector &left,PointVector &right,
		float dis,Point &from,Point &to);

};
#endif

还有的是算法的实现文件：CloestCouple.cpp

//-----------------------------------------
// CloestCouple.cpp
// 最近对算法
// 使用分治法实现
// 洪庚伟 2010.01.03
// http://www.cnblogs.com/G_Weber
//-----------------------------------------
#include "CloestCouple.h"


void CCloestCouple::find(PointVector &Points,Point &from,Point &to)
{
	if(Points.size() < 2)
		return;
	//先对输入点集合做按x轴的排序
	std::sort(Points.begin(),Points.end(),PointCmpX());

	//递归分治查找
	_find(Points,0,Points.size(),from,to);
}

void CCloestCouple::_find(PointVector &Points,int begin,int end,Point &from,Point &to)
{
	if(end - begin == 2)
	{
		//两个点，直接返回结果
		from = Points[begin];
		to = Points[begin+1];
	}
	else if(end - begin == 3)
	{
		//三个点 A、B、C

		//先找出 AB近 还是 AC近，保存在 (from,to)中
		from = Points[begin];
		to = Points[begin+1];
		float ab = distance(from,to);
		float ac = distance(from,Points[begin+2]);
		if( ab > ac )
		{
			to = Points[begin+2];
			ab = ac;//始终用ab保存最小距离
		}
		//再与BC作比较
		if(ab > distance(Points[begin+1],Points[begin+2]))
		{
			from = Points[begin+1];
			to = Points[begin+2];
		}
	}
	else /*if(end - begin >= 4)*/
	{
		//当前集合元素足够多，平分成两份分治处理
		int mid = begin + (end - begin) / 2;

		//假设这两个集合根据直线 x=c进行划分的
		float c = (Points[mid].x + Points[mid-1].x) / 2;

 		Point rfrom,rto;
		Point lfrom,lto;

		//递归找出两个集合中的最近对
		//在区间[begin,mid)中查找
		_find(Points,begin,mid,lfrom,lto);
		//在区间[mid,end)中查找
		_find(Points,mid,end,rfrom,rto);

		//计算较近的一对
		float ldis = distance(lfrom,lto);
		float rdis = distance(rfrom,rto);
		//总是用ldis,lfrom,lto保存最近的距离
		if(ldis > rdis)
		{
			lfrom = rfrom;
			lto = rto;
			ldis = rdis;
		}
		
		//在左集合中提取位于[c-d,c]中的点
		//在右集合中提取位于[c,c+d]中的点
		PointVector left;
		PointVector right;

		//dis是距离的平方
		float d = sqrt(ldis);

		int i;
		for(i = mid-1; i>=begin && c - Points[i].x <= d;i--)
			left.push_back(Points[i]);

		for(i = mid; i < end && Points[i].x - c <= d;i++)
			right.push_back(Points[i]);

		//在left、right 两个集合中查找是否有距离的平方小于
		//ldis，的两个分别分布于left和right的点对，
		//用他们的值改写 lfrom lto
		compare(left,right,ldis,lfrom,lto);
		//记录返回值然后返回
		from = lfrom;
		to = lto;
	}
}

float CCloestCouple::distance(const Point &from,const Point &to)
{
	return (from.x - to.x ) * (from.x - to.x) 
		+ (from.y - to.y) * (from.y - to.y);
}

void CCloestCouple::compare(PointVector &left,PointVector &right,
		float dis,Point &from,Point &to)
{
	//先对两个集合按照y轴排序
	std::sort(left.begin(),left.end(),PointCmpY());
	std::sort(right.begin(),right.end(),PointCmpY());

	//dis是距离的平方
	float d = sqrt(dis);

	//遍历左集合，逐个跟右集合中的点进行比较
	//考虑到集合的已按y轴排好许了
	//如果遇到right集合中的点y坐标比
	//left集合中当前点的y坐标大d，就可以结束了
	PointVector::iterator lit = left.begin();
	for(;lit!=left.end();lit++)
	{
		PointVector::iterator rit = right.begin();
		for(;rit!=right.end();rit++)
		{
			if(lit->y - d > rit->y)
				continue;
			if(lit->y + d < rit->y)
				break;

			//(x,y)为当前left集合正在处理的点
			//只有right中的点位于[y-d,y+d]
			//才有可能两个集合中的点距离小于d
			float tdis = distance(*lit,*rit);
			if(tdis < dis)
			{
				from = *lit;
				to = *rit;
				dis = tdis;//再次记录最近一对
			}
		}
	}
}

void CCloestCouple::ForceFind(PointVector& Points,Point& from,Point& to)
{
	//嗯···暴力查找，超简单的代码
	float dis=distance(*Points.begin(),*Points.end());
	PointVector::iterator lit,rit;
	lit = Points.begin();
	for(;lit!=Points.end();lit++)
	{
		rit=lit+1;
		for(;rit!=Points.end();rit++)
		{
			float newd = distance(*lit,*rit);
			if(dis > newd)
			{
				dis = newd;
				from = *lit;
				to = *rit;
			}
		}
	}
}

二、测试程序
写了个Win32窗口测试一下代码
1.首先是随便点几个点，测试结果分治法和蛮力法都是一样的，看不出区别：

2.自动生成5000个点的时候，这时分治法的效率已经开始表现出优势了

3.用100000个点进行测试的时候，分治法还能高效处理，蛮力法直接让程序用完CPU，只得中断进程

最近对算法演示程序下载

By：洪庚伟

出处：http://www.cnblogs.com/G_Weber