KD tree

Kd-树 其实是K-dimension tree的缩写，是对数据点在k维空间中划分的一种数据结构。其实，Kd-树是一种平衡二叉树。

举一示例：

假设有六个二维数据点 = {（2,3），（5,4），（9,6），（4,7），（8,1），（7,2）}，数据点位于二维空间中。为了能有效的找到最近邻，Kd-树采用分而治之的思想，即将整个空间划分为几个小部分。六个二维数据点生成的Kd-树的图为：

 对于拥有n个已知点的kD-Tree，其复杂度如下：

1. 构建：O(log²n)
2. 插入：O(log n)
3. 删除：O(log n)
4. 查询：O(n^1-1/k+m) m---每次要搜索的最近点个数

一 Kd-树的构建

Kd-树是一个二叉树，每个节点表示的是一个空间范围。下表表示的是Kd-树中每个节点中主要包含的数据结构。Range域表示的是节点包含的空间范围。Node-data域就是数据集中的某一个n维数据点。分割超面是通过数据点Node-Data并垂直于轴split的平面，分割超面将整个空间分割成两个子空间。令split域的值为i，如果空间Range中某个数据点的第i维数据小于Node-Data[i]，那么，它就属于该节点空间的左子空间，否则就属于右子空间。Left，Right域分别表示由左子空间和右子空间空的数据点构成的Kd-树。

 

域名	数据类型	描述
Node-Data	数据矢量	数据集中某个数据点，是n维矢量
Range	空间矢量	该节点所代表的空间范围
Split	整数	垂直于分割超面的方向轴序号
Left	Kd-tree	由位于该节点分割超面左子空间内所有数据点构成的Kd-树
Right	Kd-tree	由位于该节点分割超面左子空间内所有数据点构成的Kd-树
Parent	Kd-tree	父节点

构建Kd-树的伪码为：

算法：构建Kd-tree

输入：数据点集Data_Set，和其所在的空间。

输出：Kd，类型为Kd-tree

1 if data-set is null ,return 空的Kd-tree

2 调用节点生成程序

 （1）确定split域：对于所有描述子数据（特征矢量），统计他们在每个维度上的数据方差，挑选出方差中最大值，对应的维就是split域的值。数据方差大说明沿该坐标轴方向上数据点分散的比较开。这个方向上，进行数据分割可以获得最好的分辨率。

