By RaySaint 2011/10/12

动机

先前写了一篇文章《SIFT算法研究》讲了讲SIFT特征具体是如何检测和描述的，其中也提到了SIFT常见的一个用途就是物体识别，物体识别的过程如下图所示：

如上图(a)，我们先对待识别的物体的图像进行SIFT特征点的检测和特征点的描述，然后得到了SIFT特征点集合。接下来生成物体目标描述要做的就是对特征点集合进行数据组织，形成一种特殊的表示，其作用是为了加速特征点匹配的过程。所谓的特征点匹配本质上是一个通过距离函数（例如欧式距离）在高维矢量之间进行相似性检索的问题，简单来讲就是范围查询或者K近邻查询的问题。

范围查询就是给定查询点和查询距离阈值，从数据集中找出所有与查询点距离小于查询距离阈值的数据；K近邻查询就是给定查询点和正整数K，从数据集中找到距离查询点最近的K个数据，当K=1时，它就是最近邻查询。

如上图(b)我们从输入图像中进行SIFT特征点的检测和特征点的描述后，得到了一个待查询点的集合，接下来就是要找出集合中的每一个待查询点在(a)过程得到的目标物体的特征点集合中进行2近邻查询（即得到最近邻和次近邻），得到一组特征点的匹配对<待查询点，待查询点的最近邻>；得到所有匹配对后，然后通过阈值法（与最近邻的距离要小于一个常数）和比值法(与最近邻的距离比次近邻的距离要小于一个常数)进行提纯，滤去较差的匹配对。得到最终的匹配对集合。最后在计算单应性矩阵时，使用RANSAC算法再进行一次提纯，剔除错误的匹配对。关于RANSAC算法，我还会再写一篇文章讲一讲。David G.Lowe的论文说3个或以上的特征点匹配对可以确认一个正确的识别。（因为单应性矩阵的计算最少得使用4个点，并且可能会有错误匹配的情况存在，所以最好需要多一点的特征点匹配对）

本文的主要目的是讲一下如何创建k-d tree对目标物体的特征点集合进行数据组织和使用k-d tree最近邻搜索来加速特征点匹配。上面已经讲了特征点匹配的问题其实上是一个最近邻（K近邻）搜索的问题。所以为了更好的引出k-d tree，先讲一讲最近邻搜索。

最近邻搜索

先给出一个最近邻的数学形式的定义。给定一个多维空间，把中的一个向量成为一个样本点或数据点。中样本点的有限集合称为样本集。给定样本集E，和一个样本点d，d的最近邻就是任何样本点d’∈E满足None-nearer(E,d,d’)。

None-nearer如下定义：

上面的公式中距离度量是欧式距离，当然也可以是任何其他Lp-norm。

其中di是向量d的第i个分量。

现在再来说最近邻搜索，如何找到一个这样的d’，它离d的距离在E中是最近的。

很容易想到的一个方法就是线性扫描，也称为穷举搜索，依次计算样本集E中每个样本点到d的距离，然后取最小距离的那个点。这个方法又称为朴素最近邻搜索。当样本集E较大时（在物体识别的问题中，可能有数千个甚至数万个SIFT特征点），显然这种策略是非常耗时的。

因为实际数据一般都会呈现簇状的聚类形态，因此我们想到建立数据索引，然后再进行快速匹配。索引树是一种树结构索引方法，其基本思想是对搜索空间进行层次划分。k-d tree是索引树中的一种典型的方法。

k-d tree的简介及表示

k-d tree是英文K-dimension tree的缩写，是对数据点在k维空间中划分的一种数据结构。k-d tree实际上是一种二叉树。每个结点的内容如下：

样本集E由k-d tree的结点的集合表示，每个结点表示一个样本点，dom_elt就是表示该样本点的向量。该样本点根据结点的分割超平面将样本空间分为两个子空间。左子空间中的样本点集合由左子树left表示，右子空间中的样本点集合由右子树right表示。分割超平面是一个通过点dom_elt并且垂直于split所指示的方向轴的平面。举个简单的例子，在二维的情况下，一个样本点可以由二维向量(x,y)表示，其中令x维的序号为0，y维的序号为1。假设一个结点的dom_elt为(7,2) ，split的取值为0，那么分割超面就是x=dom_elt(0)=7，它垂直与x轴且过点(7,2)，如下图所示：

（红线代表分割超平面）

于是其他数据点的x维（第split=0维）如果小于7，则被分配到左子空间；若大于7，则被分配到右子空间。例如，（5,4）被分配到左子空间，（9，6）被分配到右子空间。如下图所示:

