分治算法求最近点对-优快云博客

本文链接：https://blog.youkuaiyun.com/z444_579/article/details/52214145

问题: 寻找n >=2 个平面点集中两个最近点。

应用：交通控制中寻找两个最近的交通工具。
// 完整源码：http://download.youkuaiyun.com/detail/z444_579/9603999
传统蛮力搜索算法中，需要O(n*n)次搜索，本文介绍一种分治算法，运动时间为O(nlgn);算法步骤如下：
1. 输入P(原始点集), X(P按x坐标递增), Y(P按y坐标递增) 三个点集，点集个数为n. 若有 n<4, 则对其执行蛮力搜索(两两检查), 并返回，

否则进行Step2;

2. 将P 沿垂直线l(即横坐标)分为PL 和 PR, 它们个数均等, PL上的点均在分离直线l的左侧或直线上，而PR 均在直线l 的右侧或直线上,
点集X 被划分为XL 和XR, 分别包含PL 和 PR 的点，且XL(和XR) 每个的点集均按x 排序; 同理，将Y 分为YL 和 YR, 分别包含包含PL 和 PR 的

点，均按y 坐标排序(友情提示，一种简单的方法是从头到尾遍历Y 中点，若Yi 在XL中, 则该Yi 点属于YL, 依次类推，反归并排序!)，转Step3.

3. 进行两次递归调用，第一次找出PL 中最近点和最近距离deltaL，输入为PL, XL, YL; 第二次找出PR 中最近点和最近距离deltaR，输入为
PR, XR, YR; 转Step4. (提示, 因此，最近点(2个点) 只有3种情况：a.要么都在PL 中，b.要么在PR 中, c.要么是PL 和PR 中两个点, 其中情况c

必定位于分离直线l为中心，宽度为delta = min(deltaL, deltaR)的垂直带区域内);

4. 建立点集Y', 它是所有宽度为2*delta 垂直带中的点集, Y' 也是按y 坐标排序; 巧妙方法来了，对于Y' 中的每个点p, 寻找紧随其后(排序)

的7 个点, 计算p 到这7 个点的距离, 重复此步骤直到检查完Y' 中的所有点, 最近距离为delta';转Step5.

// 注，若看不懂，直接算也是对的! 简单解释，为什么选7个点, 因为第8 个点必定位于delta距离之外了! 见图b).
5. 若delta' < delta, 则垂直带区域中的点更近, 当前最近距离为delta';否则最近距离为delta.

// LookClosetPoint.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "time.h"
#include "stdlib.h"
#include "math.h"
#include "algorithm"

#include "assert.h"
#include "iostream"

#include "windows.h"

using namespace std;
// 定义点的最大个数;
#define MAX_SIZE 500
// 定义一个点;


class Point2
{
public:
	Point2()
	{
		m_x = 0.0;
		m_y = 0.0;
		m_i = 0;
	}

	double m_x;
	double m_y;
	int    m_i;
};

// 定义一个点集;
class PointsSet
{
public:
	PointsSet()
	{
		m_iLength = 0;
		m_points = NULL;
	}

	~PointsSet()
	{
		Clear();
	}

	void Clear()
	{
		if (m_points)
		{
			delete[] m_points;
			m_points = NULL;
		}
	}

	void Alloc(int n)
	{
		Clear();

		m_points = new Point2[n + 1];
		m_iSize  = n;
		m_iLength = 0;
	}

	Point2 *m_points;
	int     m_iSize;
	int     m_iLength;
};
// 根据X升序对点集进行排序;
void SortPointsByX(PointsSet *P, int n);
// 根据Y升序对点集进行排序;
void SortPointsByY(PointsSet *P, int n);

// 计算两点距离;
double Compute2PointSqDist(double x1, double y1, double x2, double y2)
{
	double dTemp = 0.0;
	dTemp = pow(x2 - x1, 2) + pow(y2 - y1, 2);
	return dTemp;
}

// 蛮力法计算点集中距离最小的2个点;
double  ComputePointsSetSqDistBrute(const PointsSet& P, int &iClosetIdx1, int &iClosetIdx2, double dMinSqDist)
{
	int n = P.m_iLength;

	if (n < 2)
	{
		assert(n > 1);
		return -10.0;
	}

	for (int i=0; i<n; i++)
	{
		for (int j=i+1; j<n; j++)
		{
			double dTemp = Compute2PointSqDist(P.m_points[i].m_x, P.m_points[i].m_y, \
				                               P.m_points[j].m_x, P.m_points[j].m_y);

			if (dTemp < dMinSqDist)
			{
				dMinSqDist = dTemp;
				iClosetIdx1 = i;
				iClosetIdx2 = j;
			}
		}
	}

	return dMinSqDist;
}