（2）确定Node-Data域，数据点集Data-Set按照第split维的值排序，位于正中间的那个数据点 被选为Node-Data，Data-Set` =Data-Set\Node-data

3 dataleft = {d 属于Data-Set` & d[:split]<=Node-data[:split]}

   Left-Range ={Range && dataleft}

  dataright = {d 属于Data-Set` & d[:split]>Node-data[:split]}

  Right-Range ={Range && dataright}

4 :left =由（dataleft,LeftRange）建立的Kd-tree

   设置:left的parent域（父节点）为Kd

   ：right =由（dataright,RightRange）建立的Kd-tree

设置：right的parent域为kd。

如上例中，

（1）确定：split 域=x，6个数据点在x,y 维度上的数据方差为39,28.63.在x轴方向上的方差大，所以split域值为x。

（2）确定：Node-Data=（7,2），根据x维上的值将数据排序，6个数据的中值为7，所以node-data域为数据点（7,2）。这样该节点的分割超面就是通过（7,2）并垂直于：split=x轴的直线x=7.

(3)左子空间和右子空间，分割超面x=7将整个空间分为两部分。x<=7 为左子空间，包含节点（2,3），（5,4），（4，7），另一部分为右子空间。包含节点（9,6），（8,1）

这个构建过程是一个递归过程。重复上述过程，直至只包含一个节点。

转自 http://blog.csdn.net/v_july_v/article/details/8203674 

    前两日，在微博上说：“到今天为止，我至少亏欠了3篇文章待写：1、KD树：http://weibo.com/1580904460/z1PosdcKj；2、神经网络：http://weibo.com/1580904460/yBmhfrOGl；3、编程艺术第28章：http://weibo.com/1580904460/z4ZGFiDcY。你看到，blog内的文章与你于别处所见的任何都不同。于是，等啊等，等一台电脑，只好等待..”。得益于朋友田，借了我一台电脑（借他电脑的时候，我连表示感谢，他说“能找到工作全靠你的博客，这点儿小忙还说，不地道”，有的时候，稍许感受到受人信任也是一种压力，愿我不辜负大家对我的信任），于是今天开始Top 10 Algorithms in Data Mining系列第三篇文章，即本文「从K近邻算法谈到KD树、SIFT+BBF算法」的创作。

    一个人坚持自己的兴趣是比较难的，因为太多的人太容易为外界所动了，尤其当你无法从中得到多少实际性的回报时，所幸，我能一直坚持下来。毕达哥拉斯有句名言：万物为数，最近在读「微积分概念的发展史」也感受到了这一点。同时，从算法到数据挖掘、机器学习，再到数学，其中每一个领域任何一个细节都值得探索终生，或许，这就是“终生为学”的意思吧。

    本文第一部分讲K近邻算法，第二部分讲K近邻算法的一种实现--KD树的相关实现及应用，同时，你将看到，K近邻算法同本系列的前两篇文章所讲的决策树分类贝叶斯分类，及支持向量机SVM一样，也是用于解决分类问题的算法，而本数据挖掘十大算法系列也会按照分类，聚类，关联分析，预测回归等问题依次展开阐述。

    OK，行文仓促，本文若有任何漏洞，问题或者错误，欢迎朋友们随时不吝指正，各位的批评也是我继续写下去的动力之一。感谢。

 

第一部分、K近邻算法

1.1、什么是K近邻算法

    何谓K近邻算法，即K-nearest_neighbor，简称KNN算法，单从名字来猜想，可以简单粗暴的认为是：K个最近的邻居，当K=1时，算法便成了最近邻算法，即寻找最近的那个邻居。

    用官方的话来说，所谓K近邻算法，即是给定一个训练数据集，对新的输入实例，在训练数据集中找到与该实例最邻近的K个实例（也就是上面所说的K个邻居），这K个实例的多数属于某个类，就把该输入实例分类到这个类中。根据这个说法，咱们来看下引自维基百科上的一幅图：

    如上图所示，有两类不同的样本数据，分别用蓝色的小正方形和红色的小三角形表示，而图正中间的那个绿色的圆所标示的数据则是待分类的数据。也就是说，现在，我们不知道中间那个绿色的数据是从属于哪一类（蓝色小正方形or红色小三角形），下面，我们就要解决这个问题：给这个绿色的圆分类。
    我们常说，物以类聚，人以群分，判别一个人是一个什么样品质特征的人，常常可以从他/她身边的朋友入手，所谓观其友，而识其人。我们不是要判别上图中那个绿色的圆是属于哪一类数据么，好说，从它的邻居下手。但一次性看多少个邻居呢？从上图中，你还能看到：

• 如果K=3，绿色圆点的最近的3个邻居是2个红色小三角形和1个蓝色小正方形，少数从属于多数，基于统计的方法，判定绿色的这个待分类点属于红色的三角形一类。
• 如果K=5，绿色圆点的最近的5个邻居是2个红色三角形和3个蓝色的正方形，还是少数从属于多数，基于统计的方法，判定绿色的这个待分类点属于蓝色的正方形一类。

    于此我们看到，当无法判定当前待分类点是从属于已知分类中的哪一类时，我们可以依据统计学的理论看它所处的位置特征--看它周围的邻居是哪一类，而把它归为那一类。这就是K近邻算法的核心思想。

 

1.2、近邻的距离度量表示法

    上文第一节，我们看到，K近邻算法的核心在于找到实例点的邻居，这个时候，问题就接踵而至了，如何找到邻居，邻居的判定标准是什么，用什么来度量。这一系列问题便是下面要讲的距离度量表示法。但有的读者可能就有疑问了，我是要找邻居，找相似性，怎么又跟距离扯上关系了？

    这是因为特征空间中两个实例点的距离和反应出两个实例点形似程度。K近邻模型的特征空间一般是n维实属向量空间，使用的距离可以使欧式距离，也是可以是其它距离，下面，就来具体阐述下这些距离度量的表示法。

 

• 1. 欧氏距离，最常见的两点之间或多点之间的距离表示法，又称之为欧几里得度量，它定义于欧几里得空间中，如点 x = (x1,...,xn) 和 y = (y1,...,yn) 之间的距离为：

 

(1)二维平面上两点a(x1,y1)与b(x2,y2)间的欧氏距离：

 

 

(2)三维空间两点a(x1,y1,z1)与b(x2,y2,z2)间的欧氏距离：

 

(3)两个n维向量a(x11,x12,…,x1n)与 b(x21,x22,…,x2n)间的欧氏距离：

 

　　也可以用表示成向量运算的形式：

 

• 2. 曼哈顿距离，我们可以定义曼哈顿距离的正式意义为L1-距离或城市区块距离，也就是在欧几里得空间的固定直角坐标系上两点所形成的线段对轴产生的投影的距离总和。例如在平面上，座标（x1, y1）的点P1与座标（x2, y2）的点P2的曼哈顿距离为：，要注意的是，曼哈顿距离依赖座标系统的转度，而非系统在座标轴上的平移或映射。 

 

 

     通俗来讲，想象你在曼哈顿要从一个十字路口开车到另外一个十字路口，驾驶距离是两点间的直线距离吗？显然不是，除非你能穿越大楼。实际驾驶距离就是这个“曼哈顿距离”。而这也是曼哈顿距离名称的来源， 曼哈顿距离也称为城市街区距离(City Block distance)。

 

 

(1)二维平面两点a(x1,y1)与b(x2,y2)间的曼哈顿距离 

 

 

 

 

(2)两个n维向量a(x11,x12,…,x1n)与 b(x21,x22,…,x2n)间的曼哈顿距离 

                          

 

• 3. 切比雪夫距离，若二个向量或二个点p 、and q，其座标分别为及，则两者之间的切比雪夫距离定义如下：，

    这也等于以下Lp度量的极值：，因此切比雪夫距离也称为L∞度量。

 

    以数学的观点来看，切比雪夫距离是由一致范数（uniform norm）（或称为上确界范数）所衍生的度量，也是超凸度量（injective metric space）的一种。

    在平面几何中，若二点p及q的直角坐标系坐标为及，则切比雪夫距离为：。

 

    玩过国际象棋的朋友或许知道，国王走一步能够移动到相邻的8个方格中的任意一个。那么国王从格子(x1,y1)走到格子(x2,y2)最少需要多少步？。你会发现最少步数总是max( | x2-x1 | , | y2-y1 | ) 步 。有一种类似的一种距离度量方法叫切比雪夫距离。

 

 

(1)二维平面两点a(x1,y1)与b(x2,y2)间的切比雪夫距离 

 

 

 

 

(2)两个n维向量a(x11,x12,…,x1n)与 b(x21,x22,…,x2n)间的切比雪夫距离 　　

 

 

这个公式的另一种等价形式是 

 

 

 

 

• 4. 闵可夫斯基距离(Minkowski Distance)，闵氏距离不是一种距离，而是一组距离的定义。

(1) 闵氏距离的定义       

两个n维变量a(x11,x12,…,x1n)与 b(x21,x22,…,x2n)间的闵可夫斯基距离定义为： 

其中p是一个变参数。

当p=1时，就是曼哈顿距离

当p=2时，就是欧氏距离

当p→∞时，就是切比雪夫距离       

根据变参数的不同，闵氏距离可以表示一类的距离。 

 

 

1.3、K值的选择

    除了上述1.2节如何定义邻居的问题之外，还有一个选择多少个邻居，即K值定义为多大的问题。不要小看了这个K值选择问题，因为它对K近邻算法的结果会产生重大影响。如李航博士的一书「统计学习方法」上所说：

1. 如果选择较小的K值，就相当于用较小的领域中的训练实例进行预测，“学习”近似误差会减小，只有与输入实例较近或相似的训练实例才会对预测结果起作用，与此同时带来的问题是“学习”的估计误差会增大，换句话说，K值的减小就意味着整体模型变得复杂，容易发生过拟合；
2. 如果选择较大的K值，就相当于用较大领域中的训练实例进行预测，其优点是可以减少学习的估计误差，但缺点是学习的近似误差会增大。这时候，与输入实例较远（不相似的）训练实例也会对预测器作用，使预测发生错误，且K值的增大就意味着整体的模型变得简单。
3. K=N，则完全不足取，因为此时无论输入实例是什么，都只是简单的预测它属于在训练实例中最多的累，模型过于简单，忽略了训练实例中大量有用信息。

    在实际应用中，K值一般取一个比较小的数值，例如采用交叉验证法来选择最优的K值。

 

1.4、

 

第二部分、K近邻算法的实现：KD树

2.0、背景

    之前blog内曾经介绍过SIFT特征匹配算法，特征点匹配和数据库查、图像检索本质上是同一个问题，都可以归结为一个通过距离函数在高维矢量之间进行相似性检索的问题，如何快速而准确地找到查询点的近邻，不少人提出了很多高维空间索引结构和近似查询的算法。

 

    一般说来，索引结构中相似性查询有两种基本的方式：

 

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

    同样，针对特征点匹配也有两种方法：

 

• 最容易的办法就是线性扫描，也就是我们常说的穷举搜索，依次计算样本集E中每个样本到输入实例点的距离，然后抽取出计算出来的最小距离的点即为最近邻点。此种办法简单直白，但当样本集或训练集很大时，它的缺点就立马暴露出来了，举个例子，在物体识别的问题中，可能有数千个甚至数万个SIFT特征点，而去一一计算这成千上万的特征点与输入实例点的距离，明显是不足取的。
• 另外一种，就是构建数据索引，因为实际数据一般都会呈现簇状的聚类形态，因此我们想到建立数据索引，然后再进行快速匹配。索引树是一种树结构索引方法，其基本思想是对搜索空间进行层次划分。

    那这种索引树应该是一种什么样的树呢？Rudolf Bayer 和 Ed McCreight 于1972年，在Boeing Research Labs 工作时发明了B 树，后续人们在B树的基础上提出了B+树，B*树，之前blog除了介绍B树系列之外，也介绍了R树，但不论是B树家族还是R树，都不是应对空间层次划分的树。

 

    1975年，来自斯坦福大学的Jon Louis Bentley在ACM杂志上发表的一篇论文：Multidimensional Binary Search Trees Used for Associative Searching 中正式提出和阐述的了k-d树。

2.1、什么是KD树

    首先必须搞清楚的是，k-d树是一种空间划分树，说白了，就是把整个空间划分为特定的几个部分，然后在制定空间的部分内进行相关搜索操作。

    举一个简单直观的实例来介绍k-d树算法。假设有6个二维数据点{（2,3），（5,4），（9,6），（4,7），（8,1），（7,2）}，数据点位于二维空间内，如下图所示。为了能有效的找到最近邻，k-d树采用分而治之的思想，即将整个空间划分为几个小部分，首先，粗黑线将空间一分为二，然后在两个子空间中，细直线又将整个空间划分为四部分，最后虚黑直线将这四部分进一步划分。

    经过对上图所示的空间划分之后，我们可以看出，点(7,2)可以为根结点，(5,4)和(9,6)则为根结点的左右子结点，而(2,3)，(4,7)则为(5,4)的左右孩子，最后，(8,1)为(9,6)的左孩子。如此，便形成了这样一棵k-d树：

 

 

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

 

2.2、KD树的构建

 

 k-d树的数据结构

 

    代码如下所示：

1. /** a node in a k-d tree */  
2. struct kd_node  
3. {  
4.     int ki;                      /**< partition key index *///关键点直方图方差最大向量系列位置  
5.     double kv;                   /**< partition key value *///直方图方差最大向量系列中最中间模值  
6.     int leaf;                    /**< 1 if node is a leaf, 0 otherwise */  
7.     struct feature* features;    /**< features at this node */  
8.     int n;                       /**< number of features */  
9.     struct kd_node* kd_left;     /**< left child */  
10.     struct kd_node* kd_right;    /**< right child */  
11. };  

    以下是构建k-d树的代码：

1. /* 
2. A function to build a k-d tree database from keypoints in an array. 
3.  
4. @param features an array of features 
5. @param n the number of features in features 
6.  
7. @return Returns the root of a kd tree built from features or NULL on error. 
8. */  
9. struct kd_node* kdtree_build( struct feature* features, int n )  
10. {  
11.     struct kd_node* kd_root;  
12.   
13.     if( ! features  ||  n <= 0 )  
14.     {  
15.         fprintf( stderr, "Warning: kdtree_build(): no features, %s, line %d\n",  
16.                 __FILE__, __LINE__ );  
17.         return NULL;  
18.     }  
19.   
20.     //初始化  
21.     kd_root = kd_node_init( features, n );  //n--number of features,initinalize root of tree.  
22.     expand_kd_node_subtree( kd_root );  //kd tree expand  
23.   
24.     return kd_root;  
25. }  

    上面的涉及初始化操作的两个函数kd_node_init，及expand_kd_node_subtree代码分别如下所示：

 

1. /* 
2. Initializes a kd tree node with a set of features.  The node is not 
3. expanded, and no ordering is imposed on the features. 
4.  
5. @param features an array of image features 
6. @param n number of features 
7.  
8. @return Returns an unexpanded kd-tree node. 
9. */  
10. static struct kd_node* kd_node_init( struct feature* features, int n )  
11. {                                     //n--number of features  
12.     struct kd_node* kd_node;  
13.   
14.     kd_node = (struct kd_node*)(malloc( sizeof( struct kd_node ) ));  
15.     memset( kd_node, 0, sizeof( struct kd_node ) ); //0填充  
16.     kd_node->ki = -1; //???????  
17.     kd_node->features = features;  
18.     kd_node->n = n;  
19.   
20.     return kd_node;  
21. }  

1. /* 
2. Recursively expands a specified kd tree node into a tree whose leaves 
3. contain one entry each. 
4.  
5. @param kd_node an unexpanded node in a kd tree 
6. */  
7. static void expand_kd_node_subtree( struct kd_node* kd_node )  
8. {  
9.     /* base case: leaf node */  
10.     if( kd_node->n == 1  ||  kd_node->n == 0 )  
11.     {   //叶节点               //伪叶节点  
12.         kd_node->leaf = 1;  
13.         return;  
14.     }  
15.   
16.     assign_part_key( kd_node ); //get ki,kv  
17.     partition_features( kd_node ); //creat left and right children,特征点ki位置左树比右树模值小,kv作为分界模值  
18.                                  //kd_node中关键点已经排序  
19.     if( kd_node->kd_left )  
20.         expand_kd_node_subtree( kd_node->kd_left );  
21.     if( kd_node->kd_right )  
22.         expand_kd_node_subtree( kd_node->kd_right );  
23. }  

    构建完kd树之后，如今进行最近邻搜索呢？从下面的动态gif图中，你是否能看出些许端倪呢？

2.3、KD树的最近邻搜索算法

    在k-d树中进行数据的查找也是特征匹配的重要环节，其目的是检索在k-d树中与查询点距离最近的数据点。这里先以两个简单的实例来描述最邻近查找的基本思路。

　　星号表示要查询的点（2.1,3.1）。通过二叉搜索，顺着搜索路径很快就能找到最邻近的近似点，也就是叶子节点（2,3）。而找到的叶子节点并不一定就是最邻近的，最邻近肯定距离查询点更近，应该位于以查询点为圆心且通过叶子节点的圆域内。为了找到真正的最近邻，还需要进行'回溯'操作：算法沿搜索路径反向查找是否有距离查询点更近的数据点。此例中先从（7,2）点开始进行二叉查找，然后到达（5,4），最后到达（2,3），此时搜索路径中的节点为<（7,2），（5,4），（2,3）>，首先以（2,3）作为当前最近邻点，计算其到查询点（2.1,3.1）的距离为0.1414，然后回溯到其父节点（5,4），并判断在该父节点的其他子节点空间中是否有距离查询点更近的数据点。以（2.1,3.1）为圆心，以0.1414为半径画圆，如图4所示。发现该圆并不和超平面y = 4交割，因此不用进入（5,4）节点右子空间中去搜索。

　　再回溯到（7,2），以（2.1,3.1）为圆心，以0.1414为半径的圆更不会与x = 7超平面交割，因此不用进入（7,2）右子空间进行查找。至此，搜索路径中的节点已经全部回溯完，结束整个搜索，返回最近邻点（2,3），最近距离为0.1414。

　　一个复杂点了例子如查找点为（2，4.5）。同样先进行二叉查找，先从（7,2）查找到（5,4）节点，在进行查找时是由y = 4为分割超平面的，由于查找点为y值为4.5，因此进入右子空间查找到（4,7），形成搜索路径<（7,2），（5,4），（4,7）>，取（4,7）为当前最近邻点，计算其与目标查找点的距离为3.202。然后回溯到（5,4），计算其与查找点之间的距离为3.041。以（2，4.5）为圆心，以3.041为半径作圆，如图5所示。可见该圆和y = 4超平面交割，所以需要进入（5,4）左子空间进行查找。此时需将（2,3）节点加入搜索路径中得<（7,2），（2,3）>。回溯至（2,3）叶子节点，（2,3）距离（2,4.5）比（5,4）要近，所以最近邻点更新为（2，3），最近距离更新为1.5。回溯至（7,2），以（2,4.5）为圆心1.5为半径作圆，并不和x = 7分割超平面交割，如图6所示。至此，搜索路径回溯完。返回最近邻点（2,3），最近距离1.5。

    k-d树查询算法的伪代码如下所示：

 

1. 算法：k-d树最邻近查找  
2. 输入：Kd，    //k-d tree类型  
3.   
4.      target  //查询数据点  
5.   
6. 输出：nearest， //最邻近数据点  
7.   
8.      dist      //最邻近数据点和查询点间的距离  
9.   
10. 1. If Kd为NULL，则设dist为infinite并返回  
11. 2. //进行二叉查找，生成搜索路径  
12.   
13.    Kd_point = &Kd；                   //Kd-point中保存k-d tree根节点地址  
14.   
15.    nearest = Kd_point -> Node-data；  //初始化最近邻点  
16.   
17.    while（Kd_point）  
18.   
19.    　　push（Kd_point）到search_path中； //search_path是一个堆栈结构，存储着搜索路径节点指针  
20.   
21.  /*** If Dist（nearest，target） > Dist（Kd_point -> Node-data，target） 
22.  
23.    　　　　nearest  = Kd_point -> Node-data；    //更新最近邻点 
24.  
25.    　　　　Max_dist = Dist(Kd_point，target）；  //更新最近邻点与查询点间的距离  ***/  
26.   
27.    　　s = Kd_point -> split；                       //确定待分割的方向  
28.   
29.    　　If target[s] <= Kd_point -> Node-data[s]     //进行二叉查找  
30.   
31.    　　　　Kd_point = Kd_point -> left；  
32.   
33.    　　else  
34.   
35.    　　　　Kd_point = Kd_point ->right；  
36.   
37.    nearest = search_path中最后一个叶子节点； //注意：二叉搜索时不比计算选择搜索路径中的最邻近点，这部分已被注释  
38.   
39.    Max_dist = Dist（nearest，target）；    //直接取最后叶子节点作为回溯前的初始最近邻点  
40.   
41. 3. //回溯查找  
42.   
43.    while（search_path != NULL）  
44.   
45.    　　back_point = 从search_path取出一个节点指针；   //从search_path堆栈弹栈  
46.   
47.    　　s = back_point -> split；                   //确定分割方向  
48.   
49.    　　If Dist（target[s]，back_point -> Node-data[s]） < Max_dist   //判断还需进入的子空间  
50.   
51.    　　　　If target[s] <= back_point -> Node-data[s]  
52.   
53.    　　　　　　Kd_point = back_point -> right；  //如果target位于左子空间，就应进入右子空间  
54.   
55.    　　　　else  
56.   
57.    　　　　　　Kd_point = back_point -> left;    //如果target位于右子空间，就应进入左子空间  
58.   
59.    　　　　将Kd_point压入search_path堆栈；  
60.   
61.    　　If Dist（nearest，target） > Dist（Kd_Point -> Node-data，target）  
62.   
63.    　　　　nearest  = Kd_point -> Node-data；                 //更新最近邻点  
64.   
65.    　　　　Min_dist = Dist（Kd_point -> Node-data,target）；  //更新最近邻点与查询点间的距离  

    上述两次实例表明，当查询点的邻域与分割超平面两侧空间交割时，需要查找另一侧子空间，导致检索过程复杂，效率下降。研究表明N个节点的K维k-d树搜索过程时间复杂度为：t_worst=O（kN^1-1/k）。.....

    同时，以上为了介绍方便，讨论的是二维情形。像实际的应用中，如SIFT特征矢量128维，SURF特征矢量64维，维度都比较大，直接利用k-d树快速检索（维数不超过20）的性能急剧下降。假设数据集的维数为D，一般来说要求数据的规模N满足N?2^D，才能达到高效的搜索。所以这就引出了一系列对k-d树算法的改进。

2.4、kd树近邻搜索算法的改进：BBF算法

    咱们顺着上一节的思路，参考统计学习方法一书上的内容，再来总结下kd树的最近邻搜索算法：

输入：以构造的kd树，目标点x；
输出：x 的最近邻
算法步骤如下：

1. 在kd树种找出包含目标点x的叶结点：从根结点出发，递归地向下搜索kd树。若目标点x当前维的坐标小于切分点的坐标，则移动到左子结点，否则移动到右子结点，直到子结点为叶结点为止。
2. 以此叶结点为“当前最近点”。
3. 递归的向上回溯，在每个结点进行以下操作：
   （a）如果该结点保存的实例点比当前最近点距离目标点更近，则更新“当前最近点”，也就是说以该实例点为“当前最近点”。
   （b）当前最近点一定存在于该结点一个子结点对应的区域，检查子结点的父结点的另一子结点对应的区域是否有更近的点。具体做法是，检查另一子结点对应的区域是否以目标点位球心，以目标点与“当前最近点”间的距离为半径的圆或超球体相交：
   如果相交，可能在另一个子结点对应的区域内存在距目标点更近的点，移动到另一个子结点，接着，继续递归地进行最近邻搜索；
   如果不相交，向上回溯。
4. 当回退到根结点时，搜索结束，最后的“当前最近点”即为x 的最近邻点。

    如果实例点是随机分布的，那么kd树搜索的平均计算复杂度是O（NlogN），这里的N是训练实例树。所以说，kd树更适用于训练实例数远大于空间维数时的k近邻搜索，当空间维数接近训练实例数时，它的效率会迅速下降，一降降到“解放前”：线性扫描的速度。

    也正因为上述k最近邻搜索算法的第4个步骤中的所述：“回退到根结点时，搜索结束”，每个最近邻点的查询比较完成过程最终都要回退到根结点而结束，而导致了许多不必要回溯访问和比较到的结点，这些多余的损耗在高维度数据查找的时候，搜索效率将变得相当之地下，那有什么办法可以改进这个原始的kd树最近邻搜索算法呢？

    从上述标准的kd树查询过程可以看出其搜索过程中的“回溯”是由“查询路径”决定的，并没有考虑查询路径上一些数据点本身的一些性质。一个简单的改进思路就是将“查询路径”上的结点进行排序，如按各自分割超平面（也称bin）与查询点的距离排序，也就是说，回溯检查总是从西优先级最高（Best Bin）的树结点开始。

    如此，就引出了本节要讨论的kd树最近邻搜索算法的改进：BBF（Best-Bin-First）查询算法，它能确保优先检索包含最近邻点可能性较高的空间，此外，BBF机制还设置了一个运行超时限定。采用了BBF查询机制后，kd树便可以有效的扩展到高维数据集上。

    伪代码如下图所示：

 

    还是以上面的查询（2,4.5）为例，搜索的算法流程为：

 

1. 将（7,2）压人优先队列中；
2. 提取优先队列中的（7,2），由于（2,4.5）位于（7,2）分割超平面的左侧，所以检索其左子结点。同时，将其右子结点压人优先队列中，此时优先队列为{（9,6）}，最佳点为（7,2）；然后一直检索到叶子结点（4,7），此时优先队列为{（2,3），（9,6）}，“最佳点”则为（5,4）；
3. 提取优先级最高的结点（2,3），重复步骤2，直到优先队列为空。

    如你在下图所见到的那样（话说，用鼠标在图片上写字着实不好写）：

 

....    

第三部分、KD树的应用：SIFT+kd_BBF搜索算法

3.1、SIFT特征匹配算法    

    之前本blog内阐述过图像特征匹配SIFT算法，写过五篇文章，这五篇文章分别为：

• 九、图像特征提取与匹配之SIFT算法      (sift算法系列五篇文章)
• 九（续）、sift算法的编译与实现
• 九（再续）、教你一步一步用c语言实现sift算法、上
• 九（再续）、教你一步一步用c语言实现sift算法、下
• 九（三续）：SIFT算法的应用--目标识别之Bag-of-words模型

    不熟悉SIFT算法相关概念的可以看上述几篇文章，这里不再做赘述。在sift算法中，给定两幅图片图片，若要做特征匹配，一般会先提取出图片中的下列相关属性作为特征点（注，以下代码分析基于Rob Hess维护的sift库）：

 

1. /** 
2. Structure to represent an affine invariant image feature.  The fields 
3. x, y, a, b, c represent the affine region around the feature: 
4.  
5. a(x-u)(x-u) + 2b(x-u)(y-v) + c(y-v)(y-v) = 1 
6. */  
7. struct feature  
8. {  
9.     double x;                      /**< x coord */  
10.     double y;                      /**< y coord */  
11.     double a;                      /**< Oxford-type affine region parameter */  
12.     double b;                      /**< Oxford-type affine region parameter */  
13.     double c;                      /**< Oxford-type affine region parameter */  
14.     double scl;                    /**< scale of a Lowe-style feature */  
15.     double ori;                    /**< orientation of a Lowe-style feature */  
16.     int d;                         /**< descriptor length */  
17.     double descr[FEATURE_MAX_D];   /**< descriptor */  
18.     int type;                      /**< feature type, OXFD or LOWE */  
19.     int category;                  /**< all-purpose feature category */  
20.     struct feature* fwd_match;     /**< matching feature from forward image */  
21.     struct feature* bck_match;     /**< matching feature from backmward image */  
22.     struct feature* mdl_match;     /**< matching feature from model */  
23.     CvPoint2D64f img_pt;           /**< location in image */  
24.     CvPoint2D64f mdl_pt;           /**< location in model */  
25.     void* feature_data;            /**< user-definable data */  
26.     char dense;                     /*表征特征点所处稠密程度*/  
27. };  

    而后在sift.h文件中定义两个关键函数，这里，我们把它们称之为函数一，和函数二，如下所示：

 

函数一的声明：

 

1. extern int sift_features( IplImage* img, struct feature** feat );  

 

函数一的实现：

 

1. int sift_features( IplImage* img, struct feature** feat )  
2. {  
3.     return _sift_features( img, feat, SIFT_INTVLS, SIFT_SIGMA, SIFT_CONTR_THR,  
4.                             SIFT_CURV_THR, SIFT_IMG_DBL, SIFT_DESCR_WIDTH,  
5.                             SIFT_DESCR_HIST_BINS );  
6. }  

从上述函数一的实现中，我们可以看到，它内部实际上调用的是这个函数：_sift_features(..)，也就是下面马上要分析的函数二。

 

函数二的声明：

 

1. extern int _sift_features( IplImage* img, struct feature** feat, int intvls,  
2.                           double sigma, double contr_thr, int curv_thr,  
3.                           int img_dbl, int descr_width, int descr_hist_bins );  

 

函数二的实现：

1. /** 
2. Finds SIFT features in an image using user-specified parameter values.  All 
3. detected features are stored in the array pointed to by \a feat. 
4.  
5. @param img the image in which to detect features 
6. @param fea a pointer to an array in which to store detected features 
7. @param intvls the number of intervals sampled per octave of scale space 
8. @param sigma the amount of Gaussian smoothing applied to each image level 
9.     before building the scale space representation for an octave 
10. @param cont_thr a threshold on the value of the scale space function 
11.     \f$\left|D(\hat{x})\right|\f$, where \f$\hat{x}\f$ is a vector specifying 
12.     feature location and scale, used to reject unstable features;  assumes 
13.     pixel values in the range [0, 1] 
14. @param curv_thr threshold on a feature's ratio of principle curvatures 
15.     used to reject features that are too edge-like 
16. @param img_dbl should be 1 if image doubling prior to scale space 
17.     construction is desired or 0 if not 
18. @param descr_width the width, \f$n\f$, of the \f$n \times n\f$ array of 
19.     orientation histograms used to compute a feature's descriptor 
20. @param descr_hist_bins the number of orientations in each of the 
21.     histograms in the array used to compute a feature's descriptor 
22.  
23. @return Returns the number of keypoints stored in \a feat or -1 on failure 
24. @see sift_keypoints() 
25. */  
26. int _sift_features( IplImage* img, struct feature** feat, int intvls,  
27.                    double sigma, double contr_thr, int curv_thr,  
28.                    int img_dbl, int descr_width, int descr_hist_bins )  
29. {  
30.     IplImage* init_img;  
31.     IplImage*** gauss_pyr, *** dog_pyr;  
32.     CvMemStorage* storage;  
33.     CvSeq* features;  
34.     int octvs, i, n = 0,n0 = 0,n1 = 0,n2 = 0,n3 = 0,n4 = 0;  
35.     int start;  
36.   
37.     /* check arguments */  
38.     if( ! img )  
39.         fatal_error( "NULL pointer error, %s, line %d",  __FILE__, __LINE__ );  
40.   
41.     if( ! feat )  
42.         fatal_error( "NULL pointer error, %s, line %d",  __FILE__, __LINE__ );  
43.   
44.     /* build scale space pyramid; smallest dimension of top level is ~4 pixels */  
45.     start=GetTickCount();  
46.     init_img = create_init_img( img, img_dbl, sigma );  
47.     octvs = log( (float)(MIN( init_img->width, init_img->height )) ) / log((float)(2.0)) -5;  
48.     gauss_pyr = build_gauss_pyr( init_img, octvs, intvls, sigma );  
49.     dog_pyr = build_dog_pyr( gauss_pyr, octvs, intvls );  
50.     fprintf( stderr, " creat the pyramid use %d\n",GetTickCount()-start);  
51.   
52.     storage = cvCreateMemStorage( 0 );    //创建存储内存，0为默认64k  
53.     start=GetTickCount();  
54.     features = scale_space_extrema( dog_pyr, octvs, intvls, contr_thr,  
55.         curv_thr, storage );  //在DOG空间寻找极值点，确定关键点位置  
56.     fprintf( stderr, " find the extrum points in DOG use %d\n",GetTickCount()-start);  
57.   
58.     calc_feature_scales( features, sigma, intvls ); //计算关键点的尺度  
59.   
60.     if( img_dbl )  
61.         adjust_for_img_dbl( features );  //如果原始空间图扩大，特征点坐标就缩小  
62.     start=GetTickCount();  
63.     calc_feature_oris( features, gauss_pyr );  //在gaussian空间计算关键点的主方向和幅值  
64.     fprintf( stderr, " get the main oritation use %d\n",GetTickCount()-start);  
65.   
66.     start=GetTickCount();  
67.     compute_descriptors( features, gauss_pyr, descr_width, descr_hist_bins ); //建立关键点描述器  
68.     fprintf( stderr, " compute the descriptors use %d\n",GetTickCount()-start);  
69.   
70.     /* sort features by decreasing scale and move from CvSeq to array */  
71.     //start=GetTickCount();  
72.     //cvSeqSort( features, (CvCmpFunc)feature_cmp, NULL ); //?????  
73.     n = features->total;  
74.     *feat = (feature*)(calloc( n, sizeof(struct feature) ));  
75.     *feat = (feature*)(cvCvtSeqToArray( features, *feat, CV_WHOLE_SEQ )); //整条链表放在feat指向的内存  
76.   
77.     for( i = 0; i < n; i++ )  
78.     {  
79.         free( (*feat)[i].feature_data );  
80.         (*feat)[i].feature_data = NULL;  //释放ddata(r,c,octv,intvl,xi,scl_octv)  
81.         if((*feat)[i].dense == 4) ++n4;  
82.         else if((*feat)[i].dense == 3) ++n3;  
83.         else if((*feat)[i].dense == 2) ++n2;  
84.         else if((*feat)[i].dense == 1) ++n1;  
85.         else                         ++n0;  
86.     }  
87.   
88.     //fprintf( stderr, " move features from sequce to array use %d\n",GetTickCount()-start);  
89.     //start=GetTickCount();  
90.     fprintf( stderr, "In the total feature points the extent4 points is %d\n",n4);  
91.     fprintf( stderr, "In the total feature points the extent3 points is %d\n",n3);  
92.     fprintf( stderr, "In the total feature points the extent2 points is %d\n",n2);  
93.     fprintf( stderr, "In the total feature points the extent1 points is %d\n",n1);  
94.     fprintf( stderr, "In the total feature points the extent0 points is %d\n",n0);  
95.     cvReleaseMemStorage( &storage );  
96.     cvReleaseImage( &init_img );  
97.     release_pyr( &gauss_pyr, octvs, intvls + 3 );  
98.     release_pyr( &dog_pyr, octvs, intvls + 2 );  
99.     //fprintf( stderr, " free the pyramid use %d\n",GetTickCount()-start);  
100.     return n;  
101. }  

    说明：上面的函数二，包含了SIFT算法中几乎所有函数，是SIFT算法的核心。本文不打算一一分析上面所有函数，只会抽取其中涉及到BBF查询机制相关的函数。

 

3.2、KD树近邻搜索改进之BBF算法

    原理在上文第二部分已经阐述明白，代码如下：

 

1. //KD树近邻搜索改进之BBF算法  
2. int kdtree_bbf_knn( struct kd_node* kd_root, struct feature* feat, int k,  
3.                     struct feature*** nbrs, int max_nn_chks )                  //2  
4. {                                               //200  
5.     struct kd_node* expl;  
6.     struct min_pq* min_pq;  
7.     struct feature* tree_feat, ** _nbrs;  
8.     struct bbf_data* bbf_data;  
9.     int i, t = 0, n = 0;  
10.   
11.     if( ! nbrs  ||  ! feat  ||  ! kd_root )  
12.     {  
13.         fprintf( stderr, "Warning: NULL pointer error, %s, line %d\n",  
14.                 __FILE__, __LINE__ );  
15.         return -1;  
16.     }  
17.   
18.     _nbrs = (feature**)(calloc( k, sizeof( struct feature* ) )); //2  
19.     min_pq = minpq_init();   
20.     minpq_insert( min_pq, kd_root, 0 ); //把根节点加入搜索序列中  
21.   
22.     //队列有东西就继续搜，同时控制在t<200步内  
23.     while( min_pq->n > 0  &&  t < max_nn_chks )  
24.     {                   
25.         //刚进来时,从kd树根节点搜索,exp1是根节点   
26.         //后进来时,exp1是min_pq差值最小的未搜索节点入口  
27.         //同时按min_pq中父,子顺序依次检验,保证父节点的差值比子节点小.这样减少返回搜索时间  
28.         expl = (struct kd_node*)minpq_extract_min( min_pq );  
29.         if( ! expl )                                          
30.         {                                                     
31.             fprintf( stderr, "Warning: PQ unexpectedly empty, %s line %d\n",  
32.                     __FILE__, __LINE__ );                           
33.             goto fail;  
34.         }  
35.   
36.         //从根节点(或差值最小节点)搜索,根据目标点与节点模值的差值(小)  
37.         //确定在kd树的搜索路径,同时存储各个节点另一入口地址\同级搜索路径差值.  
38.         //存储时比较父节点的差值,如果子节点差值比父节点差值小,交换两者存储位置,  
39.         //使未搜索节点差值小的存储在min_pq的前面,减小返回搜索的时间.  
40.         expl = explore_to_leaf( expl, feat, min_pq );   
41.         if( ! expl )                                    
42.         {                                               
43.             fprintf( stderr, "Warning: PQ unexpectedly empty, %s line %d\n",  
44.                     __FILE__, __LINE__ );                     
45.             goto fail;                                    
46.         }                                               
47.   
48.         for( i = 0; i < expl->n; i++ )  
49.         {   
50.             //使用exp1->n原因:如果是叶节点,exp1->n=1,如果是伪叶节点,exp1->n=0.  
51.             tree_feat = &expl->features[i];  
52.             bbf_data = (struct bbf_data*)(malloc( sizeof( struct bbf_data ) ));  
53.             if( ! bbf_data )  
54.             {  
55.                 fprintf( stderr, "Warning: unable to allocate memory,"  
56.                     " %s line %d\n", __FILE__, __LINE__ );  
57.                 goto fail;  
58.             }  
59.             bbf_data->old_data = tree_feat->feature_data;  
60.             bbf_data->d = descr_dist_sq(feat, tree_feat); //计算两个关键点描述器差平方和  
61.             tree_feat->feature_data = bbf_data;  
62.   
63.             //取前k个  
64.             n += insert_into_nbr_array( tree_feat, _nbrs, n, k );//  
65.         }                                                  //2  
66.         t++;  
67.     }  
68.   
69.     minpq_release( &min_pq );  
70.     for( i = 0; i < n; i++ ) //bbf_data为何搞个old_data?  
71.     {  
72.         bbf_data = (struct bbf_data*)(_nbrs[i]->feature_data);  
73.         _nbrs[i]->feature_data = bbf_data->old_data;  
74.         free( bbf_data );  
75.     }  
76.     *nbrs = _nbrs;  
77.     return n;  
78.   
79. fail:  
80.     minpq_release( &min_pq );  
81.     for( i = 0; i < n; i++ )  
82.     {  
83.         bbf_data = (struct bbf_data*)(_nbrs[i]->feature_data);  
84.         _nbrs[i]->feature_data = bbf_data->old_data;  
85.         free( bbf_data );  
86.     }  
87.     free( _nbrs );  
88.     *nbrs = NULL;  
89.     return -1;  
90. }  

 

3.3、SIFT+BBF算法匹配效果

    之前试了下sift + KD + BBF算法，用两幅不同的图片做了下匹配（当然，运行结果显示是不匹配的），效果还不错：http://weibo.com/1580904460/yDmzAEwcV#1348475194313

    “编译的过程中，直接用的VS2010 + opencv（并没下gsl）。2012.09.24”。....