SIFT算法的Matlab实现

个人博客原文：http://www.sun11.me/blog/2016/sift-implementation-in-matlab/

这是一次作业，内容是给出两张图像，检测特征点和匹配特征点。要求不能用诸如OpenCV的现成特征点检测函数。于是就只能造轮子了，写了这个Matlab版的sift。（其实就是把c语言的opensift翻译成了matlab

以下是算法流程，其实网上的类似博文已经很多了，只不过我看的过程中也看得不很明白，只能对照着好几个看，所以干脆自己又写了一遍。下面的图均来自于参考资料中，然而参考资料的图也是来自于参考资料的参考资料中。

1. 构建尺度空间

定理1 对图像做 σ=σ1\sigma = \sigma_1 的高斯平滑，再做一次 σ=σ2\sigma = \sigma_2 的高斯平滑，等效于对原图像做一次 σ=σ21+σ22−−−−−−√\sigma = \sqrt{\sigma_1^2+\sigma_2^2} 的高斯平滑。

1.1 构建高斯金字塔

高斯卷积核是实现尺度变换的唯一线性核（Koenderink, 1984; Lindeberg, 1994）。

一幅图像的尺度空间被定义为对其做可变尺度的高斯卷积：

L(x,y,σ)其中G(x,y,σ)=G(x,y,σ)∗I(x,y)=12πσ2e−(x2+y2)/2σ2
\begin{split}
L(x,y,\sigma)&=G(x,y,\sigma) * I(x,y)\\
其中G(x,y,\sigma)&=\frac{1}{2\pi\sigma^2}e^{-(x^2+y^2)/2\sigma^2}
\end{split}

对于给定的彩色图像，转化为灰度图像，用不同大小的σ\sigma做高斯平滑（按照 3σ\sigma 准则，高斯核矩阵的大小设为 (6σ+1)⋅(6σ+1)(6\sigma+1)\cdot(6\sigma+1) ，并保证行和列为奇数），再此基础上将图像降采样得到不同大小的组(octave)，每组若干图像(interval)。详细描述如下：

为了得到更多的特征点，将图像扩大为原来的两倍。假设原图像已有 σ=0.5\sigma=0.5 的高斯平滑，而我们需要第一个octave的第一张图像的 σ=1.6\sigma=1.6 ，按照定理1，我们要对扩大两倍的图像做一次高斯平滑，σ=1.62−(0.5×2)2−−−−−−−−−−−−−√\sigma=\sqrt{1.6^2-(0.5\times2)^2} 。

上一个octave的图像的长度和宽度分别是下一个octave的图像的两倍。因此图像组数(octaves)可由图像大小决定，将其设为 log2(min(height,width))log_2(min(height,width)) − 2-\ 2 ，这样将使顶层octave图像的长度和宽度最小值在8像素左右。

设第m个octave的第n张图像相对于原始图像的参数σ\sigma为 sigma(m,n)sigma(m,n)，则sigma(1,1)=σ0=1.6sigma(1,1)=\sigma_0=1.6。每个octave有s+1张图像（即intervals），这样得到的高斯差分金字塔(DoG)每个octave将有s张图像，我们设s为3。为了满足在不同octave间尺度的连续性，并使 sigma(m,n)=sigma(m,n)= 2⋅sigma(m−1,n)2 \cdot sigma(m-1,n)，按照定理1，则：

sigma(1,n)sigma(m,n)=σ0⋅kn−1,其中k=21/s=σ0⋅2m−1⋅kn−1
\begin{split}
sigma(1,n)&=\sigma_0\cdot k^{n-1},其中k=2^{1/s}\\
sigma(m,n)&=\sigma_0\cdot 2^{m-1}\cdot k^{n-1}
\end{split}



如上图所示，在第一个octave中尺度为k3⋅σ0k^3\cdot \sigma_0的“最后”一张图像进行下采样得到第二个octave的第一张图像，尺度仍为k3⋅σ0=2⋅σ0k^3\cdot \sigma_0=2\cdot \sigma_0。

但实际上我们需要做出更多不同尺度的高斯平滑图像，这是因为在后续高斯差分金字塔的极值检测中，需要前后两级尺度都存在图像。如图中红框所示，高斯差分金字塔中每个octave有s幅图像，则需要高斯金字塔中每个octave包含s+3幅图像。其中第s+1幅图像用作下一个octave第一幅图像的降采样。

