[学习笔记]K-D Tree

推荐：

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]就好。