从上面的表也可以看出k-d tree本质上是一种二叉树，因此k-d tree的构建是一个逐级展开的递归过程。

其算法的伪代码如下：

1. 算法：createKDTree 构建一棵k-d tree
2. 输入：exm_set 样本集
3. 输出 : Kd, 类型为kd-tree
4. 1. 如果exm_set是空的，则返回空的kd-tree
5. 2.调用分裂结点选择程序（输入是exm_set），返回两个值
6. dom_elt:= exm_set中的一个样本点
7. split := 分裂维的序号
8. 3.exm_set_left = {exm∈exm_set – dom_elt && exm[split] <= dom_elt[split]}
9. exm_set_right = {exm∈exm_set – dom_elt && exm[split] > dom_elt[split]}
10. 4.left = createKDTree(exm_set_left)
11. right = createKDTree(exm_set_right)

现在来解释一下分裂结点选择程序。分裂结点的选择通常有多种方法，最常用的是一种方法是：对于所有的样本点，统计它们在每个维上的方差，挑选出方差中的最大值，对应的维就是split域的值。数据方差最大表明沿该维度数据点分散得比较开，这个方向上进行数据分割可以获得最好的分辨率；然后再将所有样本点按其第split维的值进行排序，位于正中间的那个数据点选为分裂结点的dom_elt域。

下面以一个简单的例子来解释上述k-d tree的构建过程。假设样本集为：{(2,3), (5,4), (9,6), (4,7), (8,1), (7,2)}。构建过程如下：

（1）确定split域，6个数据点在x,y维度上的数据方差分别为39, 28.63。在x轴上方差最大，所以split域值为0（x维的序号为0）

（2）确定分裂节点，根据x维上的值将数据排序，则6个数据点再排序后位于中间的那个数据点为(7,2)，该结点就是分割超平面就是通过(7,2)并垂直于split=0(x)轴的直线x=7

（3）左子空间和右子空间，分割超面x=7将整个空间氛围两部分，x<=7的部分为左子空间，包含3个数据点{(2,3), (5,4), (4,7)}；另一部分为右子空间，包含2个数据点{(9,6), (8,1)}。如下图所示

（4）分别对左子空间中的数据点和右子空间中的数据点重复上面的步骤构建左子树和右子树直到经过划分的子样本集为空。下面的图从左至右从上至下显示了构建这棵二叉树的所有步骤：

k-d tree的最近邻搜索算法

如前所述，在k-d tree树中进行数据的k近邻搜索是特征匹配的重要环节，其目的是检索在k-d tree中与待查询点距离最近的k个数据点。

最近邻搜索是k近邻的特例，也就是1近邻。将1近邻改扩展到k近邻非常容易。下面介绍最简单的k-d tree最近邻搜索算法。

基本的思路很简单：首先通过二叉树搜索（比较待查询节点和分裂节点的分裂维的值，小于等于就进入左子树分支，等于就进入右子树分支直到叶子结点），顺着“搜索路径”很快能找到最近邻的近似点，也就是与待查询点处于同一个子空间的叶子结点；然后再回溯搜索路径，并判断搜索路径上的结点的其他子结点空间中是否可能有距离查询点更近的数据点，如果有可能，则需要跳到其他子结点空间中去搜索（将其他子结点加入到搜索路径）。重复这个过程直到搜索路径为空。下面给出k-d tree最近邻搜索的伪代码：

1. 算法：kdtreeFindNearest /* k-d tree的最近邻搜索 */
2. 输入：Kd /* k-d tree类型*/
3. target /* 待查询数据点 */
4. 输出 : nearest /* 最近邻数据结点 */
5. dist /* 最近邻和查询点的距离 */
6. 1. 如果Kd是空的，则设dist为无穷大返回
7. 2. 向下搜索直到叶子结点
8. pSearch = &Kd
9. while(pSearch != NULL)
10. {
11. pSearch加入到search_path中;
12. if(target[pSearch->split] <= pSearch->dom_elt[pSearch->split]) /* 如果小于就进入左子树 */
13. {
14. pSearch = pSearch->left;
15. }
16. else
17. {
18. pSearch = pSearch->right;
19. }
20. }
21. 取出search_path最后一个赋给nearest
22. dist = Distance(nearest, target);
23. 3. 回溯搜索路径
24. while(search_path不为空)
25. {
26. 取出search_path最后一个结点赋给pBack
27. if(pBack->left为空 && pBack->right为空) /* 如果pBack为叶子结点 */
28. {
29. if( Distance(nearest, target) > Distance(pBack->dom_elt, target) )
30. {
31. nearest = pBack->dom_elt;
32. dist = Distance(pBack->dom_elt, target);
33. }
34. }
35. else
36. {
37. s = pBack->split;
38. if( abs(pBack->dom_elt[s] - target[s]) < dist) /* 如果以target为中心的圆（球或超球），半径为dist的圆与分割超平面相交， 那么就要跳到另一边的子空间去搜索 */
39. {
40. if( Distance(nearest, target) > Distance(pBack->dom_elt, target) )
41. {
42. nearest = pBack->dom_elt;
43. dist = Distance(pBack->dom_elt, target);
44. }
45. if(target[s] <= pBack->dom_elt[s]) /* 如果target位于pBack的左子空间，那么就要跳到右子空间去搜索 */
46. pSearch = pBack->right;
47. else
48. pSearch = pBack->left; /* 如果target位于pBack的右子空间，那么就要跳到左子空间去搜索 */
49. if(pSearch != NULL)
50. pSearch加入到search_path中
51. }
52. }
53. }

