K-D树算法的一些总结

在2016年ACM ICPC青岛站的比赛中，一道K-D树问题成为了金银牌的分界线，最后由我们队一个强力队友写出了该题，其实那题本来该由我负责的，都怪我学艺不精。

k-d树（k-dimensional树的简称），是一种分割k维数据空间的数据结构。主要应用于多维空间关键数据的搜索（如：范围搜索和最近邻搜索）。

应用背景

　　SIFT算法中做特征点匹配的时候就会利用到k-d树。而特征点匹配实际上就是一个通过距离函数在高维矢量之间进行相似性检索的问题。针对如何快速而准确地找到查询点的近邻，现在提出了很多高维空间索引结构和近似查询的算法，k-d树就是其中一种。

　　索引结构中相似性查询有两种基本的方式：一种是范围查询（range searches），另一种是K近邻查询（K-neighbor searches）。范围查询就是给定查询点和查询距离的阈值，从数据集中找出所有与查询点距离小于阈值的数据；K近邻查询是给定查询点及正整数K，从数据集中找到距离查询点最近的K个数据，当K=1时，就是最近邻查询（nearest neighbor searches）。

　　特征匹配算子大致可以分为两类。一类是线性扫描法，即将数据集中的点与查询点逐一进行距离比较，也就是穷举，缺点很明显，就是没有利用数据集本身蕴含的任何结构信息，搜索效率较低，第二类是建立数据索引，然后再进行快速匹配。因为实际数据一般都会呈现出簇状的聚类形态，通过设计有效的索引结构可以大大加快检索的速度。索引树属于第二类，其基本思想就是对搜索空间进行层次划分。根据划分的空间是否有混叠可以分为Clipping和Overlapping两种。前者划分空间没有重叠，其代表就是k-d树；后者划分空间相互有交叠，其代表为R树。（这里只介绍k-d树）

实例

　　先以一个简单直观的实例来介绍k-d树算法。假设有6个二维数据点{（2,3），（5,4），（9,6），（4,7），（8,1），（7,2）}，数据点位于二维空间内（如图1中黑点所示）。k-d树算法就是要确定图1中这些分割空间的分割线（多维空间即为分割平面，一般为超平面）。下面就要通过一步步展示k-d树是如何确定这些分割线的。



图1  二维数据k-d树空间划分示意图

　　k-d树算法可以分为两大部分，一部分是有关k-d树本身这种数据结构建立的算法，另一部分是在建立的k-d树上如何进行最邻近查找的算法。

再附上HDU2966我的代码：

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=2966

#include <map>
#include <set>
#include <queue>
#include <cmath>
#include <stack>
#include <cstdio>
#include <cctype>
#include <vector>
#include <string>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include <algorithm>

using namespace std;

typedef long long LL;

const int N = 100050;
struct Point {
LL x[3];
}p[N], ori[N];
int split[20], cur, dim;

bool cmp(const Point &a, const Point &b) {
return a.x[cur] < b.x[cur];
}

#define lson l, m - 1, depth + 1
#define rson m + 1, r, depth + 1

template <class T> T sqr(T x) { return x * x; }
const LL inf = 0x7777777777777777ll;

LL dist(const Point &x, const Point &y) {
LL ret = 0;
for (int i = 0; i < dim; i++) {
ret += sqr(x.x[i] - y.x[i]);
}
return ret ? ret : inf;
}

void build(const int &l, const int &r, const int &depth) {
if (l >= r) return;
int m = l + r >> 1;
cur = depth % dim;
nth_element(p + l, p + m, p + r + 1, cmp);
build(lson);
build(rson);
}

LL Find(const Point &x, const int &l, const int &r, const int &depth) {
int cur = depth % dim;
if (l >= r) {
if (l == r) return dist(x, p[l]);
return inf;
}
int m = l + r >> 1;
LL ret = dist(x, p[m]), tmp;
if (x.x[cur] < p[m].x[cur]) {
tmp = Find(x, lson);
if (tmp > sqr(x.x[cur] - p[m].x[cur])) {
tmp = min(tmp, Find(x, rson));
}
} else {
tmp = Find(x, rson);
if (tmp > sqr(x.x[cur] - p[m].x[cur])) {
tmp = min(tmp, Find(x, lson));
}
}
return min(ret, tmp);
}

int main() {
int n, T;
scanf("%d", &T);
while (T--) {
scanf("%d", &n);
dim = 2;
for (int i = 0; i < n; i++) {
for (int j = 0; j < 2; j++) {
scanf("%I64d", &ori[i].x[j]);
}
p[i] = ori[i];
}
build(0, n - 1, 0);
for (int i = 0; i < n; i++) {
printf("%I64d\n", Find(ori[i], 0, n - 1, 0));
}
}
return 0;
}