具体实现中并未对单幅图像多次进行高斯平滑，而是由上一幅图像进行高斯平滑得到下一幅图像并迭代之，按照定理1计算σ\sigma。

1.2 构建高斯差分金字塔

对两幅高斯金字塔的图像作差。

1.3 检测极值点

如上图，与前后两幅图像及自身的共26个邻域像素点比较灰度值检测极值。

2. 关键点精确定位

检测到的极值点是离散的，通过三元二次函数拟合来精确确定关键点的位置和尺度，达到亚像素精度。以某关键点为中心的尺度空间函数 D(x,y,intvl)D(x,y,intvl) 的二次泰勒展开式为：

D(X)=D+∂D∂XTX+12XT∂2D∂X2X
D(\mathbf{X}) = D + \frac{\partial D}{\partial\mathbf{X}}^T\mathbf{X} + \frac{1}{2}\mathbf{X}^T\frac{\partial ^2D}{\partial\mathbf{X}^2}\mathbf{X}

其中等号右边第一个D为某关键点处的灰度值， X=(x,y,intvl)T\mathbf{X}=(x,y,intvl)^T 是以此点为中心的偏移量，由于 D(X)D(\mathbf{X}) 是离散的，其导数用差分法求得。令 D(X)D(\mathbf{X}) 导数为零，得到精确极值位置的偏移量为：

X^=−∂2D∂X2−1∂D∂X
\mathbf{\hat{X}}=-\frac{\partial ^2D}{\partial \mathbf{X}^2}^{-1}\frac{\partial D}{\partial\mathbf{X}}

若X^\mathbf{\hat{X}}在任意一个维度大于0.5，说明极值点精确位置距离另一个点更近，应该将关键点定位于更近的那个位置。定位到新点后再进行相同操作，若迭代5次位置仍不收敛，则不认为此点为关键点。设定图像边缘img_border，若关键点落在图像边缘区域（以img_border为宽度的矩形外框）也不认为此点为关键点。

2.1 去除低反差(low contrast)的点

精确极值点处函数值：

D(X^)=D+12∂D∂XTX^
D(\mathbf{\hat{X}}) = D + \frac{1}{2}\frac{\partial D}{\partial\mathbf{X}}^T\mathbf{\hat{X}}

若 |D(X^)|<0.04/s|D(\mathbf{\hat{X}})|<0.04/s ，同样不认为此点是极值点。在此过程中保存极值点的数据ddata，为特征的构建做准备。计算出σ_octv\sigma\_octv，即位于一个相同的octave内的尺度，某个octave内第n张图像的 σ_octv=σ0⋅kintvl−1\sigma\_octv= \sigma_0\cdot k^{intvl-1} ，此处intvl为精确定位后的intvl。

2.2 消除边缘响应

高斯差分函数有较强的边缘响应，对于比较像边缘的点应该去除掉。这样的点的特征为在某个方向有较大主曲率，而在垂直的方向主曲率很小。

设r为大主曲率与小主曲率的比值，H为关键点处的Hessian矩阵，则有（具体推导可见Lowe的论文）：

Tr(H)2Det(H)=(r+1)2r\frac{Tr(H)^2}{Det(H)}=\frac{(r+1)^2}{r}

若满足：

Tr(H)2Det(H)<(rt+1)2rt,其中rt为一阈值，设为10
\frac{Tr(H)^2}{Det(H)}<\frac{(r_t+1)^2}{r_t},其中r_t为一阈值，设为10

说明此处r较小，认为此关键点不位于边缘，否则则去除该点。

3. 方向指定

根据关键点的局部特性为每个关键点指定一个方向，可以具备旋转不变性。关键点局部特性在检测到关键点的高斯差分金字塔图像临近的高斯金字塔图像中计算。在关键点3σ邻域窗口计算梯度和方向分布，计算方式如下：

m(x,y)θ(x,y)=[L(x+1,y)−L(x−1,y)]2+[L(x,y+1)−L(x,y−1)]2−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−√=tan−1{[L(x,y+1)−L(x,y−1)]/[L(x+1,y)−L(x−1,y)]}
\begin{split}
m(x,y)&=\sqrt{[L(x+1,y)-L(x-1,y)]^2+[L(x,y+1)-L(x,y-1)]^2}\\
\theta(x,y)&=tan^{-1}\{[L(x,y+1)-L(x,y-1)]/[L(x+1,y)-L(x-1,y)]\}
\end{split}

