KD樹小結

很久之前我就想過怎么快速在二維平面上查找一個區域的信息，思考許久無果，只能想到幾種優秀一點的暴力。

KD樹就是干上面那件事的。

別的不多說，趕緊把自己的理解寫下來，免得涼了。

?

KD樹的組成

以維護k維空間(x,y,……)內的KD樹為例，主要由一下三部分組成：

1. p[k]，代表樹上這個結點所儲存的點(在題目中給出的/你自己加上的點集中的一個點)。
2. ch[2]，表示它的子結點(沒錯，KD樹是一棵二叉樹)
3. mi[k]與mx[k]，mi/mx[i]代表KD樹這個結點統轄的所有點的第i-1范圍。比如說mi[1]=2,mx[1]=4，就代表這棵樹統轄的點的y坐標都在[2,4]內。

?

不看mi和mx，長得就和splay/trie樹一樣，一個p維護當前節點，一個ch[2]記錄左右兒子。

不看p[k]，長得就和線段樹一樣，有左右兒子和區間信息。

沒錯，KD樹功能如線段樹，結點維護區域信息；形態如splay/trie樹，每個結點有實際的值和意義。

?

KD樹的構建

一般題目都是二維平面。下面就以二維平面KD樹的構建為例。

讀入把點存進結構體數組a中，坐標分別為a[x].p[i]。

inline void build(int &x,int l,int r,int type){x=(l+r)>>1;now=type;nth_element(a+l,a+x,a+r+1,cmp);nd=a[x];newnode(x);if(l<x)build(ch[x][0],l,x-1,type^1);else ch[x][0]=0;if(x<r)build(ch[x][1],x+1,r,type^1);else ch[x][1]=0;pushup(x);
}build(kd.root,1,n,0);

非常優美……對type、now作用不明的同學請繼續閱讀……你要現在就明白就奇怪了

系統函數nth_element(a+l,a+x,a+r+1)，頭文件algorithm，需定義＜或cmp函數。

作用：把排序后第x大的放到第x位，比它小的放進左邊，比它大的放進右邊(兩邊無序)。

注意區間開閉：左閉右開，中間也是閉合的。

復雜度：平均，期望是O(n)？可以接受。

下面給出cmp、newnode、pushup代碼。

struct Node{int p[2],mi[2],mx[2];}a[N];
inline bool cmp(const Node &a,const Node &b){return a.p[now]<b.p[now];}
inline void Min(int &x,int y){x=x<y?x:y;}
inline void Max(int &x,int y){x=x>y?x:y;}
inline void pushup(int x){int ls=ch[x][0],rs=ch[x][1];if(ls){Min(T[x].mi[0],T[ls].mi[0]);Max(T[x].mx[0],T[ls].mx[0]);Min(T[x].mi[1],T[ls].mi[1]);Max(T[x].mx[1],T[ls].mx[1]);}if(rs){Min(T[x].mi[0],T[rs].mi[0]);Max(T[x].mx[0],T[rs].mx[0]);Min(T[x].mi[1],T[rs].mi[1]);Max(T[x].mx[1],T[rs].mx[1]);}
}inline void newnode(int x){T[x].p[0]=T[x].mi[0]=T[x].mx[0]=nd.p[0];T[x].p[1]=T[x].mi[1]=T[x].mx[1]=nd.p[1];
}

不要問我為什么辣么長，為了減常沖榜，把循環展開了……

聰明的讀者已經發現KD樹的構建巧妙之處。它不是純粹按照x維，或者某一維排序，而是先按x維，再按y維，再按z維，再……最后又回到x維……

這樣分割的區域更加整齊劃一更加均勻緊縮，不像上面的按照某一維劃分，到最后變成一條條長條，KD樹劃分到底的圖案還是很好看的。

這樣分割有什么好處呢？等你真正領悟了KD樹的精髓之后你就會發現……嘿嘿嘿……

(就是為了把這個暴力數據結構剪枝剪更多跑更快)

?

KD樹的操作

1.往KD樹上插點

插點可以分為插新點和插老點。如果有老點，特判一句，把信息覆蓋即可。

inline void insert(int &x,int type){if(!x){x=++cnt,newnode(cnt);return;}if(nd.p[0]==T[x].p[0] && nd.p[1]==T[x].p[1]){……(自行維護);return;}if(nd.p[type]<T[x].p[type])insert(ch[x][0],type^1);else insert(ch[x][1],type^1);pushup(x);
}

