KD-Tree学习小记

最近进入了赛后养生模式，感觉身体被掏空，还要补课···学点养生数据结构。

引入：正交区域查找

正交区域的精确定义我不知道…而一般来说，一个空间内查询一个特定长方体里面的信息，一个平面内查询某个长方形的查找都可以称作正交区域查找。而二维平面的不规则的区域查找我们可以转化为多个正交区域查找的并。
而许多信息的查询是可以转化为正交区域查找的，例如问一堆员工中，年龄在[a,b]，工资在[l,r]中的有几个。这个有很多做法，什么树套树之类的。而一种思路是，把员工的年龄x和工资y映射到二维平面上的点(x,y)上，这样就可以进行正交区域查找了，即查找一个矩形中点的个数。
我们需要一个数据结构，一个可以在任何维数下使用的数据结构。（树套树只适用二维，他和一维的树代码差别还是挺大的）

利器：KD树

先解释一下名字，K是维数，D是Dimension，即维。“树”表明他是树的结构。
基本地，KD树中一个节点储存了：

• K维空间域，例如三维中的一个长方体
• 一个K维点的坐标
• 两个儿子下标
  在平衡树中，我们知道：可以维护以每个节点为根的子树权值的min和max。
  K维空间域与此很类似，维护的是子树点的坐标范围。

const int K=2;
struct KD_Tree
{
    int d[K],son[2];
    int x[2],y[2];//Range[K][2];
}tr[N];

如上代码，d为节点储存的原图的点坐标，son为儿子，第二行储存了K维空间域。

构造

大体思想：
KD树是一颗平衡二叉树，其中每个非叶节点，可以想象一个超平面，用来分割其储存的空间域，其中超平面垂直于坐标轴。
主线：

• 树尽量平衡，超平面划分的两个空间内的点尽量一样多。
• 为了有扩展性，树的每一层的超平面垂直的坐标轴，要轮流来取。即第一层垂直x轴，第二层垂直y轴，第三层垂直z轴····

垂直某个轴，意味着以这个轴的坐标为关键字来操作。
例如这次要垂直x轴，我们取当前点集的x坐标的中位数，然后把它作为切分点，切分点作为父节点，即KD树中新节点储存的点；切开的两边的点分别属于左右子树的点集。

这是一个例子

下面给出构造代码

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
#define fo(i,j,k) for(i=j;i<=k;i++)
#define fd(i,j,k) for(i=j;i>=k;i--)
#define cmax(a,b) (a<b?a=b:a)
#define cmin(a,b) (a>b?a=b:a)
struct P
{
    int d[K];
}a[N];
struct KD_Tree
{
    int d[K],son[2];
    int x[2],y[2];//Range[K][2];
}tr[N];
int D,root;
int i,n;
bool operator <(P a,P b)
{
    return a.d[D]<b.d[D];
}
void mt(int x,int y)
{
    cmin(tr[x].x[0],tr[y].x[0]);
    cmax(tr[x].x[1],tr[y].x[1]);
    cmin(tr[x].y[0],tr[y].y[0]);
    cmax(tr[x].y[1],tr[y].y[1]);
}
int build(int l,int r,int d)
{
    D=d;
    d^=1;
    int m=(l+r)/2;
    nth_element(a+l,a+m,a+r+1);
    tr[m].d[0]=tr[m].x[0]=tr[m].x[1]=a[m].d[0];
    tr[m].d[1]=tr[m].y[0]=tr[m].y[1]=a[m].d[1];
    if (l<m) mt(m,tr[m].son[0]=build(l,m-1,d));
    if (m<r) mt(m,tr[m].son[1]=build(m+1,r,d));
    return m;
}
int main()
{
    scanf("%d",&n);
    fo(i,1,n) scanf("%d %d",&a[i].d[0],&a[i].d[1]);
    root=build(1,n,0);
}

nth_element 就是找出第n大的元素，把它放到第n位，同时它左边的元素都比它大，右边的都比它小。其实就是快速选择。
用这个构造以后，我们就可以很方便的正交区域查找了，二维时，时间复杂度为$O(n*sqrt(n))$，与树套树同阶。分析详见清华大学出版社的《计算几何——算法与应用》。

估价

这是用于优化查询复杂度的东西。
这里估价估的是目标点到当前查询区域距离的上界和下界。
设a为查询点，l,r表示域的某个坐标的min和max。
那么曼哈顿距离最小

max(0,lx-ax)+max(0,ax-rx)+max(0,ly-ay)+max(0,ay-ry);

每个坐标分三种情况讨论，发现是这样的。
而曼哈顿距离最大很类似。