此处的x正方向向右，y正方向向上。其中L为关键点在上述精确定位后所处尺度的灰度值，m(x,y)为梯度的幅值，θ(x,y)\theta(x,y)为关键点处梯度方向的弧度（范围为(−π,π](-\pi,\pi]）。将360度的方向划分为36个区域(bins),第一个区域的范围是[35π36,37π36)[\frac{35\pi}{36},\frac{37\pi}{36}),按逆时针方向依次划分。对m(x,y)按 σ=1.5σ_octv\sigma=1.5\sigma\_octv 的高斯分布，在 3σ=3⋅1.5σ_octv3\sigma=3\cdot1.5\sigma\_octv 的邻域窗口加权计算，得到36个方向的直方图。然后对直方图进行两次平滑处理，即按0.25,0.5,0.25的大小对每3个连续的bin加权两次。

直方图最大值的方向代表该关键点的主方向，对于其他峰值，若大于或等于主方向值的80%，则再分配一个方向。所以对于一个关键点，可能会有多个对应的方向，将带有方向的关键点定义为feature，则一个关键点可能对应多个feature。由于第一个octave是双倍大小的图像，feature的坐标和尺度应转换到原始图像所在的octave处理。最后用抛物插值精确定位feature的方向。

对于x为-1,0,1，y为l,c,r的三个点来说，抛物插值得到极值点的x为：

0.5⋅l−rl−2c+r
0.5\cdot\frac{l-r}{l-2c+r}

4. 关键点描述子

上一步已得到具有主方向的关键点，即feature，下一步是对feature的邻域进行采样，形成对该局部图像的描述，然后可用某种度量方法对描述进行匹配。

Lowe提出的sift描述子是一个 4×4×8=1284\times4\times8=128 维的向量。描述子的数学形式可定义为 h(x,y,θ)h(x,y,\theta) ，其中的x,y代表 4×4=164\times4=16 个图像区域的位置，θ\theta即梯度方向，只能取8个值。h(x,y,θ)h(x,y,\theta)的值就是在(x,y)代表的图像区域计算得到的在θ\theta方向的梯度大小。

4.1 描述子采样区域

这16个图像区域中的每一个区域均为 3σ_octv3\sigma\_octv 像素，因此16个区域的半边长为 4×3σ_octv/24\times3\sigma\_octv/2 ，考虑到后续操作需要三线性插值，采样区域半边长设为 (4+1)×3σ_octv/2(4+1)\times3\sigma\_octv/2 ，又由于旋转操作，这个值需要乘以2√\sqrt{2}，得到 radius=(4+1)×radius=(4+1)\times 2√×3σ_octv/2\sqrt{2}\times3\sigma\_octv/2 。

如下图所示，图中 m=3,m=3, Bp=4,σ=σ_octvB_p=4,\sigma=\sigma\_octv 。

4.2 旋转至主方向

为了使描述符具有旋转不变性，将坐标轴旋转至关键点主方向。设i，j分别为采样点相对关键点的行偏移量和列偏移量，i = -radius:radius，j = -radius:radius，关键点左上角i和j均为负数。关键点主方向为θ\theta，范围是(−π,π](-\pi,\pi]。

定理2 在右手平面直角坐标系中，向量(x,y)逆时针旋转θ\theta，得到的向量(x’,y’)为：

[x′y′]=[cosθsinθ−sinθcosθ][xy]
\begin{bmatrix}
x'\\
y'
\end{bmatrix}
=\begin{bmatrix}
cos\theta & -sin\theta \\
sin\theta & cos\theta
\end{bmatrix}
\begin{bmatrix}
x\\
y
\end{bmatrix}



