编程之美--O(n*log(n))--寻找最近点对

原创文章，转载请注明出处~~ 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:[e-mail]

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）对点集S的点x坐标和y坐标进行升序排序，获得点集S_x和S_y

     2）令δ=∞;   //δ为最小点位距离

     3）Divide_conquer（S_x，S_y，δ）  //分治法

             if （S_x.count=1） then δ=∞;    //如果S_x中只有一个点，则δ=∞

                  return δ;

             else if（S_x.count=2 and d（S_x.[0],S_x.[1]）<δ） //如果S_x中只有2个点，则δ为两点之间距离

                   δ=d（S_x.[0],）S_x.[1]）; 

                   return δ;

             else    //如果S_x中多于2个点，则将S_x,S_y分治，以中心点画线，将S_x分为左右两部分S_xL和S_xR，S_y分为S_yL和S_yR

                   j₁=1,j₂=1,k₁=1,k₂=1;

                   mid = S_x.count/2;   //mid为S_x中的中间点点号

                   L = S_x.[mid].x;      //L为S_x中的中间点x坐标

                   for(i=1,i<S_x.count,i++)

                   {

                         if(i<=mid)     //将S_x中间线以左地方的点存入到S_xL，新数组保持原来的升序性质

                                S_xL[k₁] = S_x[i]   k₁++;

                         else   //将S_x中间线以右的地方的点存入到S_xR，新数组保持原来的升序性质

                                S_xR.count[k₂] = S_x[i]   k₂++;

                         if(S_y[i].x <L)   //将S_y中间线以左地方的点存入到S_yL，新数组保持原来的升序性质

                                S_yL[j₁] = S_x[i]   j₁++;

                         else   //将S_y中间线以右地方的点存入到S_yR，新数组保持原来的升序性质

                                S_yR[j₂] = S_x[i]   j₂++;

                   }

              δ_L = Divide_conquer（S_xL，S_yL，δ）;    //获取S_x中的的最小点位距离δ_L

              δ_R = Divide_conquer（S_xR，S_yR，δ）;   //获取S_y中的的最小点位距离δ_R

              δ= min (δ_L,δ_R);

              δ=merge(S_yL，S_yR，δ);    //获S_x中S_y交界处的最小点位距离，并综合 δ_L和δ_R 求出点集的最小点位距离δ

      return δ;

 

      函数merge(S_yL，S_yR，δ)

      merge(S_yL，S_yR，δ)

      {

          i₁=1,i₂=1;

          for(i=1,i<S_yL.count,i++)   //获取S_yL中在左边虚框(距离小于δ)内的点，存入到S'_yL中，新数组保持原来的升序性质

          {

              if(S_yL[i].x>L-δ)

                  then S'_yL[i₁]= S_yL[i], i₁++,

           }

          for(i=1,i<S_yR.count,i++)  //获取S_yR中在右边虚框(距离小于δ)内的点，存入到S'_yR中，新数组保持原来的升序性质

          {

              if(S_yR[i].x<L+δ)

              then S'_yR[i₂]= S_yR[i], i₂++,

          }

 

          t=1;

          for(i=1,i<S'_yL.count,i++)

           {     

                while(S'_yR[t].y< S'_yL[t].y and t < S_yR.count)  //获取点集S'_yR内距离点S'_yL[t]y坐标最接近的点号

                { t++; }

                for( j= max(1,t-3), j<=min(t+3,S'_yR.count),j++)   //计算S'_yR中的点与S'_yL[t]y坐标相邻的六个点的距离

                {

                      if(d(S'_yL[i],S'_yL[j])<δ)    //如果前两点之间距离小于δ

                      then δ = d(S'_yL[i],S'_yL[j]);   //则最小点位距离δ为当前两点之间距离

                }

          return δ

      }

 

 3）算法时间复杂度：

      首先对点集S的点x坐标和y坐标进行升序排序，需要循环2nlogn次，复杂度为O(2nlogn)

      接下来在分治过程中，对于每个S'_yL中的点，都需要与S'_yR中的6个点进行比较

            O(n)= 2O(n/2) + (n/2)*6  （一个点集划分为左右两个点集，时间复杂度为左右两个点集加上中间区域运算之和）

      其解为O(n)< O(3nlogn)

     因此总的时间复杂度为O(3nlogn)，比蛮力法的O(n²)要高效。