分治法解决平面上N点最近2点距离———算法应该OK~~~

问题描述:给定平面上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个点进行比较呢?对这个问题,《编程之美》上没有说清楚,这里再详细的论证下;

如下图所示:

image

在上图中,我们有一系列的点,按照上述的算法,我们利用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所示。

image

根据上面的算思想,我们利用线性时间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:[email protected]
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;
}

 

原文链接:加载失败,请重新获取