計算幾何 平面最近點對 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
大意:給出平面中的兩類點,求平面最近點對,要求該點對中的兩點不是同類點