OK，现在举一些例子来说明上面的最近邻搜索算法会比较直观。

假设我们的k-d tree就是上面通过样本集{(2,3), (5,4), (9,6), (4,7), (8,1), (7,2)}创建的。将上面的图转化为树形图的样子如下：

我们来查找点(2.1,3.1)，在(7,2)点测试到达(5,4)，在(5,4)点测试到达(2,3)，然后search_path中的结点为<(7,2), (5,4), (2,3)>，从search_path中取出(2,3)作为当前最佳结点nearest, dist为0.141；

然后回溯至(5,4)，以(2.1,3.1)为圆心，以dist=0.141为半径画一个圆，并不和超平面y=4相交，如下图，所以不必跳到结点(5,4)的右子空间去搜索，因为右子空间中不可能有更近样本点了。

于是在回溯至(7,2)，同理，以(2.1,3.1)为圆心，以dist=0.141为半径画一个圆并不和超平面x=7相交，所以也不用跳到结点(7,2)的右子空间去搜索。

至此，search_path为空，结束整个搜索，返回nearest(2,3)作为(2.1,3.1)的最近邻点，最近距离为0.141。

再举一个稍微复杂的例子，我们来查找点(2,4.5)，在(7,2)处测试到达(5,4)，在(5,4)处测试到达(4,7)，然后search_path中的结点为<(7,2), (5,4), (4,7)>，从search_path中取出(4,7)作为当前最佳结点nearest, dist为3.202；

然后回溯至(5,4)，以(2,4.5)为圆心，以dist=3.202为半径画一个圆与超平面y=4相交，如下图，所以需要跳到(5,4)的左子空间去搜索。所以要将(2,3)加入到search_path中，现在search_path中的结点为<(7,2), (2, 3)>；另外，(5,4)与(2,4.5)的距离为3.04 < dist = 3.202，所以将(5,4)赋给nearest，并且dist=3.04。

回溯至(2,3)，(2,3)是叶子节点，直接平判断(2,3)是否离(2,4.5)更近，计算得到距离为1.5，所以nearest更新为(2,3)，dist更新为(1.5)

回溯至(7,2)，同理，以(2,4.5)为圆心，以dist=1.5为半径画一个圆并不和超平面x=7相交, 所以不用跳到结点(7,2)的右子空间去搜索。

至此，search_path为空，结束整个搜索，返回nearest(2,3)作为(2,4.5)的最近邻点，最近距离为1.5。

两次搜索的返回的最近邻点虽然是一样的，但是搜索(2, 4.5)的过程要复杂一些，因为(2, 4.5)更接近超平面。研究表明，当查询点的邻域与分割超平面两侧的空间都产生交集时，回溯的次数大大增加。最坏的情况下搜索N个结点的k维kd-tree所花费的时间为：

后记

到此为止，k-d tree相关的基本知识就说完了。关于k-d tree还有很多扩展。由于大量回溯会导致kd-tree最近邻搜索的性能大大下降，因此研究人员也提出了改进的k-d tree近邻搜索，其中一个比较著名的就是 Best-Bin-First，它通过设置优先级队列和运行超时限定来获取近似的最近邻，有效地减少回溯的次数。这里就不详细讲了，如果想知道可以查询后面的参考资料。

参考资料

1.An intoductory tutorial on kd-trees Andrew W.Moore

2.《图像局部不变特性特征与描述》王永明 王贵锦 编著 国防工业出版社

3.kdtree A simple C library for working with KD-Trees