优化的最近点问题排序法和预处理法

tcfellow
1,386 阅读6分钟

问题描述：

设平面上有N个点，求最近的两个点之间的距离（欧氏距离）。

算法分析

  • 典型的用分治思想来解决的问题。
  1. 把所有点划分成左右两部分，然后求出左右两部分的最短距离minD。
  2. 求中间交叉部分的最短距离，从而得到整段的最小距离。
  3. 问题在于怎么求交叉部分的最短距离。
  • 如果采用简单暴力的方法，也就是搜索所有左右区间的点的最短距离，那么容易得出复杂度为O(n^2)。

优化方法

  • 我们发现如果左右部分任意点P的x坐标与中间点mid的x坐标差大于minD，那就没有比较的必要了，这样我们得到了一个更小的交叉带strip。
  • 如果每次找出交叉部分strip，然后对y排序得到数组Q，然后我们按照y递增的方式求minD，很明显如果前后两个点的y的坐标差大于minD，那么Q中后面的就没有比较的必要了。但是排序是nlogn复杂度会使得整个算法的复杂度为nlognlogn。也就是下面代码的算法1。
  • 而我们的目标是进行优化，每次按照X，Y进行划分，使得能求解的最小区间按Y有序，然后再合并，排序就被预处理了，就使得处理交叉部分情况所花费的时间降到n级，这样整体的时间复杂度就可以降到O(nlogn)。也就是下面的算法2.
    • 做法是利用归并排序的思想，能求得解的最小区间（也就是左右相差1，2），我们先对y排序（由于这里最多只会有3个元素，所以效率很高，实际上也就是对y进行划分），然后在归并阶段，由于各段的y已经是有序的了，我们只需要合并即可。

C++代码

#include <iostream>
#include <cmath>
#include <algorithm>
#include <sys/time.h>
using namespace std;

#define N 300000
const double INF=9e300;
const double EPS=0.00000000001;


typedef struct Point{
    double x;
    double y;
}Point;  //点集合

//指针函数，便于统一测试
typedef double (*Funtype)(Point [],Point[],Point[],int ,int);

int cmpX(const void *lh, const void *rh)
{
    double tmp = ((Point*)lh)->x - ((Point*)rh)->x;
    return tmp > 0 ? 1 : fabs(tmp) < EPS ? 0 : -1;
}

int cmpY(const void *lh, const void *rh)
{
    double tmp = ((Point*)lh)->y - ((Point*)rh)->y;
    return tmp > 0 ? 1 : fabs(tmp) < EPS ? 0 : -1;
}

double dis(Point p1,Point p2){
    return (p1.x-p2.x)*(p1.x-p2.x) + (p1.y-p2.y)*(p1.y-p2.y);
}
//x距离
double disX(Point p1,Point p2){
    return (p1.x-p2.x)*(p1.x-p2.x);
}
//y距离
double disY(Point p1,Point p2){
    return (p1.y-p2.y)*(p1.y-p2.y);
}

//算法1
//这里第二个T参数没用，我们是为了统一测试
double minDis(Point P[],Point T[],Point Q[],int left,int right){
    //边界条件，一个点无穷大，两个点直接取就可以了
    if(left== right) return INF;
    if(left + 1 == right)
        return dis(P[left],P[right]);
    int mid = (left + right)>>1;;
 
    double min1 = minDis(P,Q,T,left,mid);
    double min2 = minDis(P,Q,T,mid + 1,right);
    double minD = min(min1,min2);
    int i,k; //带状区域的点的标志
    for(i=left,k=0;i<=right;i++){
        if(disX(P[i],P[mid]) < minD){
            Q[k++] = P[i];    //Q作为缓存数组
        }
    }
    qsort(Q, k, sizeof(Q[0]),cmpY);
    for(int i=0;i<k;i++){
        for(int j= i+1; j< k && disY(Q[j],Q[i]) < minD; j++){
            double temp = dis(Q[i],Q[j]);
            if(temp< minD) minD = temp;
        }
    }  
    return minD;
}

void merge(Point T[],Point Q[], int left,int mid,int right){
    //左右一起合并并缓存到T中
    int index =left,l_start = left;  //新数组下标
    int r_start = mid + 1;  //右边开始的位置，需要先加1
    while(l_start <=mid && r_start <=right){
        if(Q[l_start].y < Q[r_start].y)
            T[index++] = Q[l_start++];   
        else 
            T[index++] = Q[r_start++];
    }
    while(l_start<=mid)
        T[index++] = Q[l_start++];
    while(r_start<=right)
        T[index++] = Q[r_start++];
    //复制过来，保证下次归并的时候Q依然是有序的
    memcpy(&Q[left],&T[left], (right-left+1)*sizeof(Q[0]));  
}