依然非常的美妙……等等有什么不對？

我們能估計出一棵剛建好的KD樹深度是O(log)的。

但你這么隨便亂插……有道題叫HNOI2017 spaly 插入不旋轉的單旋spaly見過？T成茍。

這都不是問題！知不知道有一種數據結構叫做替罪羊樹哇？

知道替罪羊樹怎么保證復雜度的嗎？

重構！大力重構！自信重構！不爽就重構！

為了省事大概每插入10000次就重構一次好了……

if(kd.cnt==sz){for(int i=1;i<=sz;++i)a[i]=kd.T[i];kd.rebuild(kd.root,1,sz,0);sz+=10000;
}

?

2.在KD樹上查詢

• 如果是單點(給定點)查詢：
  • 太簡單啦！把插入改改就闊以辣！
• 如果是查詢距離一個點(x',y')最近的點(曼哈頓距離，|x-x'|+|y-y'|)：
  • 首先我們看暴力的剪枝：按某一維排序，如果該維的差過大就不管了。
  • 而令我們期待的KD樹呢？呃不好意思，它也是這么做的……
  • 我們維護過兩個叫做mi[]和mx[]的東西吧……這個時候就是它派上用場了。
  • 具體還請看代碼吧：

    //查詢的點(x',y')儲存在nd中。
    //這里的l,r就是mi,mx的意思。
    inline int dis(Node p,int x,int ans=0){for(int i=0;i<2;++i)ans+=max(0,t[x].l[i]-p.p[i])+max(0,p.p[i]-t[x].r[i]);return ans;
    }inline void query(int x){Ans=min(Ans,abs(t[x].p[0]-nd.p[0])+abs(t[x].p[1]-nd.p[1]));int dl=ch[x][0]?dis(nd,ch[x][0]):Inf;int dr=ch[x][1]?dis(nd,ch[x][1]):Inf;if(dl<dr){if(dl<Ans)query(ch[x][0]);if(dr<Ans)query(ch[x][1]);}else{if(dr<Ans)query(ch[x][1]);if(dl<Ans)query(ch[x][0]);}
    }

  • dis():如果當前點在這個區間內就是0，否則就是最極的點到它的距離。
  • 聰明絕頂的你已經發現了……這TM就是個暴力。
  • 其實這個數據結構就是一個暴力……
  • 當暴力有了時間復雜度證明……還叫暴力么？讀書人的事，能叫偷么？
  • 這么暴力有幾個好處：不用枚舉所有點；剪枝有效及時。
  • 復雜度有保障，大概在O(√n)級別下，主要看數據。
• 如果是區間查詢，以區間查詢點權和為例(之前就有維護好)：

  • inline bool in(int l,int r,int xl,int xr){return l<=xl && xr<=r;}
    inline bool out(int l,int r,int xl,int xr){return xr<l || r<xl;}inline int query(int x,int x1,int y1,int x2,int y2){int ans=0;if(!x)return ans;if(in(x1,x2,T[x].mi[0],T[x].mx[0]))if(in(y1,y2,T[x].mi[1],T[x].mx[1]))return T[x].sum;if(out(x1,x2,T[x].mi[0],T[x].mx[0]))return 0;if(out(y1,y2,T[x].mi[1],T[x].mx[1]))return 0;if(in(x1,x2,T[x].p[0],T[x].p[0]))if(in(y1,y2,T[x].p[1],T[x].p[1]))ans+=T[x].val;return ans+query(ch[x][0],x1,y1,x2,y2)+query(ch[x][1],x1,y1,x2,y2);
    }

  • 別看代碼長又看起來復雜，寫起來跟線段樹似的，還是一樣的暴力搞。

?

KD樹的基本姿勢大概就是這個樣子……好寫不好寫錯，基本上都是個板子。

附上學長的一(三)句話，從各方面進行了深度總結：

“能不寫最好還是不要寫吧，輕松被卡→_→，也許可以出奇制勝？如果要寫，重新構樹是個不錯的選擇。發現大數據跑不過，多半是剪枝掛了。”

還是給個鏈接……MashiroSky大爺。

upd:以當前坐標差最大的來做type應該比輪換type更優秀……

例題有"SJY擺棋子"、"簡單題"等，在此就不做贅述了。

比較有意思的應用就是【bzoj3489】 A simple rmq problem，正如上面所言，KD樹解決傳統數據結構題出奇制勝。