// 该索引是否在点集中;
bool IsIndexInPointsSet(const PointsSet &P, int iIndex)
{
	for (int i=0; i<P.m_iLength; i++)
	{
		if (P.m_points[i].m_i == iIndex)
		{
			return true;
		}
	}

	return false;
}

// 计算点集中最近的两个点;
double ComputeClosetPoint(const PointsSet& P, const PointsSet&X, const PointsSet&Y, int iClosetIdx1, int iClosetIdx2, double dMinSqDistQuick)
{
	int n = P.m_iLength;

	if (n < 1)
	{
		// 错误！
		return -10.0;
	}

	// 如果小于4个点则蛮力计算;
	if (n < 4)
	{
		return ComputePointsSetSqDistBrute(P, iClosetIdx1, iClosetIdx2, dMinSqDistQuick);
	}
	// 将点集按照X 升序分为XL, XR;
	PointsSet XL, XR; 
	
	int iSplit = n / 2;
	XL.Alloc(iSplit + 1);
	XR.Alloc(iSplit + 1);

	// 将X分为XL, XR;
	for (int i=0; i<=iSplit; i++)
	{
		XL.m_points[XL.m_iLength++] = X.m_points[i]; 
	}

	for (int i=iSplit; i<n; i++)
	{
		XR.m_points[XR.m_iLength++] = X.m_points[i]; 
	}
	// X 分割位置;
	double dSplitX = XL.m_points[XL.m_iLength - 1].m_x;

	// PL中的点与XL 一样, PR 中的点与XR 一致，但没有排序;
	PointsSet PL, PR; 
	PL.Alloc(iSplit + 1);
	PR.Alloc(iSplit + 1);

	for (int i=0; i<n; i++)
	{
		if (IsIndexInPointsSet(XL, P.m_points[i].m_i))
		{
			PL.m_points[PL.m_iLength++] = P.m_points[i];
		}

		if (IsIndexInPointsSet(XR, P.m_points[i].m_i))
		{
			PR.m_points[PR.m_iLength++] = P.m_points[i];
		}

		/*else
		{
			MessageBox(NULL, L"运行错误！", NULL, MB_ICONINFORMATION);
			break;
		}*/
	}

	// YL中的点与XL 一样, YR 中的点与XR 一致，但按照Y排序;
	PointsSet YL, YR;
	YL.Alloc(iSplit + 1);
	YR.Alloc(iSplit + 1);

	for (int i=0; i<n; i++)
	{
		if (IsIndexInPointsSet(XL, Y.m_points[i].m_i))
		{
			YL.m_points[YL.m_iLength++] = Y.m_points[i];
		}

		if (IsIndexInPointsSet(XR, Y.m_points[i].m_i))
		{
			YR.m_points[YR.m_iLength++] = Y.m_points[i];
		}

	/*	else
		{
			MessageBox(NULL, L"运行错误！", NULL, MB_ICONINFORMATION);
			break;
		}*/
	}

	if (PL.m_iLength != YL.m_iLength)
	{
		MessageBox(NULL, L"运行错误！", NULL, MB_ICONINFORMATION);
		return 0.0;
	}

	if (PR.m_iLength != YR.m_iLength)
	{
		MessageBox(NULL, L"运行错误！", NULL, MB_ICONINFORMATION);
		return 0.0;
	}

	int iClosetLIdx1 = 0, iClosetLIdx2 = 0;
	double dDeltaL = ComputeClosetPoint(PL, XL, YL, iClosetLIdx1, iClosetLIdx2, dMinSqDistQuick);
//	dMinSqDistQuick = min(dMinSqDistQuick, dDeltaL);

	PL.Clear();
	XL.Clear();
	YL.Clear();

	int iClosetRIdx1 = 0, iClosetRIdx2 = 0;
	double dDeltaR = ComputeClosetPoint(PR, XR, YR, iClosetRIdx1, iClosetRIdx2, dMinSqDistQuick);
//	dMinSqDistQuick = min(dMinSqDistQuick, dDeltaR);

	PR.Clear();
	XR.Clear();
	YR.Clear();

	double dDaltaV = min( min(dDeltaL, dDeltaR), dMinSqDistQuick);

	//double dDaltaV = min(dDeltaL, dDeltaR);

	if (dDaltaV < 0.0)
	{
		MessageBox(NULL, L"运行错误！", NULL, MB_ICONINFORMATION);
		return dDaltaV;
	}

	PointsSet Yv;
	Yv.Alloc(Y.m_iLength);
	
	for (int i=0; i<Y.m_iLength; i++)
	{
		double dTempX = Y.m_points[i].m_x;
		if (dTempX >= dSplitX - dDaltaV && dTempX <= dSplitX + dDaltaV)
		{
			Yv.m_points[Yv.m_iLength++] = Y.m_points[i];
		}
	}

	double dMinSqDist = dDaltaV;

	for (int i=0; i<Yv.m_iLength; i++)
	{
		for (int j=i+1; j<i+8 && j < Yv.m_iLength; j++)
		{
			double dTemp = Compute2PointSqDist(Yv.m_points[i].m_x, Yv.m_points[i].m_y, Yv.m_points[j].m_x, Yv.m_points[j].m_y);

			if (dTemp < dMinSqDist)
			{
				dMinSqDist = dTemp;
				iClosetIdx1 = i;
				iClosetIdx2 = j;
			}
		}
	}

	return dMinSqDist;
}