//算法2
//P按照X排序
//Q一开始设为任意值，我们设为了P，逐步归并为按X排序
//T为缓存数组，方便更新Q
double minDisWithoutSort(Point P[],Point Q[],Point T[],int left,int right){
    //边界条件，一个点无穷大，两个点直接取就可以了
    if(left== right) return INF;
    if(right - left == 1){
        qsort(Q+left, 2, sizeof(Q[0]),cmpY); /* Y排序 */ 
        return dis(P[left],P[right]);
    }
    //三个点的情况
    if (right - left == 2) 
    { 
        double d1 = dis(P[left], P[left + 1]);
        double d2 = dis(P[left + 1], P[right]);
        double d3 = dis(P[left], P[right]);
        qsort(Q+right, 3, sizeof(Q[0]),cmpY); /* Y排序 */ 
        return min(min(d1,d2),d3);
    }
    int mid = (left + right)>>1;;
    double min1 = minDisWithoutSort(P,Q,T,left,mid);
    double min2 = minDisWithoutSort(P,Q,T,mid + 1,right);
    double minD = min(min1,min2);
    //合并按照Y排序的序列Q,保证整体有序，这样在最后一个循环里就可以减少很多判断
    merge(T,Q,left,mid,right);
    int i,k; //带状区域的点的标志
    for(i=left,k=0;i<=right;i++){
        //X坐标差超过midD，明显就该舍去
        if(disX(Q[i],P[mid]) < minD){
            T[k++] = Q[i];    
        }
    }
    //因为在merge的时候我们已经对Y排序了，如果Y坐标差超过minD后面的T也没必要比较了
    for(int i=0;i<k;i++){
        for(int j= i+1; j< k && disY(T[j],T[i]) < minD; j++){
            double temp = dis(T[i],T[j]);
            if(temp< minD) minD = temp;
        }
    }    
    return minD;
}



double fun(Point P[],int n,Funtype fp){
    double minD;
    Point *Q =  (Point *)malloc(n*sizeof(Point));  //存储按Y排序的数组
    Point *T =  (Point *)malloc(n*sizeof(Point)); //用作缓存数组
    qsort(P, n, sizeof(P[0]),cmpX ); /* X排序 */ 
    memcpy(Q,P, n*sizeof(P[0]));
    minD = fp(P,Q,T,0,n-1);
    free(Q);
    free(T);
    return sqrt(minD);
}


void test(Point P[],int n,Funtype fp){
    struct timeval start, end;
    unsigned long diff;
    double minD;
    gettimeofday(&start, NULL);
    for (int i = 0; i < 10; ++i)
    {
        minD = fun(P,n,fp);
    }
    gettimeofday(&end, NULL);
    diff = 1000000*(end.tv_sec - start.tv_sec)+ end.tv_usec - start.tv_usec;
    cout<<"Result:"<<minD<<endl;
    cout<<"Time used:"<<diff<<endl;
}

//随机生成[min，max]区间的小数
double randFloat(int min,int max){
    double m = (double)(rand()%101)/101;
    double n = (double)(rand()%(max - min) + min); //得到的值域为[min,max-1]
    return m+n;
}

Point randPoint(){
    Point p;
    double x = randFloat(0,N);  //生成0-N之间的小数
    double y = randFloat(0,N);
    return  p = {x,y};
}

int main(){ 
    Point P[N];
    for(int i=0;i<N;i++){
        P[i]=randPoint();
    }
    //P存储按X排序的数组
    //修改N=12
    // Point P[12] = {{1.1, 3.0}, {1.9, 3}, {5.7,7}, {5.25,7.3}, {39,63}, {31,6}, {1.36, 113.0}, {11.9, 3}, {1.7,17}, {9.25,17.3}, {9,6.3}, {1.0,6.0}};

    test(P,N,minDis);
    test(P,N,minDisWithoutSort);
}

算法复杂度分析

我们在这里采用了300000个点（采用更多的点C++的数组已经越界了），运行了10次，发现

  • 算法1的Time used:2490154
  • 算法2的Time used:2017634，相对算法1已经节省了1/5的时间。
  • 基本上算法2的O(nlogn)的效率还是很高的。
avatar
文章
阅读
粉丝