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

2019年5月6日 128次阅读 来源: lu1012123053

原创文章，转载请注明出处~~ 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个点进行比较呢？对这个问题，《编程之美》上没有说清楚，这里再详细的论证下；

如下图所示：

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

        在上图中，我们有一系列的点，按照上述的算法，我们利用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*log(n))--寻找最近点对》

根据上面的算思想，我们利用线性时间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;
}

 

运行结果如下：

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

 代码描述：

     1）对点集S的点x坐标和y坐标进行升序排序，获得点集Sx和Sy

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

     3）Divide_conquer（Sx，Sy，δ）  //分治法

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

                  return δ;

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

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

                   return δ;

             else    //如果Sx中多于2个点，则将Sx,Sy分治，以中心点画线，将Sx分为左右两部分SxL和SxR，Sy分为SyL和SyR

                   j1=1,j2=1,k1=1,k2=1;

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

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

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

                   {

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

                                SxL[k1] = Sx[i]   k1++;

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

                                SxR.count[k2] = Sx[i]   k2++;

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

                                SyL[j1] = Sx[i]   j1++;

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

                                SyR[j2] = Sx[i]   j2++;

                   }

              δL = Divide_conquer（SxL，SyL，δ）;    //获取Sx中的的最小点位距离δL

              δR = Divide_conquer（SxR，SyR，δ）;   //获取Sy中的的最小点位距离δR

              δ= min (δL,δR);

              δ=merge(SyL，SyR，δ);    //获Sx中Sy交界处的最小点位距离，并综合 δL和δR 求出点集的最小点位距离δ

      return δ;

 

      函数merge(SyL，SyR，δ)

      merge(SyL，SyR，δ)

      {

          i1=1,i2=1;

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

          {

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

                  then S’yL[i1]= SyL[i], i1++,

           }

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

          {

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

              then S’yR[i2]= SyR[i], i2++,

          }

 

          t=1;

          for(i=1,i<S’yL.count,i++)

           {     

                while(S’yR[t].y< S’yL[t].y and t < SyR.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(n2)要高效。

 

    原文作者：lu1012123053
    原文地址: https://blog.csdn.net/lu1012123053/article/details/9825175
    本文转自网络文章，转载此文章仅为分享知识，如有侵权，请联系博主进行删除。
点赞