SIFT算法的Matlab实现

Posted 2022-10-29 sununs11

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

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

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

1. 构建尺度空间

定理1 对图像做 $\sigma ={\sigma }_1$ 的高斯平滑，再做一次 $\sigma ={\sigma }_2$ 的高斯平滑，等效于对原图像做一次 $\sigma =\sqrt{{\sigma }_1^2+{\sigma }_2^2}$ 的高斯平滑。

1.1 构建高斯金字塔

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

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

KaTeX parse error: No such environment: split at position 8: \\begin̲s̲p̲l̲i̲t̲̲ L(x,y,\\sigma)&…

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

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

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

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

KaTeX parse error: No such environment: split at position 8: \\begin̲s̲p̲l̲i̲t̲̲ sigma(1,n)&=\\s…

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-H2YGuKey-1585897571682)(http://7d9hqq.com1.z0.glb.clouddn.com/images/sift-implementation-in-matlab/pyr.png)]

如上图所示，在第一个octave中尺度为$k^3\cdot {\sigma }_0$的“最后”一张图像进行下采样得到第二个octave的第一张图像，尺度仍为$k^3\cdot {\sigma }_0=2\cdot {\sigma }_0$。

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

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

1.2 构建高斯差分金字塔

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

1.3 检测极值点

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-sklqn4D8-1585897571683)(http://7d9hqq.com1.z0.glb.clouddn.com/images/sift-implementation-in-matlab/extremum.png)]

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

2. 关键点精确定位

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

$$D(\mathbf{X})=D+{\frac{\partial D}{\partial \mathbf{X}}}^T\mathbf{X}+\frac{1}{2}{\mathbf{X}}^T\frac{{\partial }^{2}D}{\partial {\mathbf{X}}^2}\mathbf{X}$$

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

$$\hat{\mathbf{X}}=-{\frac{{\partial }^{2}D}{\partial {\mathbf{X}}^2}}^{-1}\frac{\partial D}{\partial \mathbf{X}}$$

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

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

精确极值点处函数值：

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

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

2.2 消除边缘响应

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

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

T r ( H ) 2 D e t ( H ) = ( r + 1 )