K-D tree学习笔记

翻译过来就是维护k维信息的树，是一种可以高效处理k维空间信息的数据结构。

一般在算法竞赛中，k=2的情况较多。

考虑对于一维数组，我们想要找到一个y，使得对于给定的x，有|x-y|最小。那么不妨考虑二叉搜索树（就是二分法），取数组的中位数为根，构造一棵树，使得每个点的左儿子小于它，右儿子大于它。那么当x比当前节点小时，走左子树；大，走右子树。

而 K-D tree 就具有二叉搜索树的形态。

建树

如果我们已知k维空间中若干不同点的坐标，要建一棵 K-D tree 。就选择一个维度和一个切割点，使在当前维度上比当前点值小的都归到左子树，大的归到右子树。

为了使树的形态平衡，保证复杂度，我们轮流选择k个维度，且每次在选择切割点时选择该维度上的中位数。这样树高是logn+O(1)的。

时间瓶颈是寻找k维上的中位数，如果每次sort，会 \(O(nlog^2n)\) 。实际上我们只需要考虑我们需要的东西：另中位数在中间，左边的比它小，右边的比它大，单次可以直接调用nth_element优化到 \(O(n)\)，总体 \(O(nlogn)\) .

int tot;
struct kd_tree{
	int lc,rc,v[2],mn[2],mx[2];
}t[maxn];
void pushup(int x){
	for(int i=0;i<2;i++){//取k=2的情况 
		t[x].mn[i]=t[x].mx[i]=t[x].v[i];
		if(t[x].lc){
			t[x].mn[i]=min(t[x].mn[i],t[t[x].lc].mn[i]);
			t[x].mx[i]=max(t[x].mx[i],t[t[x].lc].mx[i]);
		}
		if(t[x].rc){
			t[x].mn[i]=min(t[x].mn[i],t[t[x].rc].mn[i]);
			t[x].mx[i]=max(t[x].mx[i],t[t[x].rc].mx[i]);
		}
	}
} 
void build(int &x,int l,int r,int typ){
	x=++tot;
    int mid=(l+r)>>1;
	nth_element(a+l,a+mid,a+r+1,[typ](node p,node q){return p.v[typ]<q.v[typ];});
	//这里，由于要调用typ，把排序函数写这里比较好 
	t[x].v[0]=a[mid].v[0],t[x].v[1]=a[mid].v[1];
	if(l<mid) build(t[x].lc,l,mid-1,!typ);
	if(mid<r) build(t[x].rc,mid+1,r,!typ);
	pushup(x);
}

插入/删除

如果我们维护的点集会发生变动，此时静态建树的 K-D tree 的复杂度就无法得到保证。而常见的平衡树维护平衡的两个操作，旋转和随机优先级，都不能用到 K-D tree 上。所以我们常用以下方式。

1.根号重构

一种方法是用替罪羊树的重构套路，设置一个平衡因子对 K-D tree 重构。但实际上这样只能保证高度是O(logn)，不是严格的logn+O(1)，所以查询复杂度可能退化。（一般能过）

正确的方法是设定一个阈值B，（每次插入时直接从根节点开始和每个节点比较并向下递归？）当插入次数达到B时暴力重构整棵树。删除时仍用（惰性删除？），当树内删除数量达到B时暴力重构。

因此，当我们取到 \(B=\sqrt{nlogn}\) 时复杂度最优，单次均摊 \(O(\sqrt{nlogn})\)

2.二进制分组

如果无删除操作，这种做法是更优的。维护若干棵大小为 \(2^i\) 的 K-D tree ，满足这些树的大小之和为 n。

插入时，新增一棵大小为 \(2^0\) 的 K-D tree ，不断将相同大小的 K-D tree 合并。实际操作时，可以先将能一起合并的所有树拍扁，只需要重构一次。

总复杂度 \(O(nlog^2n)\)。对于每一位 \(2^i\) ，每次操作到时贡献的节点数为 \(2^i\) ，但又只有 \(\frac{n}{2^i}\) 次能产生贡献。相当于每一位都总共贡献过n个节点，加上建树还有一个log，故总复杂度得证。

int tot=0,top,pol[maxn];
void del(int &x){
	pol[++top]=x;
	t[x]={0,0,0,0,0,0,0,0};
	x=0;
}
int newnode(){
	return top?pol[top--]:++tot;//回收节点省空间
} 
void redo(int &x){
	if(!x) return ;
	a[++cnt]={t[x].v[0],t[x].v[1]};
	redo(t[x].lc),redo(t[x].rc);
	del(x);
}

int main(){
	for(int i=1;i<=n;i++){
		int x,y;
		scanf("%d%d",&x,&y);
		a[cnt=1]={x,y};
		for(int j=0;j<20;j++){
			if(!rt[j]){//没有当前大小的树了，建树 
				build(rt[j],1,cnt,0);
				break;
			}
			else redo(rt[j]);
		}
	}
	return 0;
}

查询

1.矩阵查询

递归。

若当前子树对应的矩形和目标矩形无交点，则不继续搜索；
全部被目标矩形包含，返回整个子树信息；
否则，先判断当前节点合法性，继续递归。

这样复杂度单次 \(O(n^{1-\frac{1}{k}})\) ，证明看oi-wiki

2.邻域查询

即求出平面上一个点的最近/远点。

以求最近点为例。暴搜，对每个子树对应的矩阵设计一个估价函数（如，查询点到这个矩阵的最短距离），启发式搜索，先搜索估价函数更优的子树的答案。如果当前节点估价函数都大于当前答案，直接退出。

这个部分时间复杂度的瓶颈就在于，我们用于判断的估价函数是查询点到这个矩阵的一个计算值。