在左图以关键点为中心建立右手平面直角坐标系o-uv，u的正方向与左图x\mathbf{x}方向相同，v的正方向与左图y\mathbf{y}方向相反。左图中x\mathbf{x}=(1,0),y\mathbf{y}=(0,-1)，将x\mathbf{x},y\mathbf{y}代入定理2的公式，得到右图中 x′=(cosθ,sinθ),\mathbf{x'}=(cos\theta,sin\theta), y′=(sinθ,−cosθ)\mathbf{y'}=(sin\theta,-cos\theta) 。其中θ\theta为左图坐标系旋转到右图坐标系的角度，在上图中为一负数。设图像中有一点按左图的x\mathbf{x},y\mathbf{y}可表示为 j⋅x+i⋅yj \cdot\mathbf{x} + i\cdot\mathbf{y} ，按右图中的x′\mathbf{x'},y′\mathbf{y'}可表示为 j′⋅x′+i′⋅y′j'\cdot\mathbf{x'}+i'\cdot\mathbf{y'} ，则有：

j[10]+i[0−1]=j′[cosθsinθ]+i′[sinθ−cosθ]
j\begin{bmatrix}
1\\
0
\end{bmatrix}
+i\begin{bmatrix}
0\\
-1
\end{bmatrix}
=j'\begin{bmatrix}
cos\theta\\
sin\theta
\end{bmatrix}
+i'\begin{bmatrix}
sin\theta\\
-cos\theta
\end{bmatrix}

解得 j′=j⋅cosθ−i⋅sinθ,j'=j\cdot cosθ-i\cdot sinθ, i′=j⋅sinθ+i⋅cosθi'=j\cdot sinθ+i\cdot cosθ 。

得到新的行、列偏移量后，将 3σ_octv3\sigma\_octv 设为单位长度，并将中心转移至最左上角区域中心，得到新的坐标r_bin和c_bin。对梯度方向弧度值减去主方向弧度，并设 2π8\frac{2\pi}{8} 为一个单位，得到o_bin。

采样点的梯度幅值按照 σ=0.5⋅4⋅3σ_octv\sigma=0.5\cdot4\cdot3\sigma\_octv （即16个区域边长的一半）的高斯函数加权：

w=m(a+j,b+i)⋅e−j′2+i′22σ2
w=m(a+j,b+i)\cdot e^{\frac{-j'^2+i'^2}{2\sigma^2}}

其中a，b为关键点在高斯金字塔图像中的位置坐标。

4.3 三线性插值

上述过程中构造了一个三维的bin空间，如4.1中右图所示，维度包括r_bin，c_bin和o_bin。注意最上层格子和最底层格子是相连的，因为0度等于360度。所有带有三维坐标的梯度幅值都将分配到三维格子里。

为了减少一个梯度幅值从一个格子漂移(shift)到另一个格子引起的描述子突变，需要对梯度值做三线性插值。也就是根据三维坐标计算距离周围格子的距离，按距离的倒数计算权重，将梯度幅值按权重分配到临近的格子里。



某点在三维bin空间的坐标为(r_bin，c_bin，o_bin)(r\_bin，c\_bin，o\_bin)，求出r=⌊r_bin⌋,r=\lfloor r\_bin\rfloor, c=⌊c_bin⌋,c=\lfloor c\_bin\rfloor, o=⌊o_bin⌋,o=\lfloor o\_bin\rfloor, dr=r_bin−r,dr=r\_bin-r, dc=c_bin−c,dc=c\_bin-c, do=o_bin−odo=o\_bin-o，它的梯度幅值最多可能分配到周围的8个格子中。计算公式如下：

=weightedValue[r+i+1][c+j+1][((o+k) mod 8)+1] w⋅dri⋅(1−dr)1−i⋅dcj⋅(1−dc)1−j⋅dok⋅(1−do)1−k
\begin{split}
&weightedValue[r+i+1][c+j+1][((o+k)\ mod\ 8)+1]\\
=&\ w\cdot dr^i\cdot(1-dr)^{1-i}\cdot dc^j\cdot(1-dc)^{1-j}\cdot do^k\cdot(1-do)^{1-k}
\end{split}

其中i，j，k均可取0或1，weightedValue下标加1的目的是使下标从1开始。

为简化计算，可改为：

=weightedValue[r+i+1][c+j+1][((o+k) mod 8)+1] w⋅(0.5+(dr−0.5)(2i−1))⋅(0.5+(dc−0.5)(2j−1))⋅ (0.5+(do−0.5)(2k−1))
\begin{split}
&weightedValue[r+i+1][c+j+1][((o+k)\ mod\ 8)+1]\\
=&\ w\cdot(0.5+(dr-0.5)(2i-1))\cdot(0.5+(dc-0.5)(2j-1))\cdot\\
&\ (0.5+(do-0.5)(2k-1))
\end{split}

