注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

BeyondEgo

Welcome to Winsolider's yard! 超越自我,谁与争锋?

 
 
 

日志

 
 
关于我

本博为记事、畅聊、交友博客,邀你共同探讨人生、探讨成长,广交天下有志之士!愿与有相同兴趣爱好的你,共同学习、一起成长、收获喜悦!

网易考拉推荐

【转摘】SIFT算法分析(草稿)  

2012-02-27 22:50:59|  分类: 程序员 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

特征提取在CV(computer vision)领域非常重要。SIFT是非常出名的特征提取算法,它来自论文IJCV'04的“Distinctive image features from scale-invariant keypoints”,在scholar.google.com上查到的引用次数一万多次,很高了!我准备在这个帖子里,根据这篇论文和SIFT算法的一个开源实现,详细描述SIFT算法。
本文中的SIFT实现在http://blogs.oregonstate.edu/hess/,http://blogs.oregonstate.edu/hess/code/sift/,Linux版本的。

0. 摘要
这里说SIFT算法如何之好,在图像有尺度变化和旋转的时候,提取的特征是不变的,在一定程度的仿射扭曲,3D视角,噪声,照明变化的情况下,匹配稳定。
1. 介绍
SIFT算法用级联式的过滤方法,降低提取特征的时间成本。也就是说,它用4个步骤提取特征,对图像上的每个像素,只有通过了着4个步骤,才算是特征点。这样做的好处就是,把大部分的像素在开始就排除掉,不在它们那里浪费时间。这4个步骤如下:1)尺度空间极值;2)关键点定位;3)方向关联;4)关键点描述子。
2. 相关研究
特征提取文献综述。
3. 尺度空间极值的检测
开始进入正题了。尺度空间极值检测是4个步骤的第1个。为什么要检测尺度空间极值呢?如果要提取一个目标的不变特征,假设是我们要识别iphone,就要先给iphone拍张照片,然后用这个照片作为基准,提取特征。提取好特征值后,将来遇到iphone的大图,小图,旋转了一定角度的图,都能识别出来,如果只是大小和角度都不变的情况下识别,没什么实际价值。既然在大图,小图,旋转图的情况下都能识别,那么,必然在图像里必然存在一些“特性”,这些特征在大图,小图,旋转图里都“不变”的,或者近似不变。这个必然的,否则没有特征提取的必要了。尺度空间极值检测,就是找到图像上的若干候选点,这些点在图像的尺度发生变化的时候,也就是图像变大或者变小,它们和周围其他的点的关系保持不变。
引用的论文说,在一定假定下,尺度空间核只可能是高斯函数。公式(1)就给出了用微分高斯函数卷积图像的计算尺度空间极值的方法。注意:这块的算法本身是很简洁,只是写的不太容易看得懂,反正我看论文的话是没看明白如何实现它,第3节的位于3.1小节上的最后那段话写的相当费解。我从源代码出发,把第3节从开始到第3.1节之间的部分,把公式(1)的DOG计算方法写成如下步骤:
[1] 假设待处理的图像是demo.jpg,是单通道的灰度图,尺寸是512*256。
[2] 计算octave
什么是octave呢?octave原来意思是音乐上的一个八度。代码里octave数量计算方法是:
octvs = log( MIN( init_img->width, init_img->height ) ) / log(2) - 2;
显然,如果图像是512*256,按照上式计算octvs=6。
[3] 计算图像高斯金字塔
[3.1] 有两个参数intvls = 3, sigma = 1.6,在代码里是常量。
[3.2] 计算高斯金字塔的函数,返回一个图像高斯金字塔的指针。图像高斯金字塔指针,指向了6个octave指针,每个octave指针指向intvls+3=6个图像指针。这也就是说,图像指针是IplImage*,那么,每个octave指针就是指针的指针,就是IplImage **,而高斯金字塔的指针是指针的指针的指针,就是IplImage ***。
[3.3] 计算sigma:在每个octave里,有intvls+3=6个图像,每个图像都是前一个图像通过高斯卷积获得的,每次卷积的sigma不同,需要分别计算。sig[0] = sigma = 1.6。k = 2^(1/3)。其他的sig如下:
sig[1] = ( (k^0*1.6*k)^2 - ( k^0*1.6 )^2 )^(0.5)
sig[2] = ( (k^1*1.6*k)^2 - ( k^1*1.6 )^2 )^(0.5)
sig[3] = ( (k^2*1.6*k)^2 - ( k^2*1.6 )^2 )^(0.5)
sig[4] = ( (k^3*1.6*k)^2 - ( k^3*1.6 )^2 )^(0.5)
sig[5] = ( (k^4*1.6*k)^2 - ( k^4*1.6 )^2 )^(0.5)
对每个octave,sig都是相同的。这个sig会用在每次octave的计算中。为什么有这么奇怪的计算方式呢?因为在源代码里,作者是根据同octave的第i个图像进行高斯卷积,计算第i+1个图像,才会有这么怪异的计算sigma,如果是用第1个计算其他的图像,就不会这么麻烦了。
[3.4] 计算过程
[3.4.1]第0层的octave第0个图像是原始图像demo.jpg。
[3.4.2]第0层的octave第1个图像是第0层octave的第0个图像的高斯卷积,卷积参数是sig[1]
[3.4.3]第0层的octave第2个图像是第0层octave的第1个图像的高斯卷积,卷积参数是sig[2]
[3.4.4]第0层的octave第3个图像是第0层octave的第2个图像的高斯卷积,卷积参数是sig[3]
[3.4.5]第0层的octave第4个图像是第0层octave的第3个图像的高斯卷积,卷积参数是sig[4]
[3.4.6]第0层的octave第5个图像是第0层octave的第4个图像的高斯卷积,卷积参数是sig[5]
[3.4.7]第1层的octave第0个图像是第0层octave的第intvls=3个图像的下抽样,宽度和高度均原来一半,即256*128。
[3.4.8]第1层的octave第1个图像是第1层octave的第0个图像的高斯卷积,卷积参数是sig[1]
... ...
以下以此类推,最终计算出整个高斯金字塔。
[4] 计算图像DOG金字塔
[4.1]DOG金字塔指针:DOG金字塔是一个指向指针的指针的指针,即IplImage ***。DOG金字塔指针,指向6个octave指针,每个octave指针,是指向指针的指针,即IplImage **。每个octave指针指向5个图像指针。
[4.2] 计算过程:
[4.2.1]第0个octave第0个图像,是高斯金字塔第0个octave第1个图像减其的第0个octave第0个图像
[4.2.2]第0个octave第1个图像,是高斯金字塔第0个octave第2个图像减其的第0个octave第1个图像
[4.2.3]第0个octave第2个图像,是高斯金字塔第0个octave第3个图像减其的第0个octave第2个图像
[4.2.4]第0个octave第3个图像,是高斯金字塔第0个octave第4个图像减其的第0个octave第3个图像
[4.2.5]第0个octave第4个图像,是高斯金字塔第0个octave第5个图像减其的第0个octave第4个图像
[4.2.6]第1个octave第0个图像,是高斯金字塔第1个octave第1个图像减其的第1个octave第0个图像
[4.2.7]第1个octave第1个图像,是高斯金字塔第1个octave第2个图像减其的第1个octave第1个图像
[4.2.8]第1个octave第2个图像,是高斯金字塔第1个octave第3个图像减其的第1个octave第2个图像
[4.2.9]第1个octave第3个图像,是高斯金字塔第1个octave第4个图像减其的第1个octave第3个图像
... ...
以下以此类推,最终计算出整个DOG金字塔。
3.1 局部极值检测
从DOG金字塔里检测尺度空间极值。DOG金字塔有6个octave。对于每个octave,分别找极值。具体步骤如下:
[5] 计算尺度空间极值点
[5.1] 第0个octave的尺度空间极值
[5.1.1]第0层空间极值:对于第1个图像上的所有像素点的值,如果它比它的8个邻域和第0个图像的9个邻域像素和第2个图像的9个邻域像素都大,或者都小,那么它是一个极值点。
[5.1.2]第1层空间极值:对于第2个图像上的所有像素点的值,如果它比它的8个邻域和第1个图像的9个邻域像素和第3个图像的9个邻域像素都大,或者都小,那么它是一个极值点。
[5.1.3]第2层空间极值:对于第3个图像上的所有像素点的值,如果它比它的8个邻域和第2个图像的9个邻域像素和第4个图像的9个邻域像素都大,或者都小,那么它是一个极值点。
[5.1.4]第3层空间极值:对于第4个图像上的所有像素点的值,如果它比它的8个邻域和第3个图像的9个邻域像素和第5个图像的9个邻域像素都大,或者都小,那么它是一个极值点。
[5.2] 第1个octave的尺度空间极值
... ...
以下以此类推,最终计算出所有的空间极值点。这些极值点是候选点,需要进一步过滤。
3.2 尺度抽样频率
这一段分析抽样的各参数对性能的影响问题。
3.3 空间域的抽样频率
这一段分析抽样的各参数对性能的影响问题。
4 精确的关键点定位
这部分内容是SIFT的第2步。从第4节到第4.1节之间部分,是精确定位关键的坐标,候选集里的所有点,都要进一步处理。每个关键点根据泰勒级数进行二阶展开,它的像素之值是横坐标,纵坐标和sigma的函数,要求它的极值点,就要将二阶泰勒展开,求一阶导数,令导数为0,求出det_r, det_c, det_sigma,那么,就能精确定位极值点的精确的intvls,纵坐标和横坐标。求出精确坐标之后,再计算它的若干属性,如对比度,要超过一定的阈值才行。这就是公式(3)。对每个关键点,计算dD和hessen_3D,然后矩阵相乘,得出det_x,得出精确值。这一块对应的interp_extremum函数,呆定的,熟悉线性代数和微积分的这块好理解,如果不好理解....翻书去看吧。注意,着一块的推导,类似ax^2+bx+c求一阶导数,将一阶偏导和二阶偏导视为常量即可,但如果将它们视为函数,结果是不同的,是(3)的二分之一。
[6]精确定位每个关键点的坐标
[6.1] 对每个关键点,用一个循环,计算关键点精确坐标,最多计算5次。
[6.1.1] 插值计算其关键点极值在intvsl, row, col上的增量xi, xr, xc。这里的计算是根据公式(3)进行的。
[6.1.2] 如果xi, xr, xc的绝对值都小于0.5,表明够精确了,跳出循环。
[6.1.3] 对intvl, r, c 进行更新。如果intvl, r, c超过正常值,就返回NULL。
[6.2] 如果循环超过给定值,返回NULL。
[6.3] 计算关键点的对比度,也就是第11页论文的公式。如果对比度小于某个常量,返回NULL。
[6.4] 生成关键点数据结构feat,根据octave的值,对关键点的坐标进行放大,计算其在原始图像上的像素位置。同时,也把关键点所在的octave, r, c, intvl, subintvl, 并返回feat。
4.1 删除边缘点。
如果关键点是边缘线条上的,那么要去除掉。对所有的关键点都这么处理,一个公式,呆定的。
这个原理是这样的,如果一个点,它在线条边缘上,那么,它的2维hessian矩阵的特征值,必然是一个特征值比另外一个特征值大很多。如果它在两个线条的交点,也就是角点,那么两个特征值在同一个数量级上。对hessian矩阵直接求特征值计算量比较大,因此,论文第12页用一个替代性的方法求解。
[7] 删除边缘上的关键点
对于所有的关键点,即[6]里头得到的每个feat,根据第12页的最后一个公式,判断是否惟边缘上的点,如果是,删除之。
[8] 计算特征属性尺度
对[7]得到的每个特征,计算特征尺度,对应的函数是calc_feature_scales。在[6.4]里,要计算subintvl,这是一个小数,实际上,每个图像的octave和intvl都是整数,但整数的intvlb并不能精确表示极值点,subintvl就是极值点所在的intvl的增量,这一点在[6.4]步骤是很清楚的。在feat的数据结构里,有一个ddata数据结构,记录suvintvl,在这里用到它了。语句:feat->scl = sigma * pow( 2.0, ddata->octv + intvl / intvls )中,signa是常量1.6,octv就是极值点所在的octave的序号,intvl是极值点所在的精确intvl,它很可能是带小数点的实数。ddata->scl_octv = sigma * pow( 2.0, intvl / intvls )就含义类似。这两个值在第3步,进行方向关联的时候,取关键点的邻域时候用得到。
5. 方向关联
每个关键点都跟关联一个常量的方向值,让特征描述子对图像旋转具有不变性。
[9] 计算特征方向
对应函数calc_feature_oris。
[9.1] 对于每个[8]中获取的特征,都要对它进行方向关联。也就是说,确定它所在的高斯金字塔图像的一定大小的邻域内的主方向。
[9.2]计算直方图函数是ori_hist,输入参数是,该特征在高斯图像金字塔所在图像的指针,该特征在这个图像上的行号和列号,直方图的桶数(bins),这里是常量36,该特征的邻域大小,注意,这里是根据[8]里计算出来的scl_octv乘以一个常量计算的,也就是说,基本上,intvl线性变化,最后一个参数是sigma,也是根据intvl线性变化。这个函数的计算是呆定的,是按照论文上写的做的。需要注意的是:源代码使用的是atan2函数,计算出来的角度,跟纵横坐标的符号相关,它的值域在[-pi, pi],这里做的处理是bin = cvRound( n * ( ori + CV_PI ) / PI2 ),ori+pi实在[0, 2pi],也就是说,所有的角度都可能,每个bin里都会可能有值。计算出直方图后,对直方图进行平滑操作。
[9.3] 将“好的”特征压栈:对于直方图的每个bin,如果它比前一个bin和后一个bin都大且超过一个阈值,那么,就复制当前的feat,将这个bin对应的方向值跟feat做关联,然后压栈。注意,在步骤[9.1]是从栈的前面pop一个feat进行处理,而这里是将特征push,这种用法效率比开两个栈要高一些。此外,这种计算方法,一个特征点可能生成多个新特征,只要它的方向值足够大。
6.图像局部描述子
这是SIFT算法的四个步骤的最后一个。它在源代码里对应函数compute_descriptors。
[10]计算局部描述子
[10.1计算描述直方图
[10.1.1]本文将特征点的近邻区分成4*4个区块,在每个区块里,有若干个像素点,这些像素点都具有不同的方向和幅值,要将这些像素点的方向和幅值计算直方图,直方图的bin是8个,也就是说,将每个bin对应45度。
[10.1.2]对于每个特征点,首先,计算radius,也就是它的计算描述直方图的邻域的半径。
[10.1.3]对于邻域内的每个点,要将这个点的当前坐标,转化到主方向上去。着一块在代码里对应的是:
c_rot = ( j * cos_t - i * sin_t ) / hist_width;
r_rot = ( j * sin_t + i * cos_t ) / hist_width;
rbin = r_rot + d / 2 - 0.5;
cbin = c_rot + d / 2 - 0.5;
if( rbin > -1.0 && rbin < d && cbin > -1.0 && cbin < d )
求c_rot和r_rot是比较费解的。从后面我们可以知道,r+i,表明i是像素纵坐标上的变量,c+j,表明j是像素横坐标上的变量。如果ori为0,那么,此时i,j不需要旋转。因为主方向不变。假设ori不是0,那么,就需要将现在的坐标系,旋转到ori的方向上去。假设有坐标为[i, j],r是极轴,sita是极角,那么,在极坐标下,i,j 点在新坐标系下面,就是x_new = r*cos(sita+ori) = r*cos(sita)cos(ori) - r*sin(sita)sin(ori) = j*cos(ori)-i*sin(ori),y_new = r*sin(sita+ori) = r*sin(sita)cos(ori) + r*cos(sita)sin(ori) = i*cos(oir) + j*sin(ori)。这个要根据图像坐标计算,也就是说,原点在左上角,x向右 为正,y向下为正。
为什么要除以hist_width呢?在转换之后的新坐标系下,要看c_rot和r_rot跟hist_width之间的关系。因为,计算区域是限制在4个hist_width的,也就是说,c_rot和r_row的值,应该在-2*hist_width和2*hist_width之间。
为什么要计算rbin和cbin呢?我们用下面方式推理:if条件,其实就是说-1 < rbin < d <==> -1 < rbin<4 <==> -1 < r_rot+2-0.5 < 4 <==> -2.5 < r_rot < 2.5。而 r_rot是新坐标值与hist_width的比值。if判断,其实就是判断r_rot是不是在指定的范围内,也就是几个hist_width。
如果if判断通过了,就是计算该像素的方向和幅值。计算方向之后,要将方向减去ori,因为此时计算出来的是当前坐标系下的方向,不是新坐标系下的方向。然后,要根据新坐标系下的坐标,计算权重。
[10.1.4]计算每个局部区域的直方图,对应函数是interp_hist_entry。在这个函数,根据计算出的rbin和cbin,算出来这个点,应该属于4个*4区块的那个区块,然后,再根据它的方向和幅值,将它累加到不同的bin里。这样,就得到了论文p15页图的右侧的小图,唯一对区别,就是这个图不是4*4,是2*2。
[10.2]将描述直方图转化描述子,对应函数hist_to_descr。这个就没什么好说啦,就是将4*4*8个值,逐次排列起来,然后做归一化,再转化成整数。
7. 目标识别的应用
这个是见仁见智,不一而足,可以在原文之上做更多的发挥。

转自:http://blog.sina.com.cn/s/blog_916b71bb0100upwx.html

  评论这张
 
阅读(516)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017