寻找最近点（快速算法)）

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

应用：交通控制中寻找两个最近的交通工具。
// 完整源码：http://download.csdn.net/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;}