以估价函数为查询点到这个矩阵的最短距离为例。然而当前点x到子树矩形边界的最短距离mindis——我们理想情况下想去把ans更新成的最优情况，不一定有一个具体存在的点y能使Dis(x,y)=mindis，自然就不能把ans更新到那么优。

所以假设先遍历理想状态下mindis的左子树，遍历后答案也可能比右子树中点的答案劣。还需要看x到右子树矩形边界的最短距离，若小于ans，还得遍历右子树更新答案。

（所以轻轻松松被卡到O(n)......

但随机数据下 K-D tree 求解这个问题为均摊 \(O(logn)\)。（骗分的曙光

//用扫描线会很好维护
//这里使用kd tree，相当于对每个点求出最近/远邻 
#include<bits/stdc++.h>
using namespace std;
const int maxn=5e5+5;
int n,tot,rt,ans=2e9,mx,mn,b[maxn];
struct node{
	int v[2];
}a[maxn];
struct kd_tree{
	int lc,rc,v[2],mn[2],mx[2];
}t[maxn];
void pushup(int x){
	for(int i=0;i<2;i++){
		t[x].mn[i]=t[x].mx[i]=t[x].v[i];
		if(t[x].lc){
			t[x].mn[i]=min(t[x].mn[i],t[t[x].lc].mn[i]);//!!查了俩小时 
			t[x].mx[i]=max(t[x].mx[i],t[t[x].lc].mx[i]);
		}
		if(t[x].rc){
			t[x].mn[i]=min(t[x].mn[i],t[t[x].rc].mn[i]);
			t[x].mx[i]=max(t[x].mx[i],t[t[x].rc].mx[i]);
		}
	}
}
void build(int &x,int l,int r,int typ){
	x=++tot; 
	int mid=(l+r)>>1;
	nth_element(a+l,a+mid,a+r+1,[typ](node p,node q){
		return p.v[typ]<q.v[typ];
	});
	t[x].v[0]=a[mid].v[0],t[x].v[1]=a[mid].v[1];
	if(l<mid) build(t[x].lc,l,mid-1,!typ);
	if(mid<r) build(t[x].rc,mid+1,r,!typ);
	pushup(x);
}
int sol(int x){
	return x==0?2e9:x;
}
int dis(int x,int y){
	return abs(t[x].v[0]-t[y].v[0])+abs(t[x].v[1]-t[y].v[1]);
}
int mxdis(int x,int y){
	int t1=abs(t[x].v[0]-t[y].mn[0])+max(abs(t[x].v[1]-t[y].mn[1]),abs(t[x].v[1]-t[y].mx[1]));
	int t2=abs(t[x].v[0]-t[y].mx[0])+max(abs(t[x].v[1]-t[y].mn[1]),abs(t[x].v[1]-t[y].mx[1]));
	return max(t1,t2);
}
int mndis(int x,int y){
	int t1=(t[x].v[0]<t[y].mn[0])?t[y].mn[0]-t[x].v[0]:((t[x].v[0]>t[y].mx[0])?t[x].v[0]-t[y].mx[0]:0);
	int t2=(t[x].v[1]<t[y].mn[1])?t[y].mn[1]-t[x].v[1]:((t[x].v[1]>t[y].mx[1])?t[x].v[1]-t[y].mx[1]:0);
	return t1+t2;
}
void mxque(int x,int p){
	mx=max(mx,dis(p,x));
	if(!t[x].lc){
		if(t[x].rc) mxque(t[x].rc,p);
	}
	else if(!t[x].rc) mxque(t[x].lc,p);
	else{
		int mx1=mxdis(p,t[x].lc),mx2=mxdis(p,t[x].rc);
		if(mx>=mx1&&mx>=mx2) return ;
		if(mx1>=mx2){
			mxque(t[x].lc,p);
			if(mx<mx2) mxque(t[x].rc,p);
		} else{
			mxque(t[x].rc,p);
			if(mx<mx1) mxque(t[x].lc,p);
		}
	}
}
void mnque(int x,int p){
	mn=min(mn,sol(dis(p,x)));//细节 
	if(!t[x].lc){
		if(t[x].rc) mnque(t[x].rc,p);
	}
	else if(!t[x].rc) mnque(t[x].lc,p);
	else{
		int mn1=mndis(p,t[x].lc),mn2=mndis(p,t[x].rc);
		if(mn<=mn1&&mn<=mn2) return ;
		if(mn1<=mn2){
			mnque(t[x].lc,p);
			if(mn>mn2) mnque(t[x].rc,p);
		} else{
			mnque(t[x].rc,p);
			if(mn>mn1) mnque(t[x].lc,p); 
		}
	}
}
int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;i++){
		scanf("%d%d",&a[i].v[0],&a[i].v[1]);
	}
	for(int i=1;i<=n;i++) b[i]=i;
	build(rt,1,n,0);
	for(int i=1;i<=n;i++){
		mx=0,mn=2e9;
		mxque(rt,i);
		mnque(rt,i);
		ans=min(ans,mx-mn);
	}
	printf("%d",ans);
	return 0;
} 

ps.

1.调 K-D tree 小技巧：阅读程序，重点关注下标和数组

2.实际上如果真要考 K-D tree 的邻域查询，那么考欧氏距离最远邻居多。（但其实这正解好像是旋转卡壳，K-D tree顶多是优秀的骗分手段）

原因如下。

若限制为曼哈顿距离：二维平面最近/远邻，可以直接用cdq分治（不带修改——二维偏序，带修改——三维偏序）解决，\(O(nlog^2n)\) 并不会因为数据影响正确性。例如[Violet] 天使玩偶/SJY摆棋子。

若限制为欧氏距离：二维平面最近邻，即平面最近点对。oiwiki上有专门的介绍，算法也是分治。