// X升序对点集排序;
void SortPointsByX(PointsSet *P, int n)
{
	int k = 0;

	for (int i=0; i<n; i++)
	{
		k = i;

		for (int j=i+1; j<n; j++)
		{
			if (P->m_points[j].m_x < P->m_points[k].m_x)
			{
				k = j;
			}
		}

		if (k != i)
		{
			swap(P->m_points[i], P->m_points[k]);
		}
	}
}

// Y升序对点集排序;
void SortPointsByY(PointsSet *P, int n)
{
	int k = 0;

	for (int i=0; i<n; i++)
	{
		k = i;

		for (int j=i+1; j<n; j++)
		{
			if (P->m_points[j].m_y < P->m_points[k].m_y)
			{
				k = j;
			}
		}

		if (k != i)
		{
			swap(P->m_points[i], P->m_points[k]);
		}
	}
}

int _tmain(int argc, _TCHAR* argv[])
{
	PointsSet  P;
	PointsSet  X;
	PointsSet  Y;

	srand((unsigned int)time(NULL));

	P.Alloc(MAX_SIZE);
	X.Alloc(MAX_SIZE);
	Y.Alloc(MAX_SIZE);

	for (int i=0; i<MAX_SIZE; i++)
	{
		P.m_points[i].m_x = rand() % 1000;
		P.m_points[i].m_y = rand() % 1000;
		P.m_points[i].m_i = i;
		P.m_iLength++;
	}

	clock_t start1 = clock();

	int iBruteIndex1 = 0, iBruteIndex2 = 0;
	double dMinSqDistBrute = DBL_MAX;
	dMinSqDistBrute = ComputePointsSetSqDistBrute(P, iBruteIndex1, iBruteIndex2, dMinSqDistBrute);

	clock_t end1 = clock();

	double x= P.m_iLength;
	double x1= P.m_points[0].m_x;

	for (int i=0; i<P.m_iLength; i++)
	{
		X.m_points[i] = P.m_points[i];
		X.m_iLength++;
		Y.m_points[i] = P.m_points[i];
		Y.m_iLength++;
	}

	SortPointsByX(&X, MAX_SIZE);
	SortPointsByY(&Y, MAX_SIZE);

	clock_t start2 = clock();
	int iQuickIndex1 = 0, iQuickIndex2 = 0;
	double dMinSqDistQuick = DBL_MAX;
	dMinSqDistQuick = ComputeClosetPoint(P, X, Y, iQuickIndex1, iQuickIndex2, dMinSqDistQuick);
	clock_t end2 = clock();

	if (fabs(dMinSqDistQuick - dMinSqDistBrute) > 1e-8)
	{
		cout<<"运算结果不对 ! "<<dMinSqDistQuick<<"不等于"<<dMinSqDistBrute<<endl;
	}
	else
	{
		cout<<"运算结果正确 ! "<<endl<<"最近距离平方为："<<dMinSqDistQuick<<endl;
	}

	cout<<"蛮力法耗费时间为："<<end1 - start1<<endl;
	cout<<"技巧法耗费时间为："<<end2 - start2<<endl;

	return 0;
}