max(abs(lx-ax),abs(rx-ax))+max(abs(ly-ay),abs(ry-ay));

接下来还有欧几里得距离，我们不开方，也分三种情况讨论。很类似就不写了。
切比雪夫距离？转45°就成了曼哈顿。

修改

离线就一样的嘛。只要在每个节点处标记一下是激活了还是未激活。
在线的话，像替罪羊树弄一个平衡因子fac，只要siz[x]*fac<=max(siz[son[x]])，就暴力重构。
这样单次修改时间复杂度均摊logn。

查询

查询与(x,y)最近的点的思路是这样的：
我们从根KD-Tree，对每个节点x，我们先用它更新答案，然后再对两个儿子做一次估价，估价值优的优先访问。当然，如果估价值比答案要劣，那就不用走了，因为估价函数是最优的情况了，最优都更劣，那就不用管了。
随机数据下，这样找一次是O(log n)的，而构造可以卡成O(sqrt(n))。
构造之后讲。

题目

代码里有一些新学的，比较优美的，优化代码排版的方法。
1,【SDOI2010】Hide and Seek

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
#define fo(i,j,k) for(i=j;i<=k;i++)
#define fd(i,j,k) for(i=j;i>=k;i--)
#define cmax(a,b) (a<b?a=b:a)
#define cmin(a,b) (a>b?a=b:a)
typedef long long ll;
typedef double db;
const int N=100005,K=2;
struct P
{
    int d[K];
}a[N];
struct KD_Tree
{
    int d[K],son[2];
    int x[2],y[2];//Range[K][2];
}tr[N];
int D,root;
int ans,tp[2],ax,ay,lx,ly,rx,ry,p;
int i,n;
bool operator <(P a,P b)
{
    return a.d[D]<b.d[D];
}
void mt(int x,int y)
{
    cmin(tr[x].x[0],tr[y].x[0]);
    cmax(tr[x].x[1],tr[y].x[1]);
    cmin(tr[x].y[0],tr[y].y[0]);
    cmax(tr[x].y[1],tr[y].y[1]);
}
int build(int l,int r,int d)
{
    D=d;
    d^=1;
    int m=(l+r)/2;
    nth_element(a+l,a+m,a+r+1);
    tr[m].d[0]=tr[m].x[0]=tr[m].x[1]=a[m].d[0];
    tr[m].d[1]=tr[m].y[0]=tr[m].y[1]=a[m].d[1];
    if (l<m) mt(m,tr[m].son[0]=build(l,m-1,d));
    if (m<r) mt(m,tr[m].son[1]=build(m+1,r,d));
    return m;
}
int calc(int tar,int x,int sig)
{
    if (!x) return -1;
    lx=tr[x].x[0];rx=tr[x].x[1];
    ly=tr[x].y[0];ry=tr[x].y[1];
    if (sig>0) return max(0,lx-ax)+max(0,ax-rx)+max(0,ly-ay)+max(0,ay-ry);
    else      return max(abs(lx-ax),abs(rx-ax))+max(abs(ly-ay),abs(ry-ay));
}
void get(int tar,int x,int d,int sig)
{
    D=d;
    d^=1;
    if (x!=tar) 
    {
        int tmp=abs(ax-tr[x].d[0])+abs(ay-tr[x].d[1]);
        if (sig>0) tp[0]=min(tp[0],tmp);
        else       tp[1]=max(tp[1],tmp);
    }
    int f1=calc(tar,tr[x].son[0],sig),t1=tr[x].son[0];
    int f2=calc(tar,tr[x].son[1],sig),t2=tr[x].son[1];
    if (f1*sig>f2*sig||!t1) swap(f1,f2),swap(t1,t2);
    if (t1&&f1*sig<tp[p]*sig)
        get(tar,t1,d,sig);
    if (t2&&f2*sig<tp[p]*sig)
        get(tar,t2,d,sig);
}
int main()
{
    freopen("1519.in","r",stdin);
//  freopen("1519.out","w",stdout);
    scanf("%d",&n);
    fo(i,1,n) scanf("%d %d",&a[i].d[0],&a[i].d[1]);
    root=build(1,n,0);
    ans=1e9;
    fo(i,1,n)// enumerating indexes of tr
    {
        tp[0]=1e9;
        tp[1]=0;
        ax=a[i].d[0];
        ay=a[i].d[1];
        p=0;
        get(i,root,0,1);
        p=1;
        get(i,root,0,-1);
        ans=min(ans,tp[1]-tp[0]);
    }
    printf("%d\n",ans);
}

参考

Crazy的计算几何PPT
n+e的教学PDF《K-D Tree在信息学竞赛中的应用》
清华大学出版社《计算几何——算法与应用》