附上"BZOJ4066簡單題"代碼一份，操作是單點修改+矩形求和在線。

?

#include    <iostream>
#include    <cstdio>
#include    <cstdlib>
#include    <algorithm>
#include    <vector>
#include    <cstring>
#include    <queue>
#include    <complex>
#include    <stack>
#define LL long long int
#define dob double
#define FILE "bzoj_4066"
//#define FILE "簡單題"
using namespace std;const int N = 200010;
int n,lst,now,sz=10000;inline int gi(){int x=0,res=1;char ch=getchar();while(ch>'9'||ch<'0'){if(ch=='-')res*=-1;ch=getchar();}while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();return x*res;
}inline void Min(int &x,int y){x=x<y?x:y;}
inline void Max(int &x,int y){x=x>y?x:y;}
struct Node{int p[2],mi[2],mx[2],val,sum;}a[N];
inline bool cmp(const Node &a,const Node &b){return a.p[now]<b.p[now];}
struct KD_Tree{int ch[N][2],root,cnt;Node T[N],nd;inline void pushup(int x){int ls=ch[x][0],rs=ch[x][1];if(ls){Min(T[x].mi[0],T[ls].mi[0]);Max(T[x].mx[0],T[ls].mx[0]);Min(T[x].mi[1],T[ls].mi[1]);Max(T[x].mx[1],T[ls].mx[1]);}if(rs){Min(T[x].mi[0],T[rs].mi[0]);Max(T[x].mx[0],T[rs].mx[0]);Min(T[x].mi[1],T[rs].mi[1]);Max(T[x].mx[1],T[rs].mx[1]);}T[x].sum=T[x].val;if(ls)T[x].sum+=T[ls].sum;if(rs)T[x].sum+=T[rs].sum;}inline void newnode(int x){T[x].p[0]=T[x].mi[0]=T[x].mx[0]=nd.p[0];T[x].p[1]=T[x].mi[1]=T[x].mx[1]=nd.p[1];T[x].sum=T[x].val=nd.val;}inline void insert(int &x,int type){if(!x){x=++cnt,newnode(cnt);return;}if(nd.p[0]==T[x].p[0] && nd.p[1]==T[x].p[1]){T[x].val+=nd.val;T[x].sum+=nd.val;return;}if(nd.p[type]<T[x].p[type])insert(ch[x][0],type^1);else insert(ch[x][1],type^1);pushup(x);}inline void rebuild(int &x,int l,int r,int type){x=(l+r)>>1;now=type;nth_element(a+l,a+x,a+r+1,cmp);nd=a[x];newnode(x);if(l<x)rebuild(ch[x][0],l,x-1,type^1);else ch[x][0]=0;if(x<r)rebuild(ch[x][1],x+1,r,type^1);else ch[x][1]=0;pushup(x);}inline bool in(int l,int r,int xl,int xr){return l<=xl && xr<=r;}inline bool out(int l,int r,int xl,int xr){return xr<l || r<xl;}inline int query(int x,int x1,int y1,int x2,int y2){int ans=0;if(!x)return ans;if(in(x1,x2,T[x].mi[0],T[x].mx[0]))if(in(y1,y2,T[x].mi[1],T[x].mx[1]))return T[x].sum;if(out(x1,x2,T[x].mi[0],T[x].mx[0]))return 0;if(out(y1,y2,T[x].mi[1],T[x].mx[1]))return 0;if(in(x1,x2,T[x].p[0],T[x].p[0]))if(in(y1,y2,T[x].p[1],T[x].p[1]))ans+=T[x].val;return ans+query(ch[x][0],x1,y1,x2,y2)+query(ch[x][1],x1,y1,x2,y2);}}kd;int main()
{freopen(FILE".in","r",stdin);freopen(FILE".out","w",stdout);n=gi();while(1){int type=gi();if(type==3)break;int x1=gi()^lst,y1=gi()^lst;if(type==1){int A=gi()^lst;kd.nd.p[0]=x1;kd.nd.p[1]=y1;kd.nd.sum=kd.nd.val=A;kd.insert(kd.root,0);if(kd.cnt==sz){for(int i=1;i<=sz;++i)a[i]=kd.T[i];kd.rebuild(kd.root,1,sz,0);sz+=10000;}}if(type==2){int x2=gi()^lst,y2=gi()^lst;lst=kd.query(kd.root,x1,y1,x2,y2);printf("%d\n",lst);}}fclose(stdin);fclose(stdout);return 0;
}

?

?

?