寻找最近点对

原创文章，转载请注明出处~~ 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 plain copy

#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. }

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. }

其中，寻找矩形带内的点可以通过二分搜索来寻找，为了简单，用了循环来寻找。