KD树与SKD树

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的点, 且我们已知Abox(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表示.
    • 最小的体素是Box,划分空间时会用新的大Box来包络这些有体积的体素Box
  • SKD建树过程中,一般用box的中心作为比较的基准.
    • 在确定了分割面所在的坐标axis_point后,一般认为,只要box的中心cent <= axis_point,那么都划分到左子树,否则为右子树.
    • SKD树中,左子树和右子树的Box一般总会有重叠的部分,这和KD树是不同的.但这并不影响效率.
    • SKD树中,非叶节点一般仅存储划分点的坐标,划分点并不是一个已经存在的体素,而是一个抽象的点. 也就是说,非叶节点并不会对应到某个体素.

附上: SKD 树的原始论文
Spatial KD tree A Data Structure for Geographic Database ( 访问密码:gZixPB )

Read more

Coroutine

从一般概念上说, 协程是特殊的函数调用: 被调用的函数可以在可控的位置被中断,然后在下一次调用时,继续从上次中断的位置继续执行。 本文主要通过Python的协程来介绍协程, 这是我唯一熟悉的一种协程实现. Classic Coroutine 下面的python代码很好的说明了协程的核心功能 def co_routine(): recv0 = yield 996 # hangs here after first coro.send assert recv0 == "Second" yield 711 # hangs here after second coro.send return def main(): coro = co_routine() # Create a new coroutine object value = coro.send(None)

By Edimetia3D

GDB with Python

这篇文章的主要应用场景是调试Python的C/C++ Extension 1. 同时使用pdb / gdb 进行调试. 通俗点说, 既可以break在 .py 文件中,也可以break在 .cc 文件中 2. 在gdb中不但可以获得常规的调试信息, 还可以获得python VM 的调试信息, 例如获得python的调用栈, 访问Python局部变量等. 这将会在调试exception时(如Segmentfalut)非常有用, 这种场景下, 定位 Python VM 正运行到哪一行代码往往可以提供一些直观的重要信息. 第一步: 编译源码以获得一些辅助数据. 我们并不真的需要使用从源码编译的Python, 但是一些调试相关的辅助文件需要从源码中获得, 包括 python-gdb.py及debug symbol等. 在 https://www.python.org/ftp/python/ 或 https://github.com/python/cpython

By Edimetia3D

Bazel Notes

这是一篇2019年左右的记录, 内容可能过时, 也不太全面 杂谈 Bazel是Google为Monorepo服务而开发的构建工具. 首先是巨大,当问题的规模变大,事情总是会变得更复杂. 而Google面对的"巨大Monorepo",应该是世间罕有的. 然后是Monorepo,这极大的影响了代码的组织风格.例如,你要写一个操作系统内核ProjectOS,还要写一个游戏ProjectGame.在传统的开发习惯中,这两个项目会组织到两个不同的Repo里,PorjectOS和ProjectGame之间无法直接相互引用,例如,你在ProjectOS里写了一个高级的数据结构,想要在Game里也使用,要么直接复制粘贴,要么是创建一个新的CommonRepo,把可公用的代码都放在Common里,然后两个项目各自引入Common作为依赖. 使用MonoRepo则不存在这个问题,Game可以直接依赖OS内的组件,按照Bazel的语法描述,就是在Game中可以直接使用@ProjectOS//path/to/package:AdvancedStruct.当然,你仍然可以选择重构一

By Edimetia3D

Unix related things

这是一篇2017年左右的记录, 仅用作分享 杂 * 在shell内能干的事,我们都可以比较简单地通过系统调用实现. * `称为反引号,^称为脱字符,常用来表示CTRL * windows的系统调用是不开放的,windows下只能直接使用windows.h里的windows API. * /dev目录下的设备是供用于程序直接使用的,主要由block,char,pipe,socket类型 * 并不是所有设备都能映射为这种形式 * /sys/device/目录称为sysfs,他下面存放了所有设备的信息.(不能直接从/dev获得任何设备信息) * udevadm info --query=all --name="/dev/sda1"可以用于查询/dev下某个设备对应的sysfs路径 权限系统 * 权限系统由两部分组成 * 文件属性:用于标注文件owner,所属组,以及权限的设定(默认只有owner和root可以修改权限设置) *

By Edimetia3D