[學習筆記]K-D Tree
推薦:
對於D維的點若幹,多次查詢距離某個點第K大的點是什麽。
處理這一類問題的一個數據結構,叫K-D Tree
基本思想是對點進行區域分塊處理。
圖示:
K-D Tree是一個二叉樹。
每個點維護的信息是,
split :分裂坐標軸
ls、rs:左右兒子
node:該節點存儲的真實點
建樹:
遞歸建樹。類似線段樹(但是每個點有實際的點)
選擇當前區域的點的各個維度的方差最大的維度(傳說如果方差大,數據分散,復雜度或者精度有所保證??),把這個維度當做split
這個節點的真實點就是c[mid]
然後,把這個維度[s]小於c[mid][s]的放在左邊,大於的放在右邊。
(實現時,用一個nth_element,再重載小於號,可以O(n)實現把中間的放在mid位置上,並且這個維度[s]小於c[mid][s]的放在左邊,大於的放在右邊。)
然後遞歸建樹即可。
x,y是split
這樣,整個K-D Tree就把一些點分成了若幹個塊。
我們一塊一塊處理會比較容易剪枝。
查詢:最近的點(即K=1)
本質是爆搜+剪枝。。。
設查詢距離點st的最近的點to
設距離為now
法一:
不斷通過當前split維度和st這個維度的大小比較,
我們先走st所屬的塊,
回溯回來之後,
由於可能在另一半有更近的點。
如果分界線到st的該維度距離小於now
那麽再走另外一個塊搜索。
法二:
上面那個剪枝比較粗糙。
我們發現,一個塊的所有點,其實可以用一個矩形框住。
那麽,如果st到這個矩形可能的最近點距離小於now的話,再搜下去。
具體來說,我們每個節點維護這個節點代表的塊內,最大最小的x,y坐標。(其實就是矩形四個頂點)
這個最短距離:
x差為:dx=max(st.x-x.mxx,0)+max(x.min-st.x,0)
畫個圖理解下,如果st的x在mix,mxx之間的話,那麽x差就認為是0
如果在mix左邊,那麽就是x-mix,
如果在mxx右邊,那麽就是mxx-x
y同理。
dis=sqrt(dx*dx+dy*dy)
發現,當st在所屬的塊內時,dis一定是0
然後就可以剪枝了。
對於兩個兒子,選擇估價距離較小的那個先搜,回溯時,如果另一個距離還比now小的話,再搜另一個。
理論上,應該比法一多減一些枝。
總之,復雜度不明。
傳說最差O(k*n^(1-1/k))每次(k是維度)
例題:
平面最近點對(加強版)
這個題用分治是最好的。
我們用KD樹來試試。
枚舉所有的點,找到與它最近的點距離,然後所有距離取min即可
直接做就好了。
註意:
1.如果寫法二,那麽對於0號節點的哨兵必須mi=inf,mx=-inf。否則剪枝就掛了。
2.建樹的時候,build返回節點編號,不能返回計數器tot。。。。
這個題法一更快??
可能數據水,然後法二常數大吧。。。
法一:
#include<bits/stdc++.h>
#define reg register int
#define il inline
#define numb (ch^‘0‘)
using namespace std;
typedef long long ll;
il void rd(int &x){
char ch;x=0;bool fl=false;
while(!isdigit(ch=getchar()))(ch==‘-‘)&&(fl=true);
for(x=numb;isdigit(ch=getchar());x=x*10+numb);
(fl==true)&&(x=-x);
}
namespace Miracle{
const int N=200000+5;
const double inf=2333333333.00;
int n;
struct po{
double x,y;
int id;
po(){}
po(double xx,double yy){
x=xx;y=yy;
}
}a[N],c[N],st,to;
bool cmp1(po a,po b){
return a.x<b.x;
}
bool cmp2(po a,po b){
return a.y<b.y;
}
double ans;
double now;
double dis(po a,po b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
struct tr{
int sp;
po O;
int ls,rs;
}t[2*N];
int tot;
int rt;
int build(int l,int r){
if(l>r){
return 0;
}
if(l==r){
++tot;
t[tot].O=c[l];
t[tot].ls=t[tot].rs=0;
t[tot].sp=1;
return tot;
}
int mid=(l+r)>>1;
double ax=0,ay=0;
for(reg i=l;i<=r;++i) ax+=c[i].x,ay+=c[i].y;
ax/=(r-l+1);ay/=(r-l+1);
double fx=0,fy=0;
for(reg i=l;i<=r;++i) fx+=(c[i].x-ax)*(c[i].x-ax),fy+=(c[i].y-ay)*(c[i].y-ay);
fx/=(r-l+1);fy/=(r-l+1);
int ret=++tot;
if(fx>fy){//choose x;
nth_element(c+l,c+mid,c+r+1,cmp1);
t[ret].sp=1;
}else{
nth_element(c+l,c+mid,c+r+1,cmp2);
t[ret].sp=2;
}
t[ret].O=c[mid];
t[ret].ls=build(l,mid-1);
t[ret].rs=build(mid+1,r);
return ret;
}
void dfs(int x){
if(!x) return;
if(st.id!=t[x].O.id&&dis(st,t[x].O)<now){
now=dis(st,t[x].O);
}
if(t[x].sp==1){
double d=fabs(t[x].O.x-st.x);
if(st.x<=t[x].O.x){
dfs(t[x].ls);
if(d<now) dfs(t[x].rs);
}
else{
dfs(t[x].rs);
if(d<now) dfs(t[x].ls);
}
}
else{
double d=fabs(t[x].O.y-st.y);
if(st.y<=t[x].O.y){
dfs(t[x].ls);
if(d<now) dfs(t[x].rs);
}
else{
dfs(t[x].rs);
if(d<now) dfs(t[x].ls);
}
}
}
int main(){
scanf("%d",&n);
for(reg i=1;i<=n;++i){
scanf("%lf%lf",&a[i].x,&a[i].y);
a[i].id=i;
c[i]=a[i];
}
rt=build(1,n);
ans=inf;
for(reg i=1;i<=n;++i){
st=a[i];
now=inf;
to=po(inf,inf);
dfs(1);
ans=min(ans,now);
}
printf("%.4lf",ans);
return 0;
}
}
int main(){
Miracle::main();
return 0;
}
/*
Author: *Miracle*
Date: 2018/11/26 8:43:17
*/
法一
法二:
#include<bits/stdc++.h>
#define reg register int
#define il inline
#define numb (ch^‘0‘)
using namespace std;
typedef long long ll;
il void rd(int &x){
char ch;x=0;bool fl=false;
while(!isdigit(ch=getchar()))(ch==‘-‘)&&(fl=true);
for(x=numb;isdigit(ch=getchar());x=x*10+numb);
(fl==true)&&(x=-x);
}
namespace Miracle{
const int N=200000+5;
const double inf=2333333333.00;
int n;
struct po{
double x,y;
int id;
po(){}
po(double xx,double yy){
x=xx;y=yy;
}
}a[N],c[N],st,to;
bool cmp1(po a,po b){
return a.x<b.x;
}
bool cmp2(po a,po b){
return a.y<b.y;
}
double ans;
double now;
double dis(po a,po b){
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
struct tr{
double mxx,mix,mxy,miy;
int sp;
po O;
int ls,rs;
}t[2*N];
int tot;
int rt;
int build(int l,int r){
if(l>r){
return 0;
}
if(l==r){
++tot;
t[tot].mxx=t[tot].mix=c[l].x;
t[tot].mxy=t[tot].miy=c[l].y;
t[tot].O=c[l];
t[tot].ls=t[tot].rs=0;
t[tot].sp=1;
return tot;
}
int mid=(l+r)>>1;
double ax=0,ay=0;
for(reg i=l;i<=r;++i) ax+=c[i].x,ay+=c[i].y;
ax/=(r-l+1);ay/=(r-l+1);
double fx=0,fy=0;
for(reg i=l;i<=r;++i) fx+=(c[i].x-ax)*(c[i].x-ax),fy+=(c[i].y-ay)*(c[i].y-ay);
fx/=(r-l+1);fy/=(r-l+1);
int ret=++tot;
if(fx>fy){//choose x;
nth_element(c+l,c+mid,c+r+1,cmp1);
t[ret].sp=1;
}else{
nth_element(c+l,c+mid,c+r+1,cmp2);
t[ret].sp=2;
}
t[ret].O=c[mid];
t[ret].ls=build(l,mid-1);
t[ret].rs=build(mid+1,r);
t[ret].mxx=max(t[t[ret].rs].mxx,t[t[ret].ls].mxx);
t[ret].mix=min(t[t[ret].rs].mix,t[t[ret].ls].mix);
t[ret].mxy=max(t[t[ret].rs].mxy,t[t[ret].ls].mxy);
t[ret].miy=min(t[t[ret].rs].miy,t[t[ret].ls].miy);
//cout<<" ret "<<ret<<" "<<l<<" "<<r<<endl;
return ret;
}
void dfs(int x){
if(st.id!=t[x].O.id&&dis(st,t[x].O)<now){
now=dis(st,t[x].O);
to=t[x].O;
}
if(t[x].ls&&t[x].rs){
double lx=max(st.x-t[t[x].ls].mxx,0.0)+max(t[t[x].ls].mix-st.x,0.0);
double ly=max(st.y-t[t[x].ls].mxy,0.0)+max(t[t[x].ls].miy-st.y,0.0);
double len1=sqrt(lx*lx+ly*ly);
double rx=max(st.x-t[t[x].rs].mxx,0.0)+max(t[t[x].rs].mix-st.x,0.0);
double ry=max(st.y-t[t[x].rs].mxy,0.0)+max(t[t[x].rs].miy-st.y,0.0);
double len2=sqrt(rx*rx+ry*ry);
if(len1<=len2&&len1<now){
dfs(t[x].ls);
if(len2<now)
dfs(t[x].rs);
}
else if(len2<=len1&&len2<now){
dfs(t[x].rs);
if(len1<now)
dfs(t[x].ls);
}
}
else if(t[x].ls){
double lx=max(st.x-t[t[x].ls].mxx,0.0)+max(t[t[x].ls].mix-st.x,0.0);
double ly=max(st.y-t[t[x].ls].mxy,0.0)+max(t[t[x].ls].miy-st.y,0.0);
double len1=sqrt(lx*lx+ly*ly);
if(len1<now)
dfs(t[x].ls);
}
else if(t[x].rs){
double rx=max(st.x-t[t[x].rs].mxx,0.0)+max(t[t[x].rs].mix-st.x,0.0);
double ry=max(st.y-t[t[x].rs].mxy,0.0)+max(t[t[x].rs].miy-st.y,0.0);
double len2=sqrt(rx*rx+ry*ry);
if(len2<now)
dfs(t[x].rs);
}
else return;
}
int main(){
scanf("%d",&n);
t[0].mix=inf;t[0].mxx=-inf;
t[0].miy=inf;t[0].mxy=-inf;
for(reg i=1;i<=n;++i){
scanf("%lf%lf",&a[i].x,&a[i].y);
a[i].id=i;
c[i]=a[i];
}
rt=build(1,n);
ans=inf;
for(reg i=1;i<=n;++i){
// cout<<" ii "<<i<<" : "<<a[i].x<<" "<<a[i].y<<" ------------------ "<<endl;
st=a[i];
now=inf;
to=po(inf,inf);
dfs(1);
// cout<<" after "<<now<<endl;
ans=min(ans,now);
}
printf("%.4lf",ans);
return 0;
}
}
int main(){
Miracle::main();
return 0;
}
/*
Author: *Miracle*
Date: 2018/11/26 8:43:17
*/
法二
對了,KD-Tree其實也可以不記錄左右兒子,以及代表實際點
因為,每次我們選擇的是mid位置的點,之後這個點的位置也不會再動了。
而左右兒子區間也是定值。
所以,query時記錄(l,r)即可,訪問實際點的話,直接取c[mid]就好。
[學習筆記]K-D Tree