KD树与SKD树
首先给出两个搜到的有点内容的KD树文章,论述的比我说的更完整(更冗长),可以先看,也可以看完本文再看.
* https://www.zybuluo.com/l1ll5/note/967681
* http://www.whudj.cn/?p=920
主体思想
- KD树和SKD树都使用坐标轴对齐的最小包围盒来描述空间.
- 例如,平面内,一堆点的点集对应的空间可以用点
A
=(min(all_x),min(all_y))
,B
=(max(all_x),max(all_y))
对应的矩形空间来描述. 当点的个数变为1时,这个矩形空间也会自然地退化为一个点.
- 例如,平面内,一堆点的点集对应的空间可以用点
- 构建时的主要思想: 每个节点
node
都对应了一个空间box(node)
,node->plane
用一个超平面把这个空间box(node)
一分为二,记为sub_left
,sub_right
- 左子树
node->left
对应的空间box(node->left)
包含于sub_left
- 右子树
node->left
对应的空间box(node->right)
包含于sub_right
- 左子树
- 访问时的主要思想: 由于我们知道
box(node->left)
和box(node->right)
,那么我们就可以利用这些box
进行快捷的加速.- 例如, 如果我们需要找到距离
A
距离为5的点, 且我们已知A
到box(node->left)
的最近距离为10,那么,距离为5的点就一定不在node->left
中.
- 例如, 如果我们需要找到距离
注意: KD树和SKD的主体思想基本都是这样,但是具体的实现/应用中会有一些细节上的变化.
补充: 算法中使用的超平面都是垂直于坐标轴的超平面,很容易用数学语言描述,
x_i=K
就是穿过x_i=K
处且与x_i
轴垂直的超平面,例如三维空间中的z=8
KD树
直接上代码, It tells everything. 注意,这份代码主要关注的是可读性,而非性能.
#define DIM_NUM 3
struct Point {
double d[DIM_NUM];
double &operator[](int i) {
return d[i];
}
const double &operator[](int i) const {
return d[i];
}
};
struct Box {
Point lo;
Point hi;
Box() {
for (int i = 0; i < DIM_NUM; ++i) {
lo[i] = INFINITY;
hi[i] = -1.0 * INFINITY;
}
}
};
static std::function<bool(const Point &, const Point &)> GetCmpAtAxis(int axis_ind) {
return [axis_ind](const Point &l, const Point &r) -> bool { return l[axis_ind] < r[axis_ind]; };
}
struct Node {
Node *left;
Node *right;
Point p; // 既用于存储p,又用于作为分割平面
int axis_id;
Box node_box;
};
Box BuildNodeBoxFromPoints(const std::vector<Point> &points,int left,int right) {
Box ret;
for (int j = left; j <= right; ++j>) {
auto & p = points[j];
for (int i = 0; i < DIM_NUM; ++i) {
if (p[i] < ret.lo[i]) {
ret.lo[i] = p[i];
}
if (p[i] > ret.hi[i]) {
ret.hi[i] = p[i];
}
}
}
return ret;
}
int FindPlaneAxis(std::vector<Point> &points, int left, int right) {
// There are many way to determine the axis, here just shows a simple one.
static int i = -1;
++i;
return i % DIM_NUM;
}
// Note: points[right] is access-able
Node *BuildKDTree(std::vector<Point> &points, int left, int right) {
if (left > right) {
return nullptr;
}
int axis_id = FindPlaneAxis(points, left, right);
auto cmp = GetCmpAtAxis(axis_id);
int mid = left + (right - mid) / 2;
std::nth_element(points.begin() + left, points.begin() + mid, points.begin() + right + 1, cmp);
Node *node = new Node;
node->p = points[mid];
node->axis_id = axis_id;
node->node_box = BuildNodeBoxFromPoints(points,left,right);
node->left = BuildKDTree(points, left, mid - 1);
node->right = BuildKDTree(points, mid + 1, right);
return node;
}
- 补充:
- 显然,KD树是一个完全二叉树
- 因此,它是平衡树
- 必要时,可以使用数组来存储节点(第i个节点的子节点分别位于2i和2i+1处)
FindPlaneAxis
的较好实现是选择点集中标准差最大的方向.
- 寻找距离
target
最近点的一种方案:- 先设置
d=inf
作为初始距离. - 如果
PointToBoxDist(target,root->right) <= d
,那么最近点可能会在左子树中,否则一定不在左子树中 - 如果
PointToBoxDist(target,root->right) <= d
,那么最近点可能会在右子树中.否则一定不在右子树中 - 求解过程中不断更新
d
即可,很容易用递归实现.
- 先设置
SKD树: KD树的简单变种.
- 在KD树中,
Point
是树中存储的元素,而在SKD树中,变为了Box
- SKD树中,
Point
可以用退化的Box
表示.
- SKD树中,
- SKD建树过程中,一般用每个
box
的中心作为比较的基准.- 在确定了分割面所在的坐标
axis_point
后,会出现一个问题:axis_point
对应的box
肯定有一半位于分割面左侧,一半位于右侧.对此,SKD的策略是,只要box的中心cent <= axis_point
,那么都划分到左子树,否则为右子树. - SKD树建树的过程中,node的node_box需要能包含所有元素,也就是
node.node_box.lo=min(all_box.lo);node.node_box.hi=max(all_box.hi)
- 显然,按照这种策略实现时,只有叶节点会存储元素,非叶节点仅用于分割
- 在确定了分割面所在的坐标
- SKD树中,左子树和右子树对应的
node_box
有所重叠,但这并不影响效率.
附上: SKD 树的原始论文
Spatial KD tree A Data Structure for Geographic Database