前言前两日,在微博上说:“到今天为止,我至少亏欠了3篇文章待写:1、KD树;2、神经网络;3、编程艺术第28章。你看到,blog内的文章与你于别处所见的任何都不同。于是,等啊等,等一台电脑,只好等待..”。得益于田,借了我一台电脑(借他电脑的时候,我连表示感谢,他说“能找到工作全靠你的博客,这点儿小忙还说,不地道”,有的时候,稍许感受到受人信任也是一种压力,愿我不辜负大家对我的信任),于是今天开始Top 10 Algorithms in Data Mining系列第三篇文章,即本文「从K近邻算法谈到KD树、SIFT+BBF算法」的创作。 一个人坚持自己的兴趣是比较难的,因为太多的人太容易为外界所动了,而尤其当你无法从中得到多少实际性的回报时,所幸,我能一直坚持下来。毕达哥拉斯学派有句名言:“万物皆数”,最近读完「微积分概念发展史」后也感受到了这一点。同时,从算法到数据挖掘、机器学习,再到数学,其中每一个领域任何一个细节都值得探索终生,或许,这就是“终生为学”的意思。 本文各部分内容分布如下:
同时,你将看到,K近邻算法同本系列的前两篇文章所讲的决策树分类贝叶斯分类,及支持向量机SVM一样,也是用于解决分类问题的算法, 而本数据挖掘十大算法系列也会按照分类,聚类,关联分析,预测回归等问题依次展开阐述。 OK,行文仓促,本文若有任何漏洞,问题或者错误,欢迎朋友们随时不吝指正,各位的批评也是我继续写下去的动力之一。感谢。 第一部分、K近邻算法1.1、什么是K近邻算法何谓K近邻算法,即K-Nearest Neighbor algorithm,简称KNN算法,单从名字来猜想,可以简单粗暴的认为是:K个最近的邻居,当K=1时,算法便成了最近邻算法,即寻找最近的那个邻居。为何要找邻居?打个比方来说,假设你来到一个陌生的村庄,现在你要找到与你有着相似特征的人群融入他们,所谓入伙。 用官方的话来说,所谓K近邻算法,即是给定一个训练数据集,对新的输入实例,在训练数据集中找到与该实例最邻近的K个实例(也就是上面所说的K个邻居),这K个实例的多数属于某个类,就把该输入实例分类到这个类中。根据这个说法,咱们来看下引自维基百科上的一幅图: 如上图所示,有两类不同的样本数据,分别用蓝色的小正方形和红色的小三角形表示,而图正中间的那个绿色的圆所标示的数据则是待分类的数据。也就是说,现在,我们不知道中间那个绿色的数据是从属于哪一类(蓝色小正方形or红色小三角形),下面,我们就要解决这个问题:给这个绿色的圆分类。
于此我们看到,当无法判定当前待分类点是从属于已知分类中的哪一类时,我们可以依据统计学的理论看它所处的位置特征,衡量它周围邻居的权重,而把它归为(或分配)到权重更大的那一类。这就是K近邻算法的核心思想。 1.2、近邻的距离度量表示法上文第一节,我们看到,K近邻算法的核心在于找到实例点的邻居,这个时候,问题就接踵而至了,如何找到邻居,邻居的判定标准是什么,用什么来度量。这一系列问题便是下面要讲的距离度量表示法。但有的读者可能就有疑问了,我是要找邻居,找相似性,怎么又跟距离扯上关系了? 这是因为特征空间中两个实例点的距离可以反应出两个实例点之间的相似性程度。K近邻模型的特征空间一般是n维实数向量空间,使用的距离可以使欧式距离,也是可以是其它距离,既然扯到了距离,下面就来具体阐述下都有哪些距离度量的表示法,权当扩展。
其上,二维平面上两点欧式距离,代码可以如下编写: //unixfy:计算欧氏距离double euclideanDistance(const vector
通俗来讲,想象你在曼哈顿要从一个十字路口开车到另外一个十字路口,驾驶距离是两点间的直线距离吗?显然不是,除非你能穿越大楼。而实际驾驶距离就是这个“曼哈顿距离”,此即曼哈顿距离名称的来源, 同时,曼哈顿距离也称为城市街区距离(City Block distance)。
这也等于以下Lp度量的极值:,因此切比雪夫距离也称为L∞度量。 以数学的观点来看,切比雪夫距离是由一致范数(uniform norm)(或称为上确界范数)所衍生的度量,也是超凸度量(injective metric space)的一种。 在平面几何中,若二点p及q的直角坐标系坐标为及,则切比雪夫距离为:。 玩过国际象棋的朋友或许知道,国王走一步能够移动到相邻的8个方格中的任意一个。那么国王从格子(x1,y1)走到格子(x2,y2)最少需要多少步?。你会发现最少步数总是max( | x2-x1 | , | y2-y1 | ) 步 。有一种类似的一种距离度量方法叫切比雪夫距离。
//动态规划: //f[i,j]表示s[0...i]与t[0...j]的最小编辑距离。 f[i,j] = min { f[i-1,j]+1, f[i,j-1]+1, f[i-1,j-1]+(s[i]==t[j]?0:1) } //分别表示:添加1个,删除1个,替换1个(相同就不用替换)。
相关系数 ( Correlation coefficient )的定义是: (其中,E为数学期望或均值,D为方差,D开根号为标准差,E{ [X-E(X)] [Y-E(Y)]}称为随机变量X与Y的协方差,记为Cov(X,Y),即Cov(X,Y) = E{ [X-E(X)] [Y-E(Y)]},而两个变量之间的协方差和标准差的商则称为随机变量X与Y的相关系数,记为) 相关系数衡量随机变量X与Y相关程度的一种方法,相关系数的取值范围是[-1,1]。相关系数的绝对值越大,则表明X与Y相关度越高。当X与Y线性相关时,相关系数取值为1(正线性相关)或-1(负线性相关)。 具体的,如果有两个变量:X、Y,最终计算出的相关系数的含义可以有如下理解:
相关距离的定义是:
以上方程定义了总体相关系数, 一般表示成希腊字母ρ(rho)。基于样本对协方差和方差进行估计,可以得到样本标准差, 一般表示成r: 一种等价表达式的是表示成标准分的均值。基于(Xi, Yi)的样本点,样本皮尔逊系数是
(2)皮尔逊相关系数的适用范围 (3)如何理解皮尔逊相关系数
简单说来,各种“距离”的应用场景简单概括为,空间:欧氏距离,路径:曼哈顿距离,国际象棋国王:切比雪夫距离,以上三种的统一形式:闵可夫斯基距离,加权:标准化欧氏距离,排除量纲和依存:马氏距离,向量差距:夹角余弦,编码差别:汉明距离,集合近似度:杰卡德类似系数与距离,相关:相关系数与相关距离。 1.3、K值的选择除了上述1.2节如何定义邻居的问题之外,还有一个选择多少个邻居,即K值定义为多大的问题。不要小看了这个K值选择问题,因为它对K近邻算法的结果会产生重大影响。如李航博士的一书「统计学习方法」上所说:
在实际应用中,K值一般取一个比较小的数值,例如采用交叉验证法(简单来说,就是一部分样本做训练集,一部分做测试集)来选择最优的K值。 第二部分、K近邻算法的实现:KD树2.0、背景之前blog内曾经介绍过SIFT特征匹配算法,特征点匹配和数据库查、图像检索本质上是同一个问题,都可以归结为一个通过距离函数在高维矢量之间进行相似性检索的问题,如何快速而准确地找到查询点的近邻,不少人提出了很多高维空间索引结构和近似查询的算法。 一般说来,索引结构中相似性查询有两种基本的方式:
同样,针对特征点匹配也有两种方法:
而关于R树本blog内之前已有介绍(同时,关于基于R树的最近邻查找,还可以看下这篇文章:http://blog.sina.com.cn/s/blog_72e1c7550101dsc3.html),本文着重介绍k-d树。 1975年,来自斯坦福大学的Jon Louis Bentley在ACM杂志上发表的一篇论文:Multidimensional Binary Search Trees Used for Associative Searching 中正式提出和阐述的了如下图形式的把空间划分为多个部分的k-d树。 2.1、什么是KD树Kd-树是K-dimension tree的缩写,是对数据点在k维空间(如二维(x,y),三维(x,y,z),k维(x1,y,z..))中划分的一种数据结构,主要应用于多维空间关键数据的搜索(如:范围搜索和最近邻搜索)。本质上说,Kd-树就是一种平衡二叉树。 首先必须搞清楚的是,k-d树是一种空间划分树,说白了,就是把整个空间划分为特定的几个部分,然后在特定空间的部分内进行相关搜索操作。想像一个三维(多维有点为难你的想象力了)空间,kd树按照一定的划分规则把这个三维空间划分了多个空间,如下图所示: 2.2、KD树的构建kd树构建的伪代码如下图所示(伪代码来自《图像局部不变特性特征与描述》王永明 王贵锦 编著): 再举一个简单直观的实例来介绍k-d树构建算法。假设有6个二维数据点{(2,3),(5,4),(9,6),(4,7),(8,1),(7,2)},数据点位于二维空间内,如下图所示。为了能有效的找到最近邻,k-d树采用分而治之的思想,即将整个空间划分为几个小部分,首先,粗黑线将空间一分为二,然后在两个子空间中,细黑直线又将整个空间划分为四部分,最后虚黑直线将这四部分进一步划分。 6个二维数据点{(2,3),(5,4),(9,6),(4,7),(8,1),(7,2)}构建kd树的具体步骤为:
如上算法所述,kd树的构建是一个递归过程,我们对左子空间和右子空间内的数据重复根节点的过程就可以得到一级子节点(5,4)和(9,6),同时将空间和数据集进一步细分,如此往复直到空间中只包含一个数据点。 与此同时,经过对上面所示的空间划分之后,我们可以看出,点(7,2)可以为根结点,从根结点出发的两条红粗斜线指向的(5,4)和(9,6)则为根结点的左右子结点,而(2,3),(4,7)则为(5,4)的左右孩子(通过两条细红斜线相连),最后,(8,1)为(9,6)的左孩子(通过细红斜线相连)。如此,便形成了下面这样一棵k-d树: k-d树的数据结构 针对上表给出的kd树的数据结构,转化成具体代码如下所示(注,本文以下代码分析基于Rob Hess维护的sift库): /** a node in a k-d tree */struct kd_node{ int ki; /**< partition="" key="" index="" *///关键点直方图方差最大向量系列位置="" double="" kv;="">< partition="" key="" value="" *///直方图方差最大向量系列中最中间模值="" int="" leaf;="">< 1="" if="" node="" is="" a="" leaf,="" 0="" otherwise="" */="" struct="" feature*="" features;="">< features="" at="" this="" node="" */="" int="" n;="">< number="" of="" features="" */="" struct="" kd_node*="" kd_left;="">< left="" child="" */="" struct="" kd_node*="" kd_right;="">< right="" child=""> 也就是说,如之前所述,kd树中,kd代表k-dimension,每个节点即为一个k维的点。每个非叶节点可以想象为一个分割超平面,用垂直于坐标轴的超平面将空间分为两个部分,这样递归的从根节点不停的划分,直到没有实例为止。经典的构造k-d tree的规则如下:
对于n个实例的k维数据来说,建立kd-tree的时间复杂度为O(k*n*logn)。 以下是构建k-d树的代码: struct kd_node* kdtree_build( struct feature* features, int n ){ struct kd_node* kd_root; if( ! features || n <= 0="" )="" {="" fprintf(="" stderr,="" 'warning:="" kdtree_build():="" no="" features,="" %s,="" line="" %d\n',="" __file__,="" __line__="" );="" return="" null;="" }="" 初始化="" kd_root="kd_node_init(" features,="" n="" );="" n--number="" of="" features,initinalize="" root="" of="" tree.="" expand_kd_node_subtree(="" kd_root="" );="" kd="" tree="" expand="" return="">=> 上面的涉及初始化操作的两个函数kd_node_init,及expand_kd_node_subtree代码分别如下所示: static struct kd_node* kd_node_init( struct feature* features, int n ){ //n--number of features struct kd_node* kd_node; kd_node = (struct kd_node*)(malloc( sizeof( struct kd_node ) )); memset( kd_node, 0, sizeof( struct kd_node ) ); //0填充 kd_node->ki = -1; //??????? kd_node->features = features; kd_node->n = n; return kd_node;} static void expand_kd_node_subtree( struct kd_node* kd_node ){ /* base case: leaf node */ if( kd_node->n == 1 || kd_node->n == 0 ) { //叶节点 //伪叶节点 kd_node->leaf = 1; return; } assign_part_key( kd_node ); //get ki,kv partition_features( kd_node ); //creat left and right children,特征点ki位置左树比右树模值小,kv作为分界模值 //kd_node中关键点已经排序 if( kd_node->kd_left ) expand_kd_node_subtree( kd_node->kd_left ); if( kd_node->kd_right ) expand_kd_node_subtree( kd_node->kd_right );} 构建完kd树之后,如今进行最近邻搜索呢?从下面的动态gif图中,你是否能看出些许端倪呢? k-d树算法可以分为两大部分,除了上部分有关k-d树本身这种数据结构建立的算法,另一部分是在建立的k-d树上各种诸如插入,删除,查找(最邻近查找)等操作涉及的算法。下面,咱们依次来看kd树的插入、删除、查找操作。 2.3、KD树的插入元素插入到一个K-D树的方法和二叉检索树类似。本质上,在偶数层比较x坐标值,而在奇数层比较y坐标值。当我们到达了树的底部,(也就是当一个空指针出现),我们也就找到了结点将要插入的位置。生成的K-D树的形状依赖于结点插入时的顺序。给定N个点,其中一个结点插入和检索的平均代价是O(log2N)。 下面4副图(来源:中国地质大学电子课件)说明了插入顺序为(a) Chicago, (b) Mobile, (c) Toronto, and (d) Buffalo,建立空间K-D树的示例: 应该清楚,这里描述的插入过程中,每个结点将其所在的平面分割成两部分。因比,Chicago 将平面上所有结点分成两部分,一部分所有的结点x坐标值小于35,另一部分结点的x坐标值大于或等于35。同样Mobile将所有x坐标值大于35的结点以分成两部分,一部分结点的Y坐标值是小于10,另一部分结点的Y坐标值大于或等于10。后面的Toronto、Buffalo也按照一分为二的规则继续划分。 2.4、KD树的删除 KD树的删除可以用递归程序来实现。我们假设希望从K-D树中删除结点(a,b)。如果(a,b)的两个子树都为空,则用空树来代替(a,b)。否则,在(a,b)的子树中寻找一个合适的结点来代替它,譬如(c,d),则递归地从K-D树中删除(c,d)。一旦(c,d)已经被删除,则用(c,d)代替(a,b)。假设(a,b)是一个X识别器,那么,它得替代节点要么是(a,b)左子树中的X坐标最大值的结点,要么是(a,b)右子树中x坐标最小值的结点。 也就是说,跟普通二叉树(包括如下图所示的红黑树)结点的删除是同样的思想:用被删除节点A的左子树的最右节点或者A的右子树的最左节点作为替代A的节点(比如,下图红黑树中,若要删除根结点26,第一步便是用23或28取代根结点26)。 当(a,b)的右子树为空时,找到(a,b)左子树中具有x坐标最大的结点,譬如(c,d),将(a,b)的左子树放到(c,d)的右子树中,且在树中从它的上一层递归地应用删除过程(也就是(a,b)的左子树) 。 下面来举一个实际的例子(来源:中国地质大学电子课件,原课件错误已经在下文中订正),如下图所示,原始图像及对应的kd树,现在要删除图中的A结点,请看一系列删除步骤: 要删除上图中结点A,选择结点A的右子树中X坐标值最小的结点,这里是C,C成为根,如下图: 从C的右子树中找出一个结点代替先前C的位置, 这里是D,并将D的左子树转为它的右子树,D代替先前C的位置,如下图: 在D的新右子树中,找X坐标最小的结点,这里为H,H代替D的位置, 在D的右子树中找到一个Y坐标最小的值,这里是I,将I代替原先H的位置,从而A结点从图中顺利删除,如下图所示: 从一个K-D树中删除结点(a,b)的问题变成了在(a,b)的子树中寻找x坐标为最小的结点。不幸的是寻找最小x坐标值的结点比二叉检索树中解决类似的问题要复杂得多。特别是虽然最小x坐标值的结点一定在x识别器的左子树中,但它同样可在y识别器的两个子树中。因此关系到检索,且必须注意检索坐标,以使在每个奇数层仅检索2个子树中的一个。 从K-D树中删除一个结点是代价很高的,很清楚删除子树的根受到子树中结点个数的限制。用TPL(T)表示树T总的路径长度。可看出树中子树大小的总和为TPL(T)+N。 以随机方式插入N个点形成树的TPL是O(N*log2N),这就意味着从一个随机形成的K-D树中删除一个随机选取的结点平均代价的上界是O(log2N) 。 2.5、KD树的最近邻搜索算法现实生活中有许多问题需要在多维数据的快速分析和快速搜索,对于这个问题最常用的方法是所谓的kd树。在k-d树中进行数据的查找也是特征匹配的重要环节,其目的是检索在k-d树中与查询点距离最近的数据点。在一个N维的笛卡儿空间在两个点之间的距离是由下述公式确定: 2.5.1、k-d树查询算法的伪代码 k-d树查询算法的伪代码如下所示: 算法:k-d树最邻近查找输入:Kd, //k-d tree类型 target //查询数据点输出:nearest, //最邻近数据点 dist //最邻近数据点和查询点间的距离1. If Kd为NULL,则设dist为infinite并返回2. //进行二叉查找,生成搜索路径 Kd_point = &Kd; //Kd-point中保存k-d tree根节点地址 nearest = Kd_point -> Node-data; //初始化最近邻点 while(Kd_point) push(Kd_point)到search_path中; //search_path是一个堆栈结构,存储着搜索路径节点指针 If Dist(nearest,target) > Dist(Kd_point -> Node-data,target) nearest = Kd_point -> Node-data; //更新最近邻点 Min_dist = Dist(Kd_point,target); //更新最近邻点与查询点间的距离 ***/ s = Kd_point -> split; //确定待分割的方向 If target[s] <= kd_point="" -=""> Node-data[s] //进行二叉查找 Kd_point = Kd_point -> left; else Kd_point = Kd_point ->right; End while3. //回溯查找 while(search_path != NULL) back_point = 从search_path取出一个节点指针; //从search_path堆栈弹栈 s = back_point -> split; //确定分割方向 If Dist(target[s],back_point -> Node-data[s]) < max_dist="" 判断还需进入的子空间="" if="" target[s]=""><= back_point="" -=""> Node-data[s] Kd_point = back_point -> right; //如果target位于左子空间,就应进入右子空间 else Kd_point = back_point -> left; //如果target位于右子空间,就应进入左子空间 将Kd_point压入search_path堆栈; If Dist(nearest,target) > Dist(Kd_Point -> Node-data,target) nearest = Kd_point -> Node-data; //更新最近邻点 Min_dist = Dist(Kd_point -> Node-data,target); //更新最近邻点与查询点间的距离的 End while =>=> 读者来信点评@yhxyhxyhx,在“将Kd_point压入search_path堆栈;”这行代码后,应该是调到步骤2再往下走二分搜索的逻辑一直到叶结点,我写了一个递归版本的二维kd tree的搜索函数你对比的看看: void innerGetClosest(NODE* pNode, PT point, PT& res, int& nMinDis){ if (NULL == pNode) return; int nCurDis = abs(point.x - pNode->pt.x) + abs(point.y - pNode->pt.y); if (nMinDis < 0="" ||="" ncurdis="">< nmindis)="" {="" nmindis="nCurDis;" res="pNode-">pt; } if (pNode->splitX && point.x <= pnode-="">pt.x || !pNode->splitX && point.y <= pnode-="">pt.y) innerGetClosest(pNode->pLft, point, res, nMinDis); else innerGetClosest(pNode->pRgt, point, res, nMinDis); int rang = pNode->splitX ? abs(point.x - pNode->pt.x) : abs(point.y - pNode->pt.y); if (rang > nMinDis) return; NODE* pGoInto = pNode->pLft; if (pNode->splitX && point.x > pNode->pt.x || !pNode->splitX && point.y > pNode->pt.y) pGoInto = pNode->pRgt; innerGetClosest(pGoInto, point, res, nMinDis);}=>=> 下面,以两个简单的实例(例子来自图像局部不变特性特征与描述一书)来描述最邻近查找的基本思路。 2.5.2、举例:查询点(2.1,3.1)星号表示要查询的点(2.1,3.1)。通过二叉搜索,顺着搜索路径很快就能找到最邻近的近似点,也就是叶子节点(2,3)。而找到的叶子节点并不一定就是最邻近的,最邻近肯定距离查询点更近,应该位于以查询点为圆心且通过叶子节点的圆域内。为了找到真正的最近邻,还需要进行相关的‘回溯'操作。也就是说,算法首先沿搜索路径反向查找是否有距离查询点更近的数据点。 以查询(2.1,3.1)为例:
2.5.3、举例:查询点(2,4.5)一个复杂点了例子如查找点为(2,4.5),具体步骤依次如下:
上述两次实例表明,当查询点的邻域与分割超平面两侧空间交割时,需要查找另一侧子空间,导致检索过程复杂,效率下降。 一般来讲,最临近搜索只需要检测几个叶子结点即可,如下图所示: 但是,如果当实例点的分布比较糟糕时,几乎要遍历所有的结点,如下所示: 研究表明N个节点的K维k-d树搜索过程时间复杂度为:tworst=O(kN1-1/k)。 同时,以上为了介绍方便,讨论的是二维或三维情形。但在实际的应用中,如SIFT特征矢量128维,SURF特征矢量64维,维度都比较大,直接利用k-d树快速检索(维数不超过20)的性能急剧下降,几乎接近贪婪线性扫描。假设数据集的维数为D,一般来说要求数据的规模N满足N?2D,才能达到高效的搜索。所以这就引出了一系列对k-d树算法的改进:BBF算法,和一系列M树、VP树、MVP树等高维空间索引树(下文2.6节kd树近邻搜索算法的改进:BBF算法,与2.7节球树、M树、VP树、MVP树)。 2.6、kd树近邻搜索算法的改进:BBF算法咱们顺着上一节的思路,参考《统计学习方法》一书上的内容,再来总结下kd树的最近邻搜索算法: 输入:以构造的kd树,目标点x;输出:x 的最近邻 算法步骤如下:
如果实例点是随机分布的,那么kd树搜索的平均计算复杂度是O(logN),这里的N是训练实例树。所以说,kd树更适用于训练实例数远大于空间维数时的k近邻搜索,当空间维数接近训练实例数时,它的效率会迅速下降,一降降到“解放前”:线性扫描的速度。 也正因为上述k最近邻搜索算法的第4个步骤中的所述:“回退到根结点时,搜索结束”,每个最近邻点的查询比较完成过程最终都要回退到根结点而结束,而导致了许多不必要回溯访问和比较到的结点,这些多余的损耗在高维度数据查找的时候,搜索效率将变得相当之地下,那有什么办法可以改进这个原始的kd树最近邻搜索算法呢? 从上述标准的kd树查询过程可以看出其搜索过程中的“回溯”是由“查询路径”决定的,并没有考虑查询路径上一些数据点本身的一些性质。一个简单的改进思路就是将“查询路径”上的结点进行排序,如按各自分割超平面(也称bin)与查询点的距离排序,也就是说,回溯检查总是从优先级最高(Best Bin)的树结点开始。 针对此BBF机制,读者Feng&书童点评道:
如此,就引出了本节要讨论的kd树最近邻搜索算法的改进:BBF(Best-Bin-First)查询算法,它是由发明sift算法的David Lowe在1997的一篇文章中针对高维数据提出的一种近似算法,此算法能确保优先检索包含最近邻点可能性较高的空间,此外,BBF机制还设置了一个运行超时限定。采用了BBF查询机制后,kd树便可以有效的扩展到高维数据集上。 伪代码如下图所示(图取自图像局部不变特性特征与描述一书): 还是以上面的查询(2,4.5)为例,搜索的算法流程为:
2.7、球树、M树、VP树、MVP树2.7.1、球树咱们来针对上文内容总结回顾下,针对下面这样一棵kd树: 现要找它的最近邻。 通过上文2.5节,总结来说,我们已经知道: 1、为了找到一个给定目标点的最近邻,需要从树的根结点开始向下沿树找出目标点所在的区域,如下图所示,给定目标点,用星号标示,我们似乎一眼看出,有一个点离目标点最近,因为它落在以目标点为圆心以较小长度为半径的虚线圆内,但为了确定是否可能还村庄一个最近的近邻,我们会先检查叶节点的同胞结点,然叶节点的同胞结点在图中所示的阴影部分,虚线圆并不与之相交,所以确定同胞叶结点不可能包含更近的近邻。 2、于是我们回溯到父节点,并检查父节点的同胞结点,父节点的同胞结点覆盖了图中所有横线X轴上的区域。因为虚线圆与右上方的矩形(KD树把二维平面划分成一个一个矩形)相交... 如上,我们看到,KD树是可用于有效寻找最近邻的一个树结构,但这个树结构其实并不完美,当处理不均匀分布的数据集时便会呈现出一个基本冲突:既邀请树有完美的平衡结构,又要求待查找的区域近似方形,但不管是近似方形,还是矩形,甚至正方形,都不是最好的使用形状,因为他们都有角。 什么意思呢?就是说,在上图中,如果黑色的实例点离目标点星点再远一点,那么势必那个虚线圆会如红线所示那样扩大,以致与左上方矩形的右下角相交,既然相交了,那么势必又必须检查这个左上方矩形,而实际上,最近的点离星点的距离很近,检查左上方矩形区域已是多余。于此我们看见,KD树把二维平面划分成一个一个矩形,但矩形区域的角却是个难以处理的问题。 解决的方案就是使用如下图所示的球树: 先从球中选择一个离球的中心最远的点,然后选择第二个点离第一个点最远,将球中所有的点分配到离这两个聚类中心最近的一个上,然后计算每个聚类的中心,以及聚类能够包含它所有数据点所需的最小半径。这种方法的优点是分裂一个包含n个殊绝点的球的成本只是随n呈线性增加。 使用球树找出给定目标点的最近邻方法是,首先自上而下贯穿整棵树找出包含目标点所在的叶子,并在这个球里找出与目标点最靠近的点,这将确定出目标点距离它的最近邻点的一个上限值,然后跟KD树查找一样,检查同胞结点,如果目标点到同胞结点中心的距离超过同胞结点的半径与当前的上限值之和,那么同胞结点里不可能存在一个更近的点;否则的话,必须进一步检查位于同胞结点以下的子树。 如下图,目标点还是用一个星表示,黑色点是当前已知的的目标点的最近邻,灰色球里的所有内容将被排除,因为灰色球的中心点离的太远,所以它不可能包含一个更近的点,像这样,递归的向树的根结点进行回溯处理,检查所有可能包含一个更近于当前上限值的点的球。 球树是自上而下的建立,和KD树一样,根本问题就是要找到一个好的方法将包含数据点集的球分裂成两个,在实践中,不必等到叶子结点只有两个胡数据点时才停止,可以采用和KD树一样的方法,一旦结点上的数据点打到预先设置的最小数量时,便可提前停止建树过程。 也就是上面所述,先从球中选择一个离球的中心最远的点,然后选择第二个点离第一个点最远,将球中所有的点分配到离这两个聚类中心最近的一个上,然后计算每个聚类的中心,以及聚类能够包含它所有数据点所需的最小半径。这种方法的优点是分裂一个包含n个殊绝点的球的成本只是随n呈线性增加(注:本小节内容主要来自参考条目19:数据挖掘实用机器学习技术,[新西兰]Ian H.Witten 著,第4章4.7节)。 2.7.2、VP树与MVP树简介高维特征向量的距离索引问题是基于内容的图像检索的一项关键技术,目前经常采用的解决办法是首先对高维特征空间做降维处理,然后采用包括四叉树、kd树、R树族等在内的主流多维索引结构,这种方法的出发点是:目前的主流多维索引结构在处理维数较低的情况时具有比较好的效率,但对于维数很高的情况则显得力不从心(即所谓的维数危机) 。 实验结果表明当特征空间的维数超过20 的时候,效率明显降低,而可视化特征往往采用高维向量描述,一般情况下可以达到10^2的量级,甚至更高。在表示图像可视化特征的高维向量中各维信息的重要程度是不同的,通过降维技术去除属于次要信息的特征向量以及相关性较强的特征向量,从而降低特征空间的维数,这种方法已经得到了一些实际应用。 然而这种方法存在不足之处采用降维技术可能会导致有效信息的损失,尤其不适合于处理特征空间中的特征向量相关性很小的情况。另外主流的多维索引结构大都针对欧氏空间,设计需要利用到欧氏空间的几何性质,而图像的相似性计算很可能不限于基于欧氏距离。这种情况下人们越来越关注基于距离的度量空间高维索引结构可以直接应用于高维向量相似性查询问题。 度量空间中对象之间的距离度量只能利用三角不等式性质,而不能利用其他几何性质。向量空间可以看作由实数坐标串组成的特殊度量空间,目前针对度量空间的高维索引问题提出的索引结构有很多种大致可以作如下分类,如下图所示: 其中,VP树和MVP树中特征向量的举例表示为: 读者点评:
更多内容请参见论文1:DIST ANCE-BASED INDEXING FOR HIGH-DIMENSIONAL METRIC SP ACES,作者:Tolga Bozkaya & Meral Ozsoyoglu,及论文2:基于度量空间高维索引结构VP-tree及MVP-tree的图像检索,王志强,甘国辉,程起敏。 当然,如果你觉得上述论文还不够满足你胃口的话,这里有一大堆nearest neighbor algorithms相关的论文可供你看:http://scholar./scholar?q=nearest+neighbor+algorithms&btnG=&hl=zh-CN&as_sdt=0&as_vis=1(其中,这篇可以看下:Spill-Trees,An investigation of practical approximate nearest neighbor algorithms)。 第三部分、KD树的应用:SIFT+KD_BBF搜索算法3.1、SIFT特征匹配算法之前本blog内阐述过图像特征匹配SIFT算法,写过五篇文章,这五篇文章分别为:
不熟悉SIFT算法相关概念的可以看上述几篇文章,这里不再做赘述。与此同时,本文此部分也作为十五个经典算法研究系列里SIFT算法的九之四续。 OK,我们知道,在sift算法中,给定两幅图片图片,若要做特征匹配,一般会先提取出图片中的下列相关属性作为特征点: /**Structure to represent an affine invariant image feature. The fieldsx, y, a, b, c represent the affine region around the feature:a(x-u)(x-u) + 2b(x-u)(y-v) + c(y-v)(y-v) = 1*/struct feature{ double x; /**< x="" coord="" */="" double="" y;="">< y="" coord="" */="" double="" a;="">< oxford-type="" affine="" region="" parameter="" */="" double="" b;="">< oxford-type="" affine="" region="" parameter="" */="" double="" c;="">< oxford-type="" affine="" region="" parameter="" */="" double="" scl;="">< scale="" of="" a="" lowe-style="" feature="" */="" double="" ori;="">< orientation="" of="" a="" lowe-style="" feature="" */="" int="" d;="">< descriptor="" length="" */="" double="" descr[feature_max_d];="">< descriptor="" */="" int="" type;="">< feature="" type,="" oxfd="" or="" lowe="" */="" int="" category;="">< all-purpose="" feature="" category="" */="" struct="" feature*="" fwd_match;="">< matching="" feature="" from="" forward="" image="" */="" struct="" feature*="" bck_match;="">< matching="" feature="" from="" backmward="" image="" */="" struct="" feature*="" mdl_match;="">< matching="" feature="" from="" model="" */="" cvpoint2d64f="" img_pt;="">< location="" in="" image="" */="" cvpoint2d64f="" mdl_pt;="">< location="" in="" model="" */="" void*="" feature_data;="">< user-definable="" data="" */="" char="" dense;=""> 而后在sift.h文件中定义两个关键函数,这里,我们把它们称之为函数一,和函数二,如下所示:函数一的声明: extern int sift_features( IplImage* img, struct feature** feat ); 函数一的实现: int sift_features( IplImage* img, struct feature** feat ){ return _sift_features( img, feat, SIFT_INTVLS, SIFT_SIGMA, SIFT_CONTR_THR, SIFT_CURV_THR, SIFT_IMG_DBL, SIFT_DESCR_WIDTH, SIFT_DESCR_HIST_BINS );} 从上述函数一的实现中,我们可以看到,它内部实际上调用的是这个函数:_sift_features(..),也就是下面马上要分析的函数二。函数二的声明: extern int _sift_features( IplImage* img, struct feature** feat, int intvls, double sigma, double contr_thr, int curv_thr, int img_dbl, int descr_width, int descr_hist_bins ); 函数二的实现: int _sift_features( IplImage* img, struct feature** feat, int intvls, double sigma, double contr_thr, int curv_thr, int img_dbl, int descr_width, int descr_hist_bins ){ IplImage* init_img; IplImage*** gauss_pyr, *** dog_pyr; CvMemStorage* storage; CvSeq* features; int octvs, i, n = 0,n0 = 0,n1 = 0,n2 = 0,n3 = 0,n4 = 0; int start; /* check arguments */ if( ! img ) fatal_error( 'NULL pointer error, %s, line %d', __FILE__, __LINE__ ); if( ! feat ) fatal_error( 'NULL pointer error, %s, line %d', __FILE__, __LINE__ ); /* build scale space pyramid; smallest dimension of top level is ~4 pixels */ start=GetTickCount(); init_img = create_init_img( img, img_dbl, sigma ); octvs = log( (float)(MIN( init_img->width, init_img->height )) ) / log((float)(2.0)) -5; gauss_pyr = build_gauss_pyr( init_img, octvs, intvls, sigma ); dog_pyr = build_dog_pyr( gauss_pyr, octvs, intvls ); fprintf( stderr, ' creat the pyramid use %d\n',GetTickCount()-start); storage = cvCreateMemStorage( 0 ); //创建存储内存,0为默认64k start=GetTickCount(); features = scale_space_extrema( dog_pyr, octvs, intvls, contr_thr, curv_thr, storage ); //在DOG空间寻找极值点,确定关键点位置 fprintf( stderr, ' find the extrum points in DOG use %d\n',GetTickCount()-start); calc_feature_scales( features, sigma, intvls ); //计算关键点的尺度 if( img_dbl ) adjust_for_img_dbl( features ); //如果原始空间图扩大,特征点坐标就缩小 start=GetTickCount(); calc_feature_oris( features, gauss_pyr ); //在gaussian空间计算关键点的主方向和幅值 fprintf( stderr, ' get the main oritation use %d\n',GetTickCount()-start); start=GetTickCount(); compute_descriptors( features, gauss_pyr, descr_width, descr_hist_bins ); //建立关键点描述器 fprintf( stderr, ' compute the descriptors use %d\n',GetTickCount()-start); /* sort features by decreasing scale and move from CvSeq to array */ //start=GetTickCount(); //cvSeqSort( features, (CvCmpFunc)feature_cmp, NULL ); //????? n = features->total; *feat = (feature*)(calloc( n, sizeof(struct feature) )); *feat = (feature*)(cvCvtSeqToArray( features, *feat, CV_WHOLE_SEQ )); //整条链表放在feat指向的内存 for( i = 0; i < n;="" i++="" )="" {="" free(="" (*feat)[i].feature_data="" );="" (*feat)[i].feature_data="NULL;" 释放ddata(r,c,octv,intvl,xi,scl_octv)="" if((*feat)[i].dense="=" 4)="" ++n4;="" else="" if((*feat)[i].dense="=" 3)="" ++n3;="" else="" if((*feat)[i].dense="=" 2)="" ++n2;="" else="" if((*feat)[i].dense="=" 1)="" ++n1;="" else="" ++n0;="" }="" fprintf(="" stderr,="" '="" move="" features="" from="" sequce="" to="" array="" use="" %d\n',gettickcount()-start);="" start="GetTickCount();" fprintf(="" stderr,="" 'in="" the="" total="" feature="" points="" the="" extent4="" points="" is="" %d\n',n4);="" fprintf(="" stderr,="" 'in="" the="" total="" feature="" points="" the="" extent3="" points="" is="" %d\n',n3);="" fprintf(="" stderr,="" 'in="" the="" total="" feature="" points="" the="" extent2="" points="" is="" %d\n',n2);="" fprintf(="" stderr,="" 'in="" the="" total="" feature="" points="" the="" extent1="" points="" is="" %d\n',n1);="" fprintf(="" stderr,="" 'in="" the="" total="" feature="" points="" the="" extent0="" points="" is="" %d\n',n0);="" cvreleasememstorage(="" &storage="" );="" cvreleaseimage(="" &init_img="" );="" release_pyr(="" &gauss_pyr,="" octvs,="" intvls="" +="" 3="" );="" release_pyr(="" &dog_pyr,="" octvs,="" intvls="" +="" 2="" );="" fprintf(="" stderr,="" '="" free="" the="" pyramid="" use="" %d\n',gettickcount()-start);="" return=""> 说明:上面的函数二,包含了SIFT算法中几乎所有函数,是SIFT算法的核心。本文不打算一一分析上面所有函数,只会抽取其中涉及到BBF查询机制相关的函数。3.2、K个最小近邻的查找:大顶堆优先级队列上文中一直在讲最近邻问题,也就是说只找最近的那唯一一个邻居,但如果现实中需要我们找到k个最近的邻居。该如何做呢?对的,之前blog内曾相近阐述过寻找最小的k个数的问题,显然,寻找k个最近邻与寻找最小的k个数的问题如出一辙。 回忆下寻找k个最小的数中关于构造大顶堆的解决方案: “寻找最小的k个树,更好的办法是维护k个元素的最大堆,即用容量为k的最大堆存储最先遍历到的k个数,并假设它们即是最小的k个数,建堆费时O(k)后,有k1<><><><>O(n*logk)。此方法得益于在堆中,查找等各项操作时间复杂度均为logk。” 根据上述方法,咱们来写大顶堆最小优先级队列相关代码实现: void* minpq_extract_min( struct min_pq* min_pq ){ void* data; if( min_pq->n < 1="" )="" {="" fprintf(="" stderr,="" 'warning:="" pq="" empty,="" %s="" line="" %d\n',="" __file__,="" __line__="" );="" return="" null;="" }="" data="min_pq-">pq_array[0].data; //root of tree min_pq->n--; //0 min_pq->pq_array[0] = min_pq->pq_array[min_pq->n]; restore_minpq_order( min_pq->pq_array, 0, min_pq->n ); //0 return data;} 上述函数中,restore_minpq_order的实现如下:static void restore_minpq_order( struct pq_node* pq_array, int i, int n ){ //0 //0 struct pq_node tmp; int l, r, min = i; l = left( i ); //2*i+1,???????????? r = right( i );//2*i+2,???????????? if( l < n="" )="" if(="" pq_array[l].key="">< pq_array[i].key="" )="" min="l;" if(="" r="">< n="" )="" if(="" pq_array[r].key="">< pq_array[min].key="" )="" min="r;" if(="" min="" !="i" )="" {="" tmp="pq_array[min];" pq_array[min]="pq_array[i];" pq_array[i]="tmp;" restore_minpq_order(="" pq_array,="" min,="" n="" );=""> 3.3、KD树近邻搜索改进之BBF算法原理在上文第二部分已经阐述明白,结合大顶堆找最近的K个近邻思想,相关主体代码如下: //KD树近邻搜索改进之BBF算法int kdtree_bbf_knn( struct kd_node* kd_root, struct feature* feat, int k, struct feature*** nbrs, int max_nn_chks ) //2{ //200 struct kd_node* expl; struct min_pq* min_pq; struct feature* tree_feat, ** _nbrs; struct bbf_data* bbf_data; int i, t = 0, n = 0; if( ! nbrs || ! feat || ! kd_root ) { fprintf( stderr, 'Warning: NULL pointer error, %s, line %d\n', __FILE__, __LINE__ ); return -1; } _nbrs = (feature**)(calloc( k, sizeof( struct feature* ) )); //2 min_pq = minpq_init(); minpq_insert( min_pq, kd_root, 0 ); //把根节点加入搜索序列中 //队列有东西就继续搜,同时控制在t<200步内 while(="" min_pq-="">n > 0 && t < max_nn_chks="" )="" {="" 刚进来时,从kd树根节点搜索,exp1是根节点="" 后进来时,exp1是min_pq差值最小的未搜索节点入口="" 同时按min_pq中父,子顺序依次检验,保证父节点的差值比子节点小.这样减少返回搜索时间="" expl="(struct" kd_node*)minpq_extract_min(="" min_pq="" );="" if(="" !="" expl="" )="" {="" fprintf(="" stderr,="" 'warning:="" pq="" unexpectedly="" empty,="" %s="" line="" %d\n',="" __file__,="" __line__="" );="" goto="" fail;="" }="" 从根节点(或差值最小节点)搜索,根据目标点与节点模值的差值(小)="" 确定在kd树的搜索路径,同时存储各个节点另一入口地址\同级搜索路径差值.="" 存储时比较父节点的差值,如果子节点差值比父节点差值小,交换两者存储位置,="" 使未搜索节点差值小的存储在min_pq的前面,减小返回搜索的时间.="" expl="explore_to_leaf(" expl,="" feat,="" min_pq="" );="" if(="" !="" expl="" )="" {="" fprintf(="" stderr,="" 'warning:="" pq="" unexpectedly="" empty,="" %s="" line="" %d\n',="" __file__,="" __line__="" );="" goto="" fail;="" }="" for(="" i="0;" i="">< expl-="">n; i++ ) { //使用exp1->n原因:如果是叶节点,exp1->n=1,如果是伪叶节点,exp1->n=0. tree_feat = &expl->features[i]; bbf_data = (struct bbf_data*)(malloc( sizeof( struct bbf_data ) )); if( ! bbf_data ) { fprintf( stderr, 'Warning: unable to allocate memory,' ' %s line %d\n', __FILE__, __LINE__ ); goto fail; } bbf_data->old_data = tree_feat->feature_data; bbf_data->d = descr_dist_sq(feat, tree_feat); //计算两个关键点描述器差平方和 tree_feat->feature_data = bbf_data; //取前k个 n += insert_into_nbr_array( tree_feat, _nbrs, n, k );// } //2 t++; } minpq_release( &min_pq ); for( i = 0; i < n;="" i++="" )="" bbf_data为何搞个old_data?="" {="" bbf_data="(struct" bbf_data*)(_nbrs[i]-="">feature_data); _nbrs[i]->feature_data = bbf_data->old_data; free( bbf_data ); } *nbrs = _nbrs; return n;fail: minpq_release( &min_pq ); for( i = 0; i < n;="" i++="" )="" {="" bbf_data="(struct" bbf_data*)(_nbrs[i]-="">feature_data); _nbrs[i]->feature_data = bbf_data->old_data; free( bbf_data ); } free( _nbrs ); *nbrs = NULL; return -1;}200步内> 依据上述函数kdtree_bbf_knn从上而下看下来,注意几点: 1、上述函数kdtree_bbf_knn中,explore_to_leaf的代码如下: static struct kd_node* explore_to_leaf( struct kd_node* kd_node, struct feature* feat, struct min_pq* min_pq ) //kd tree //the before { struct kd_node* unexpl, * expl = kd_node; double kv; int ki; while( expl && ! expl->leaf ) { //0 ki = expl->ki; kv = expl->kv; if( ki >= feat->d ) { fprintf( stderr, 'Warning: comparing imcompatible descriptors, %s' \ ' line %d\n', __FILE__, __LINE__ ); return NULL; } if( feat->descr[ki] <= kv="" )="" 由目标点描述器确定搜索向kd="" tree的左边或右边="" {="" 如果目标点模值比节点小,搜索向tree的左边进行="" unexpl="expl-">kd_right; expl = expl->kd_left; } else { unexpl = expl->kd_left; //else expl = expl->kd_right; } //把未搜索数分支入口,差值存储在min_pq, if( minpq_insert( min_pq, unexpl, ABS( kv - feat->descr[ki] ) ) ) { fprintf( stderr, 'Warning: unable to insert into PQ, %s, line %d\n', __FILE__, __LINE__ ); return NULL; } } return expl;}=> 2、上述查找函数kdtree_bbf_knn中的参数k可调,代表的是要查找近邻的个数,即number of neighbors to find,在sift特征匹配中,k一般取2 k = kdtree_bbf_knn( kd_root_0, feat, 2, &nbrs, KDTREE_BBF_MAX_NN_CHKS_0 );//点匹配函数(核心) if( k == 2 ) //只有进行2次以上匹配过程,才算是正常匹配过程 3、上述函数kdtree_bbf_knn中“bbf_data->d = descr_dist_sq(feat, tree_feat); //计算两个关键点描述器差平方和”,使用的计算方法是本文第一部分1.2节中所述的欧氏距离。 3.3、SIFT+BBF算法匹配效果之前试了下sift + KD + BBF算法,用两幅不同的图片做了下匹配(当然,运行结果显示是不匹配的),效果还不错:http://weibo.com/1580904460/yDmzAEwcV#1348475194313。 “编译的过程中,直接用的VS2010 + opencv(并没下gsl)。2012.09.24”。.... 本文完整源码有pudn账号的朋友可以前去这里下载:http://www./downloads340/sourcecode/graph/texture_mapping/detail1486667.html(没有pudn账号的同学请加群:169056165,至群共享下载,验证信息:sift)。感谢诸位。 参考文献及推荐阅读
后记 从当天下午4点多一直写,一直写,直接写到了隔日凌晨零点左右,或参考或直接间接引用了很多人的作品及一堆wikipedia页面(当然,已经注明在上面参考文献及推荐阅读中前半部分条目,后半部分条目则为行文之后看到的一些好的资料文章,可以课外读读),但本文还是有诸多地方有待修补完善,也欢迎广大读者不吝赐教 & 指正,感谢大家。完。July、二零一二年十一月二十一日零点五分。
|
|
来自: shiyiyuting > 《公共》