算法思想
本实验使用了分治法的算法思想，所谓分治，就是把一个复杂的问题分成两个或更多的相同或相似的子问题，再把子问题分成更小的子问题……直到最后子问题可以简单的直接求解，原问题的解即子问题的解的合并。
因此使用分治法解决问题一般有以下三个步骤：
分解：将问题划分为一些子问题，子问题的形式与原问题一样，只是规模更小。
解决：递归的求解子问题，当子问题的规模足够小时，停止递归直接求解。
合并：将子问题的解组合成原问题的解。
逻辑上看，分治法是一个自顶向下求解子问题然后自底向上合并的算法，整体算法结构是一个树形结构。
问题描述
给定n个二维平面上的点，求一组点对，它们之间的几何距离最短。
算法描述
对于平面上给定分布的n个点，若使用暴力穷举所有情况，则有 Cn2=2n(n+1) 种情况，对应的时间复杂度是O(n2)。但我们可以应用分治法给出一个时间复杂度为O(nlog n)的算法。该算法在 1975 年由 Franco P. Preparata 提出，Preparata 和 Michael Ian Shamos 证明了该算法在决策树模型下是最优的。
该算法将整个平面按点的x坐标平均分成两部分，那么原本问题中的n个节点中的最短距离就被划分成了3部分：左边 n/2 个节点中的最短距离，右边 n/2 个节点中的最短距离，节点对跨越左右的最短距离。其中前面两个问题是规模为原问题一半的类型相同的子问题，显然的，原问题的解由上面三个问题解的最小者构成。也即如果我们求出了上面三个问题的解d_l、d_r、d_c，那么原问题的解就是 min(d_l,d_r,d_c)。至此我们已经完成了分治法的思路部分。接下来问题所在是：我们应当考虑递归停止的边界条件与第三个子问题（跨越左右节点对的最短距离）的解决方法。
递归停止的边界条件：
（1）当分治得到的子图中只有一个节点时，该图中的最近点对距离为无穷大。
（2）当分治得到的子图中只有两个节点时，该图中最近点对距离就是这两个点之间的几何距离。
跨越左右节点对的最短距离求法：
首先我们先用策略缩小搜索范围，可以知道，令 d = min(d_l,d_r)，则跨越左右点对的最小距离若要比d还小，则一定需要在划分线的左右各长度为d的狭带区域内。
同样的道理，如果跨越左右点对的最小距离若要比d还小，则一定需要在下图中长为2d，宽为d的矩形区域内，因此我们把搜索范围限制在了这样一个矩形范围内。
那么在上述范围内的点的个数是否是有上界的呢？答案是肯定的。我们把矩形范围划分成下图中6个区域，那么一定有：每个区域最多存在一个节点。
原因是若下图中的6个区域中若有一个区域含有2个节点，则其距离必然小于等于(2d)2+(32d)2=65d，然而我们有：两边的最近点对的距离为d，此时我们得到了一个距离为65d的点对，这与已知条件产生矛盾，因此下图6个区域中，每个区域内最多有一个节点。
我们可以知道，最坏情况下，左边狭带内有n/2个节点，而我们需要对左边狭带内所有节点寻找它的上面说的六个节点。因此我们除了子问题外额外的开销是O(n)的代价。由主定理，我们可以得到这是一个时间复杂度为O(nlog n)的算法。
现在的问题是，我们知道对于每个节点，我们最多要找6个节点，检查其距离是否更优，但是我们找到对于每个左边狭带内的节点找到这六个节点是否也花费时间O(n)呢？
并不一定，但我们尝试找到了一种方法：
我们需要对右边狭带内的点按y排序，这个时间复杂度当然是O(nlog n)的，若我们把它加在分治过程中，额外开销将从O(n)增加到O(nlog n)，算法总体的时间复杂度将涨到O（nlog2n），但我们可以将这个排序放在程序最初，由于在程序中我们要按照x将节点平均分成两部分，此时又有对x排序的需求。我们将使用两个不同的数组存储节点，一个按x排序，另一个按y排序，开销都是O(nlog n)，算法的其他部分开销也是O(nlog n)，总体仍然能保证O(nlog n)的时间复杂度。
排序后，我们对于狭带左边内的点开始进行遍历，对于开始的第一个（y最小的节点）找该结点对应的处在矩形范围内的点（最多6个）。然后对于下一个左边狭带内的点，我们不必再考虑那些y坐标要小于上一左边狭带的点的那些点，因为左边狭带内的点距离必然大于等于d，所以左边狭带内的下一点与y坐标要小于左边狭带的上一点的那些点之间的距离必然大于d。因此不必考虑那些点。为了排除掉那些点，我们有必要将左边狭带内的点也加入数组，这样我们每次只需要遍历上一节点的下面至多6个节点即可，至此我们实现了对于任意一个给定的左边节点，我们都能在常数时间内找到它对应的矩形区域内的至多6个节点，也即我们实现了扫描一遍就能找到所有左边狭带内节点的对应矩形区域节点。因此在最后的表达式中，虽然我们采用2层循环的方式，但事实上我们这部分操作只花费了O(n)的时间复杂度。
算法复杂度分析
可以看出，我们把原来的问题分成了两个规模更小的相同性质子问题与一个新的问题，这个问题所需要的额外开销从O(nlog n)被我们降低到了O(n)。因此我们给出该算法的递归方程：
T(n)={O(1)2T(2n)+O(n)n<4n≥4
由主定理，我们可以得到，该算法的时间复杂度为O(nlog n)。
算法实现
#include<iostream>
#include<ctime>
#include<random>
#include<vector>
#include<algorithm>
#include<math.h>
#include<stdlib.h>
using namespace std;
class Point
{
public:
Point(double x, double y);
void printPoint();
double x;
double y;
int order;
bool operator<(Point rhs)const {
return x < rhs.x;
};
Point& operator=(const Point rhs) {
x = rhs.x;
y = rhs.y;
order = rhs.order;
return *this;
}
};
Point::Point(double x, double y)
{
this->x = x;
this->y = y;
}
void Point::printPoint()
{
cout << "(" << x << "," << y << ")" << endl;
}
vector<Point> getRandomPoint(int num)
{
default_random_engine e(time(NULL));
uniform_real_distribution<double> u(-10, 10);
vector<Point> vec;
for (int i = 0; i < num; i++)
{
vec.push_back(Point(u(e),u(e)));
}
sort(vec.begin(), vec.end());
cout << "这些点的坐标为：" << endl;
for (int i = 0; i < num; i++)
{
vec[i].printPoint();
vec[i].order = i;
}
return vec;
}
double getDistance(Point &a, Point &b)
{
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
double getClosestPoint(vector<Point> &xVec, vector<Point> &yVec,int begin,int end,int &lpos,int &rpos)//分治法
{
if (begin>=end)//边界条件
{
return 99999;
}
else if (begin==end-1)//边界条件
{
lpos =xVec[begin].order;
rpos = xVec[end].order;
return getDistance(xVec[begin],xVec[end]);
}
else if (begin==end-2)
{
double minlm = getDistance(xVec[begin], xVec[begin + 1]);
double minlr = getDistance(xVec[begin], xVec[begin + 2]);
double minmr = getDistance(xVec[begin+1], xVec[begin + 2]);
lpos = begin;
rpos = begin+1;
if (minlm>minlr)
{
minlm = minlr;
rpos = begin+2;
}
if (minlm>minmr)
{
minlm = minmr;
lpos = begin + 1;
rpos = begin + 2;
}
return minlm;
}
else
{
int mid = (begin + end) / 2 ,lpos1,lpos2,rpos1,rpos2;
double midx = xVec[mid].x;
//分治子问题
double leftMin = getClosestPoint(xVec, yVec,begin, mid, lpos1, rpos1);
double rightMin = getClosestPoint(xVec, yVec, mid + 1, end, lpos2, rpos2);
double patitionLen;
if (leftMin<rightMin)
{
patitionLen = leftMin;
lpos = lpos1;
rpos = rpos1;
}
else
{
patitionLen = rightMin;
lpos = lpos2;
rpos = rpos2;
}
double minCrossLen = patitionLen;
vector<Point> ctrVec;
for (int i = 0; i < yVec.size(); i++)//将狭带内节点加入数组
{
if (fabs(yVec[i].x-midx)<patitionLen)
{
ctrVec.push_back(yVec[i]);
}
}
for (int i = 0; i < ctrVec.size(); i++)//判断跨越左右的节点（最多6个）是否产生更小的几何距离
{
if (ctrVec[i].x < midx)
{
for (int j = i + 1; j < ctrVec.size() && ctrVec[j].y - ctrVec[i].y < patitionLen; j++)
{
double distance = getDistance(ctrVec[i], ctrVec[j]);
if (distance < minCrossLen)
{
minCrossLen = distance;
lpos = ctrVec[i].order;
rpos = ctrVec[j].order;
}
}
}
}
return minCrossLen;
}
}
int main(void)
{
cout << "请输入平面上点的个数（系统将会随机生成点的坐标）：" << endl;
int num;
cin >> num;
vector<Point> xVec=getRandomPoint(num);
vector<Point> yVec = xVec;
sort(yVec.begin(), yVec.end(), [](Point lhs, Point rhs) { return lhs.y < rhs.y; });
int lpos, rpos;
cout << "这些点间的最小距离为：" << getClosestPoint(xVec, yVec, 0, num - 1, lpos, rpos)<< endl;
cout << "这两个点的序号为：" << lpos << "," << rpos << endl;
system("pause");
}
进一步讨论
我们试讨论3维情况下的最近点对问题，也即空间最近点对问题。
情况如同2维那样，我们拿x做一个垂面，将空间内的点都分成两部分，一部分在该垂面左边，另一部分在垂面右边，划分时要求左右两边的点个数相同。我们把左右两端空间最短距离称作d，则跨越左右空间的最近最距离的点对一定在d×2d×2d的空间内（如图）
对于平面情况下，我们认为每个点最多找到6个点需要判断是否取得更小的距离。对于空间距离，我们采用下面的分割方法：d/2×d/2×2d/3
d×2d×2d的空间将被分割成24块，由于每块内的最大距离为
(2d)2+(2d)2+(32d)2=1817d
，仍然不到d，因此每块内最多只能有一个节点，这意味对于左边的每个节点，最多只需检查24个节点即可，这个时间复杂度仍然是个O（n）的。因此我们给出递归关系式：
T(n)={O(1)2T(2n)+O(n)n<4n≥4
由主定理，最终时间复杂度仍然为O(nlog n)，但该算法最差情况下前面的常数相当大，但同时也应考虑到，平均情况下，点的分布情况是均匀的，我们一般很少会检查全部24个节点，因为稀疏情况是更常见的。
我们给出它的一个伪代码：
参考文献
1 Thomas H.Cormen, Charles E.Leiserson, Ronald L.Rivest, Clifford Stein. Introduction to Algorithms, The MIT Press, 2009
2 胡金初,张晓红. 空间最近点对的计算机算法研究. 计算机科学, Vol35, 2008