原创文章，转载请注明出处~~ http://www.cnblogs.com/justinzhang/
问题描述：给定平面上N个点的坐标，找出距离最近的两个点。
这是编程之美2.11的一道题目，从昨天到现在就一直在设法解决它；如果用常规的解法，只需要将N个点两两计算距离，然后找出最小距离的两个点就可以了；但是这种解法的算法复杂度为O(N^2); 为了降低算法的复杂度，我们需要有更好的方法。这里我们找到的解法是分治法。
设点集为S，|S|=N，S的横坐标集合为Sx，纵坐标集合为Sy;
算法的步骤如下：
1.找出Sx的中位数:median_Sx;用median_Sx对点集S进行划分，左边的为S1，右边的为S2;
2.分别求出S1和S2中的最近点对，设S1和S2中最近点对的距离分别为:delta(S1), delta(S2);
3.求出S1和S2最近点对距离的较小值：Delta = min{delta(S1), delta(S2)};
4.找出S2中y值前6大的点，设为S2’.(此处不用排序，用O(n)时间找出第i大的点，因为排序的时间复杂度为O(n*logn), 网上绝大部分的算法是用排序，显然用排序的方法来找出前6大的点效率低下）;
5.对于S1中的每一个点，与S2’中的每一个点计算距离delta’, 如果delta’ < Delta，令Delta=delta’;
现在我们分析一下上述算法的时间复杂度，由于我们采用中位数进行划分，每一次划分都能确保点集被分为大小相当的两个部分，所以：
T(n)= 2*T(n/2) + 合并复杂度。 由于合并的时候，是用S1中的n/2个点和S2中最多6个点进行比较，比较的次数为n/2 * 6 = 3n. 所以
T(n)= 2*T(n/2)+O(n). 由《算法导论》中文第二版44页的主定理，可知T(n) = O(n*log(n));
算法第4步中，为什么只用选取S2中最多前6大的6个点进行比较呢？对这个问题，《编程之美》上没有说清楚，这里再详细的论证下；
如下图所示：
在上图中，我们有一系列的点，按照上述的算法，我们利用median(中位数)将点分为了S1和S2两个部分，安装算法的第2、3步我们可以求得δ，在图中median线的左边和右边，距离其δ的距离，分别作两条直线，和直线median一起，我们得到了两个区域，分别为a1、a2，如图所示。由于δ是S1和S2中距离较小者，由此可知，a1和a2区域内不可能再存在距离<δ的两个点对。那么，最近点对要么是S1，S2各自内部的某两个点，或者是一个点在S1中、一个点在S2中。在求得δ后，我们只需要查看是否存在两个点，他们满足条件：一个在a1中，一个在a2中；并且它们之间的距离d<δ; 仔细观察上图，我们来证明区域a2中最多只可能存在6个点；
根据δ定义，我们知道，a2中的任意两点的距离都必须大于δ，否则与δ的定义相矛盾。如下图所示，我们将δ*2δ的矩形区域，分成边长为1/2δ *2/3δ的六个小矩形，如图中的红色编号1..6所示。我们假设该区域内的点大于6个，那么根据鸽笼原理，那么必然有一个小矩形中存在至少两个点，这两个点的距离必然小于边长为1/2δ *2/3δ矩形的对角线，即D^2<（1/2δ)^2 + (2/3δ)^2=(5/6δ)^2<δ^2;所以a2矩形区域当中的点数不可能大于6.最多就6个点，这种情况是如上图中红色的A,B…EF所示。
根据上面的算思想，我们利用线性时间O(n)查找中位数的算法，来找到中位数，并且对集合进行划分，这样就确保了每次都是T(n)=T(n/2)+O(n).关于中位数和顺序统计学可以参考：http://www.cnblogs.com/justinzhang/archive/2012/04/24/2469009.html 。
本问题的代码如下所示（有点长），有错误的地方还望大家批评指正~~
/*
Author:JustinZhang
Time:2012年4月24日11:29:25
End: 2012年4月25日16:56:35
Email:uestczhangchao@gmail.com
Desc:平面上有N个点，找出N个点中距离最近的两个点；如果用穷举法，那么算法的复杂度将达到O(n^2)；
本算法采用分治法，可以将复杂度降低到O(n*log(n));
*/
#include <iostream>
#include <vector>
#include <iomanip>
#include <algorithm>
#include <cmath>
using namespace std;
//delta=min{min(dest(S1), min(dest(S2}}.
double delta = INT_MAX;
void swap(int &x, int &y)
{
int tmp = x;
x = y;
y = tmp;
}
//插入排序算法
void insert_sort(vector<int>&A, int p, int r)
{
int i,j;
int key = 0;
for (i=p+1; i<=r; i++)
{
key = A[i];
j = i - 1;
while (j>=p && A[j]>key)
{
A[j+1] = A[j];//将A[j]往前移动
j--;
}
A[j+1] = key;
}
}
//将整型容器划分为两个部分
int partion(vector<int> &A,int p, int r)
{
int count = p - 1;
int key = A[r];
for (int i=p; i<=r-1; i++)
{
if (A[i] < key)
{
count++;
swap(A[count],A[i]);
}
}
swap(A[count+1],A[r]);
return count+1;
}
//根据算法导论上的思想，取得中位数的中位数
int get_median(vector<int>&data, int p, int r)
{
int i=0, j=0;
int n = r - p + 1; //获得原始数组数据个数
int remains = n%5;
int int_count = n - remains;
int int_group = int_count/5;
int group_count = (n+4)/5; //计算出五个元素一组的组数；
//int *group_median = new int[group_count];
vector<int> group_median(group_count);
if (p==r)
{
return data[p];
}
//以下代码求出五个元素为一组的中位数
if (0 == remains) //如果元素的个数正好可以分为以5个元素为单位的整数组
{
for (i=p; i<=r-4; i=i+5)
{
insert_sort(data, i, i+4);
group_median[(i-p)/5] = data[p+2];
}
}
else
{
for (i=p; i<=(r-remains-4);i=i+5)
{
insert_sort(data, i, i+4);
group_median[(i-p)/5] = data[i+2];
}
//处理不够5个元素的组，改组开始的序号为r-remains+1，结束序号为r
insert_sort(data,r-remains+1,r);
group_median[group_count-1] = data[r-remains+1+remains/2];
}
if (group_count==1)
{
return group_median[0];
}
else
{
return get_median(group_median,0,group_count-1);
}
return 0;
}
/*就是因为get_position函数写错了，导致很久都没有能够发现错误，要仔细啊~~*/
int get_position(vector<int> &A, int p, int r, int key)
{
for (int i=p; i<=r; i++)
{
if (A[i]==key)
{
return i;
}
}
return -1;
}
//该函数是为了找到数组A中，第seq小的数
int select(vector<int> &A,int p, int r, int seq)
{
//获得数组A的中位数的中位数，将其作为划分数组的支点
int pivot_key = get_median(A, p, r);
int pos = get_position(A,p,r,pivot_key);
swap(A[pos],A[r]);
int i = partion(A, p, r);
int x = i - p + 1;
if (seq == x)
{
return A[i];
}
else if (seq < x)
{
select(A, p, i-1, seq);
}
else
{
select(A, i+1, r, seq-x);
}
}
class Point
{
public:
int x;
int y;
Point(int x=0, int y=0)
{
this->x = x;
this->y = y;
}
Point& operator=(const Point &p)
{
this->x = p.x;
this->y = p.y;
return *this;
}
Point(const Point &pp)
{
this->x = pp.x;
this->y = pp.y;
}
friend ostream &operator<<(ostream &os, const Point & p)
{
os<<"("<<p.x<<"," << p.y << ")";
return os;
}
};
//将两个点交换
void swap_point(Point &p1, Point &p2)
{
Point tmp = p1;
p1 = p2;
p2 = tmp;
}
//根据点容器的最后一个点，将点集合划分为两个部分
int partion_Point_Set(vector<Point> &p,int l, int r)
{
int count = -1;
int key = p[r].x;
for (unsigned i=0; i<=p.size()-2; i++)
{
if (p[i].x<key)
{
count++;
swap_point(p[i],p[count]);
}
}
swap_point(p[count+1],p[r]);
return count+1;
}
//获得两点间的距离
double Distance(const Point &p1, const Point &p2)
{
int x = (p1.x - p2.x);
int y= (p1.y - p2.y);
double distance = sqrt((double)(x*x+y*y));
return distance;
}
//找到两个数的最小值
double min(double x, double y)
{
if (x>y)
{
return y;
}
else
{
return x;
}
}
//按照点的x坐标来检索点在容器中的位置
int get_point_position_x(const vector<Point> &p, int l, int r, int x_key)
{
for (int i=l; i<=r; i++)
{
if (p[i].x == x_key)
{
return i;
}
}
return -1;
}
//按照点的y坐标来检索点在容器中的位置
int get_point_position_y(const vector<Point> &p, int l, int r, int y_key)
{
for (int i=l; i<=r; i++)
{
if (p[i].y == y_key)
{
return i;
}
}
return -1;
}
//找到p中y坐标第order大的点
Point select_point(vector<Point> &p,int l, int r, int order)
{
vector<int> point_y;
for (int i=l; i<=r; i++)
{
point_y.push_back(p[i].y);
}
int key_y = select(point_y,0,point_y.size()-1,order);
int position_y = get_point_position_y(p, 0, p.size()-1,key_y);
return p[position_y];
}
double get_minimum_distance(vector<Point> &p,int l, int r,vector<Point> &result)
{
int i=0,j=0;
if ((r-l+1)==2)//如果点数为两个
{
double ret = Distance(p[l],p[r]);
if (ret < delta)
{
result[0]=p[l];
result[1]=p[r];
}
return ret;
}
else if ((r-l+1)==3) //如果点数为3个
{
double tmp_dis1 = Distance(p[l],p[l+1]);
double tmp_dis2 = Distance(p[l],p[l+2]);
double tmp_dis3 = Distance(p[l+1],p[l+2]);
double ret = min(min(tmp_dis1,tmp_dis2),tmp_dis3);
if (ret <delta )
{
if (tmp_dis1 == ret)
{
result[0] = p[l];
result[1] = p[l+1];
}
else if (tmp_dis2==ret)
{
result[0]=p[l];
result[1]=p[l+2];
}
else
{
result[0]= p[l+1];
result[1]= p[l+2];
}
}
return ret;
}
else //大于三个点的情况
{
//求出点集p的x坐标和y坐标
vector<int> Pointx;
vector<int> Pointy;
for (i=l; i<=r; i++)
{
Pointx.push_back(p[i].x);
Pointy.push_back(p[i].y);
}
//找到点x坐标的中位数
int x_median = select(Pointx,0,Pointx.size()-1,(Pointx.size()+1)/2);
int point_median_pos = get_point_position_x(p,l,r,x_median);
swap_point(p[point_median_pos],p[r]);
//利用找到的中位数对点集进行划分
int point_pivot_index = partion_Point_Set(p,l,r);
//利用x坐标作为中位数，将点集合划分好后，递归的求中位数左边点集S1和右边点集合S2距离最小值；
//考虑如何将距离最小的两个点保存下来
double min_s1 = get_minimum_distance(p,l,point_pivot_index,result);
double min_s2 = get_minimum_distance(p,point_pivot_index+1,r,result);
if (min_s1>min_s2)
{
delta = min_s2;
}
else
{
//result[0] = tmp_result1;
//result[1] = tmp_result2;
delta = min_s1;
}
//现在已经递归的求到了S1和S2中点集合最短距离，下面开始求S1和S2之间点之间的距离
//找出点集合S2中，y坐标前六大的点,如果|S2|<6，则只需找出|S2|个点
int size_s2 = (r-point_pivot_index<6)?r-point_pivot_index:6;
vector<Point> S2;
Point tmp;
for (i=1;i<=size_s2;i++)
{
tmp = select_point(p,point_pivot_index+1,r,i);
S2.push_back(tmp);
}
for (i=l; i<=point_pivot_index; i++)
{
for (j=0; j<size_s2; j++)
{
double d = Distance(p[i],S2[j]);
if (d < delta)
{
result[0]= p[i];
result[1]= S2[j];
delta = d;
}
}
}
return delta;
}
}
void getinput(vector<Point> &p)
{
int i;
Point pp;
int Pointnumber = 0;
cout << "please input Point numbers" <<endl;
cin >> Pointnumber;
for (i=0; i<Pointnumber; i++)
{
cin >> pp.x;
cin >> pp.y;
p.push_back(pp);
}
}
int main()
{
vector<Point> result(2);
vector<Point> input;
getinput(input);
double minimum_distance = get_minimum_distance(input,0,input.size()-1,result);
cout << "The nearest point is: "<<result[0] << " and " << result[1] << endl;
cout << "their distance is : "<<minimum_distance << endl;
return 0;
}
运行结果如下：
题目：平面中有若干个点，寻找距离最近的两个点。
分析：
方法1：两两点比较，寻找最近的两个点对，复杂度O(N^2)，优点代码简单不容易出错
方法2：观察两两比较的方法，发现有很多无用的比较，对于每一个点只要计算到它最近的点的距离就可以了，枚举所有的点，最后得出距离最近的一对点，但对于一个给定的点，如何找到距离它最近的点呢？可以使用一些启发式规则，减少比较的次数，例如：对于(x,y)取x轴投影上最近的点p1（x1,y1）和取y轴上投影最近的点p2(x2,y2)，以（x，y）为一个定点，取x1，x2中较大者x'，y1，y2中较大者y',构成（x‘，y’）为矩形的另一个顶点，计算矩形内所有的点与(x,y)的距离，从而减少了比较的次数。再以(x,y)为中心，判断4个项限内的所有点，最后得出距离（x，y）最小的点。复杂度是O(N*K)，K为要判断的点的平均个数
方法3：采用分治法，即从x轴将数据不断分成两部分，两部分最近的两点取其小，在进行两部分中间比较，通过两部分的最小距离，缩小了中间带状的点数量，详见编程之美的分析。
按照方法3的程序如下：
[cpp] view plaincopy
- #include <stdio.h>
- #include <algorithm>
- #include <vector>
- #include <math.h>
- class Point {
- public:
- Point(int x, int y) : x_(x), y_(y) {}
- Point() : x_(0), y_(0) {}
- static bool OrderByX(const Point& left, const Point& right) {
- 10. return left.x_ < right.x_;
- 11. }
- 12. static bool OrderByY(const Point& left, const Point& right) {
- 13. return left.y_ < right.y_;
- 14. }
- 15. int x_;
- 16. int y_;
17. };
18. float Distance(const Point& left, const Point& right) {
- 19. return sqrt(pow(left.x_ - right.x_, 2) + pow(left.y_ - right.y_, 2));
20. }
21. int NearestPoints(const std::vector<Point>& points, int start, int end, Point* point1, Point* point2) {
- 22. if (end > start) {
- 23. int middle = (start + end) / 2;
- 24. int left_min_distance = NearestPoints(points, start, middle, point1, point2);
- 25. int right_min_distance = NearestPoints(points, middle + 1, end, point1, point2);
- 26. int min_distance = left_min_distance > right_min_distance ? right_min_distance : left_min_distance;
- 27. std::vector<Point> left_part_points;
- 28. for (int i = start; i <= middle; ++i) {
- 29. if (points[middle].x_ - points[i].x_ <= min_distance) {
- 30. left_part_points.push_back(points[i]);
- 31. }
- 32. }
- 33. sort(left_part_points.begin(), left_part_points.end(), Point::OrderByY);
- 34. std::vector<Point> right_part_points;
- 35. for (int i = middle + 1; i <= end; ++i) {
- 36. if (points[i].x_ - points[middle].x_ <= min_distance) {
- 37. right_part_points.push_back(points[i]);
- 38. }
- 39. }
- 40. sort(right_part_points.begin(), right_part_points.end(), Point::OrderByY);
- 41. int distance_y = 0;
- 42. int point_distance = 0;
- 43. for(int i = 0; i < left_part_points.size(); ++i) {
- 44. for(int j = 0; j < right_part_points.size(); ++j) {
- 45. distance_y = left_part_points[i].y_ > right_part_points[j].y_ ? left_part_points[i].y_ - right_part_points[j].y_ :
- 46. right_part_points[j].y_ - left_part_points[i].y_;
- 47. if (distance_y <= min_distance) {
- 48. point_distance = Distance(left_part_points[i], right_part_points[j]);
- 49. if (point_distance < min_distance) {
- 50. min_distance = point_distance;
- 51. *point1 = left_part_points[i];
- 52. *point2 = right_part_points[j];
- 53. }
- 54. }
- 55. }
- 56. }
- 57. return min_distance;
- 58. } else {
- 59. return 0x7FFFFFFF;
- 60. }
61. }
- 62.
- 63.
64. int main(int argc, char** argv) {
- 65. std::vector<Point> points;
- 66. points.push_back(Point(2,3));
- 67. points.push_back(Point(1,4));
- 68. points.push_back(Point(3,0));
- 69. points.push_back(Point(5,0));
- 70. points.push_back(Point(5,1));
- 71. sort(points.begin(), points.end(), Point::OrderByX);
- 72. Point point1;
- 73. Point point2;
- 74. NearestPoints(points, 0, points.size() - 1, &point1, &point2);
- 75. printf("Point1: (%d, %d) <--> Point2: (%d, %d) ", point1.x_, point1.y_, point2.x_, point2.y_);
76. }
其中，寻找矩形带内的点可以通过二分搜索来寻找，为了简单，用了循环来寻找。