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表示.
  • 最小的体素是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 )