1. 程式人生 > >計算幾何 平面最近點對 nlogn分治演算法 求平面中距離最近的兩點

計算幾何 平面最近點對 nlogn分治演算法 求平面中距離最近的兩點

本文全文原創 轉載請註明出處

http://blog.csdn.net/lytning/article/details/25370169

平面最近點對,即平面中距離最近的兩點

分治演算法:

int SOLVE(int left,int right)//求解點集中區間[left,right]中的最近點對

{

    double ans; //answer

    0)    呼叫前的預處理:對所有點排序,以x為第一關鍵詞y為第二關鍵字 , 從小到大;


    1)    將所有點按x座標分成左右兩部分;

    /*      分析當前集合[left,right]中的最近點對,有兩種可能:

          1. 當前集合中的最近點對,點對的兩點同屬於集合[left,mid]同屬於集合[mid,right]

             則ans = min(集合1中所有點的最近距離, 集合2中所有點的最近距離)

          2. 當前集合最近點對中的兩點分屬於不同集合:[left,mid][mid,right]

              則需要對兩個集合進行合併,找出是否存在p∈[left,mid],q∈[mid,right],使得distance(p,q)小於當前ans(即步驟1中求得的ans);

    */

    2)    Mid = (left+right)/2;

          ans = min( SOLVE(left,mid), SOLVE(mid,right) );

          即:遞迴求解左右兩部分中的最近距離,並取最小值;

          //此步驟實現上文分析中的第一種情況

    /*      


          再次進行分析

          我們將集合[left,right]用x = mid這條直線分割成兩部分

          則如果畫出直線l1:x=mid-ans 和 l2:x=mid+ans,顯然如果有p∈[left,mid], q∈[mid,right]且distance(p,q) < ans則p,q一定在直線l1和直線l2之間,否則distance(p,q)必定大於ans。

          於是掃描出在l1和l2之間的點

    */

    3)    建立快取陣列temp[];

          for i = left TO right

          {

               如果 abs(Point[i].x - Point[mid].x) <= ans

               則向temp中加入點Point[i];

           }

    /*

            對於temp中的點,列舉求所有點中距離最近兩點的距離,然後與ans比較即可。

            列舉的時候不必兩兩列舉。觀察下圖中的點p

           不難發現,若有q∈[mid,mid+ans]使得distance(p,q) < ans,則q點的位置一定在圖中畫出的一個2ans×ansd的矩形中。可以證明點集[mid,mid+ans]中的、矩形外的點與p點的距離一定大於ans。
           於是我們可以對temp以y為唯一關鍵字從小到大排序,進行列舉, 更新ans,然後在列舉時判斷:一旦列舉到的點與p點y值之差大於ans,停止列舉。最後就能得到該區間的最近點對。

    */

    4)    sort(temp);

          for i = 0 TO k-1

          {

                for j = i+1 TO k-1

                    如果 temp[j].y - temp[i].y >= ans  break;

                    ans = min( ans, distance(temp[i], temp[j]) );

           }

    5)    return ans;

}


演算法的時間複雜度

        由鴿巢原理,程式碼中第四步的列舉實際上最多隻會列舉6個點,效率極高(一種蒟蒻的證明請看下方的評論)

        本演算法時間複雜度為O(n log n)


程式碼:

//是vijos1012的ac程式碼

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;

const double eps = 1e-8;
const int INF = 0x7fffffff;

int n;

struct Point
{
	double x,y;
	Point(double x=0, double y=0):x(x),y(y) {}
	bool operator < (const Point& p) const
	{
		if(x != p.x)	return x < p.x;
		else	return y < p.y;
	}
}p[200000+5],temp[200000+5];

bool cmpy(Point a, Point b)
{
	return a.y < b.y;
}

double Dis(Point a, Point b)
{
	return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));
}

double Closest_Pair(int left, int right)
{
	double d = INF;
	if(left == right)
		return d;
	if(left +1 == right)
		return Dis(p[left],p[right]);
	int mid = (left+right)>>1;
	double d1 = Closest_Pair(left,mid);
	double d2 = Closest_Pair(mid,right);
	d = min(d1,d2);
	int k = 0;
	for(int i = left; i <= right; i++)
	{
		if(fabs(p[mid].x - p[i].x) <= d)
			temp[k++] = p[i];
	}
	sort(temp,temp+k,cmpy);
	for(int i = 0; i < k; i++)
	{
		for(int j = i+1; j < k && temp[j].y - temp[i].y < d; j++)
		{
			double d3 = Dis(temp[i],temp[j]);
			d = min(d,d3);
		}
	}
	return d;
}


int main()
{
	cin>>n;
	for(int i=0; i<n; i++)
	{
		double a,b;
		scanf("%lf%lf",&a,&b);
		p[i] = Point(a,b);
	}
	sort(p,p+n);
	printf("%.3f",Closest_Pair(0,n-1));
}

練習題:

vijos1012  https://vijos.org/p/1012

    裸題,平面最近點對(雖然本題暴力也能過=。=)

poj 3714   http://poj.org/problem?id=3714

    大意:給出平面中的兩類點,求平面最近點對,要求該點對中的兩點不是同類點