4.4 生成描述子

将上述直方图数组按顺序排列可转换为一个128维的向量。

为了减少光照变化的影响，对该向量进行归一化处理。非线性光照变化仍可能导致梯度幅值的较大变化，然而影响梯度方向的可能性较小。因此对于超过阈值0.2的梯度幅值设为0.2，然后再进行一次归一化。最后将描述子按照对应高斯金字塔图像的尺度大小排序。

5. 匹配

描述子向量已经归一化，所以可直接用向量之间的夹角进行匹配，相当于球面距离。图像A 的描述子匹配图像B最近的两个描述子点积之比小于0.6，则认为匹配成功。

6. 一些废话

6.1 性能优化

因为用的是Matlab，所以不注重性能。然而又不得不注重性能，因为第一次跑通程序的时候跑了一晚上都没跑完一半！也就是一张图片的描述子都没算完。后来发现是因为在运行次数最多的for循环(描述子计算中的梯度计算)里用到了cell数组。把对这个cell数组的查询操作提到两重循环前以后，这个程序好像跑了半个小时左右跑出结果了。然而还是太慢，于是我又用Matlab的计时分析工具分析了程序最耗时的地方:

把cell数组的查询尽可能减少

充分利用Matlab的向量操作

一些没用的参数给去掉了（如计算梯度时的三个返回值合并到了一个二维数组）

一个三维数组折叠成了一维的（hist）

程序里用了很多全局变量，是因为我把函数分成了文件而不是放在一个文件，为了节省点内存（以及方便）只能这么做（虽然据说Matlab在不改变变量的情况下函数传值等于引用，然而我并不清楚究竟是怎样的）。把for循环换成parfor的时候提示，parfor里似乎不推荐用全局变量，而且实际运行的时候全局变量似乎也会影响性能，于是我把全局变量复制成了局部的再放进parfor里。

我还发现一个奇葩的问题，在运行次数最多的计算梯度的函数里用zeros(1,2)创建一个数组竟然也耗时非常多，改成[0 0]就好了。

经过这些修改后，在开启parallel pool的情况下运行时间缩短到了7分钟左右。（然而Lowe的C语言版本只要十几秒

6.2 运行结果

这次作业老师给的是两张768x1024的图片，分别检测到5288和4798个特征点，最后匹配了906对点。用Lowe的siftDemoV4跑出来的结果是1252对匹配。

这个程序的参数基本都是参照opensift，但最后的匹配用的是Lowe的方案。Lowe的实现毕竟不太一样，运行的结果和opensift有一些差异。以下是匹配

siftDemoV4.zip

里的

scene.pgm

和

book.pgm

的结果：

Item	siftDemoV4	opensift	sw-sift
scene.pgm	1021	746	766
book.pgm	882	740	741
Matched	98	84	58

sw-sift和opensift的区别主要是在高斯平滑和匹配算法上。opensift的高斯平滑用的是OpenCV的CVSmooth函数，匹配用的是欧式距离（而且把描述子乘以512从double类型转成了int）。和opensift相比，sw-sift检测到的特征点数量很接近，但是匹配数量较少，所以可改进的地方主要是匹配算法（然而我不想改了==）。另外，我发现高斯平滑的核矩阵大小对结果有很大影响，根据3σ3\sigma准则它的宽度应该是 (6σ+1)⋅(6σ+1)(6\sigma+1)\cdot(6\sigma+1) ，然而有人设成 (3σ+1)⋅(3σ+1)(3\sigma+1)\cdot(3\sigma+1) 却取得了更多特征点，因此调整这个参数再用其它参数限制错误数量或许可以得到更好的结果。

6.3 源代码

代码发布在github: sw-sift，请注意sift是有专利的。

参考资料

David G. Lowe, “Distinctive image features from scale-invariant keypoints,” International Journal of Computer Vision, 60, 2 (2004), pp. 91-110. [PDF] [CODE]

Rob Hess, OpenSIFT源码: https://github.com/robwhess/opensift

zddhub, SIFT算法详解: /article/2226726.html

Rachel Zhang, SIFT特征提取分析: /article/1820926.html

JiePro, SIFT算法：特征描述子: /article/7125337.html