Problem Description:

Suppose there are N points on the plane and find the distance between the nearest two points (Euclidean distance).

Algorithm analysis

  • A typical divide-and-conquer solution.
  1. Divide all points into left and right parts, and then work out the shortest distance minD between the left and right parts.
  2. Find the shortest distance of the middle intersection, so as to get the minimum distance of the whole section.
  3. The question is how to find the shortest distance between the intersecting parts.
  • If you use the simple brute force method, which is to search the shortest distance between all the points in the left and right interval, you can easily get O(n^2).

An optimization method

  • We find that if the difference between the x coordinate of any point P on the left and right part and the x coordinate of the middle point mid is greater than minD, there is no need to compare, so we get a smaller cross strip.
  • If strip is found every time and array Q is obtained by sorting y, then minD is obtained by increasing y. Obviously, if the difference of y coordinates between the two points before and after is greater than minD, there is no need to compare the last points in Q. But sorting is nlogn which makes the whole algorithm nlognlogn. That’s algorithm 1 in the code below.
  • And our goal is to optimize, divide by X and Y each time, so that the smallest interval that can be solved is ordered by Y, and then merge, so that the sorting is preprocessed, so that the time to deal with the intersections is reduced to order N, so that the overall time complexity can be reduced to ORDER O(nlogn). Which is algorithm 2 down here.
    • Approach is to use the merge sort of thought, can get solution of the minimum interval around (that is, a difference of 1, 2), we first to y sorting (because here most will only have three elements, so the efficiency is very high, in fact, that is, to divide y), and then in the merge stage, because of the paragraphs of y is the orderly, we only need to merge.

C + + code

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

#define N 300000const 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;returntmp > 0 ? 1 : fabs(tmp) < EPS ? 0:1; } int cmpY(const void *lh, const void *rh) { double tmp = ((Point*)lh)->y - ((Point*)rh)->y;returntmp > 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); } double disX(Point p1,Point p2){return(p1.x-p2.x)*(p1.x-p2.x); } double disY(Point p1,Point p2){return(p1.y-p2.y)*(p1.y-p2.y); Double minDis(Point P[],Point T[],Point Q[],int left,int right){// Double minDis(Point P[],Point T[],Point Q[],int left,int right)if(left== right) return INF;
    if(left + 1 == right)
        returndis(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; // The mark of the point in the strip areafor(i=left,k=0; i<=right; i++){if(disX(P[i],P[mid]) < minD){ Q[k++] = P[i]; 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; }}returnminD; } void merge(Point T[],Point Q[], int left,int mid,int right){int index =left,l_start =left; Int r_start = mid + 1; // For the position starting from the right, we need to increment by 1while(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++]; Memcpy (&q [left], &t [left], (right-left+1)*sizeof(Q[0])); } // algorithm 2 //P sort by X //Q is set to any value at first, we set it to P, step by step merge to sort by X //T is cache array, 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 sort */returndis(P[left],P[right]); } // Three pointsif(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 sort */returnmin(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); Merge (T,Q,left,mid,right); merge(T,Q,left,mid,right); int i,k; // The mark of the point in the strip areafor(i=left,k=0; i<=right; I++){// the X coordinate difference exceeds midD, obviously should be removedif(disX(Q[i],P[mid]) < minD){ T[k++] = Q[i]; }} // Because we already sort Y when we merge, there is no need to compare Y if the difference between Y coordinates is greater than T after minDfor(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; }}returnminD; } double fun(Point P[],int n,Funtype fp){ double minD; Point *Q = (Point *)malloc(n*sizeof(Point)); Point *T = (Point *)malloc(n*sizeof(Point)); Qsort (P, n, sizeof(P[0]),cmpX); */ 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; } double randFloat(int min,int Max){double m = (double)(rand()%101)/101; double n = (double)(rand()%(max - min) + min); // Get the range [min,max-1]return m+n;
}

Point randPoint(){ Point p; double x = randFloat(0,N); Double y = randFloat(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 storage array of ordered by X / / modify N = 12 / / Point P [12] = {{1.1, 3.0}, {1.9, 3}, {5.7, 7}, {5.25, 7.3}, {39 12:3261}, {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);
}
Copy the code

Algorithm complexity analysis

We took 300,000 points here (the array with more points in C++ has crossed the line), ran it 10 times and found

  • Time Used of algorithm 1 :2490154
  • The Time used of Algorithm 2 is 2017634, which has saved 1/5 Time compared with algorithm 1.
  • Basically, algorithm 2 O(nlogn) is very efficient.