Local Invariant Feature Detectors

Table of Contents

1 图像局部不变性特征

什么是局部特征?A local feature is an image pattern which differs from its immediate neighborhood.

参考:
图像局部不变性特征与描述,王永明/王贵锦(本文主要摘自该书)
Local Invariant Feature Detectors - A Survey, 2008

1.1 特征的不变性

特征的不变性有很多,比如旋转不变性(物体旋转后也能检测出来),亮度不变性(物体变暗或变亮都能检测出来),尺度不变性(物体缩放后也能检测出来),仿射不变性(后文将介绍)等等。

1.1.1 仿射变换

仿射变换(Affine transformation)的定义为:

\begin{equation} \vec{y} = A \vec{x} + \vec{b} \tag{Affine transformation} \end{equation}

用于图像处理时, \(\vec{x}\) 是原图像像素, \(\vec{y}\) 是仿射变换后的图像像素, \(A\) 是一个 \(2 \times 2\) 的矩阵。从定义可知, 仿射变换可表示为两个步骤:乘以一个矩阵(即线性变换)接着再加上一个向量(即平移)。

1 (摘自Linear transformations)是乘以一个矩阵(即线性变换)的例子。

img_linear_transformations.jpg

Figure 1: 乘以一个矩阵(即线性变换)所对应图像变换的例子

图像处理中,仿射变换还可以写为:
\[\begin{bmatrix} \vec{y} \\ 1 \end{bmatrix} = \underbrace{\begin{bmatrix} A & \vec{b} \\ \boldsymbol{0} & 1 \end{bmatrix}}_{3 \times 3} \begin{bmatrix} \vec{x} \\ 1 \end{bmatrix}\]

从而,仿射变换可以 \(3 \times 3\) 的矩阵表示。图 2 (摘自Affine transformation)是仿射变换的一些例子。

img_affine_transformation.png

Figure 2: 仿射变换的例子

由于仿射变换的 \(3 \times 3\) 矩阵中最后一行固定为 \(\begin{bmatrix} 0 & 0 & 1 \end{bmatrix}\) ,所以可以去掉这行,故仿射变换也可以用 \(2 \times 3\) 的矩阵表示。

参考:
http://www.opencv.org.cn/opencvdoc/2.3.2/html/doc/tutorials/imgproc/imgtrans/warp_affine/warp_affine.html

2 尺度空间理论简介

图像尺度空间(Scale space)分析技术是一种性能突出的多尺度分析(Multiple-scale analysis)技术,它在图像处理中占据非常重要的地位。

参考:
Tony Lindeberg. (1998), "Feature Detection with Automatic Scale Selection" (pdf)

2.1 尺度空间定义

图像的尺度空间记为 \(L(x,y,\sigma)\) ,它的定义是二维高斯函数 \(G(x,y,\sigma) = \frac{1}{2 \pi \sigma^2} e^{- \frac{x^2 + y^2}{2 \sigma^2}}\) (注:有时把 \(\sigma^2\) 记为 \(t\) )和图像 \(I(x,y)\) 的卷积。 即:
\[L(x,y,\sigma) = G(x,y,\sigma) \star I(x,y)\]
其实,上式就是高斯模糊(Gaussian blur)运算。当参数 \(\sigma\) (称为尺度空间因子,或尺度因子)取不同值时,可以得到一系列模糊程度不同的图像( \(\sigma\) 越大图像越模糊),这就是图像的尺度空间。如图 3 所示:

img_scale_space.jpg

Figure 3: 上左(原图),上右(高斯模糊 \(\sigma=4\) ),下左(高斯模糊 \(\sigma=10\) ),下右(高斯模糊 \(\sigma=20\) )

在不同尺度上,可以检测出不同的特征。大尺度图像上可以检测概貌特征,小尺度图像上可以检测细节特征。

2.2 尺度空间理论应用

尺度空间理论不是一个具体的图像处理算法,而是一个通用的“图像处理框架”。 比如,可以在这个框架下检测角点(Corner)、斑点(Blob)等等。Harris-Laplace角点检测、SIFT斑点检测、SURF斑点检测等算法都是尺度空间理论的应用。

2.3 特征检测和定位


图像特征往往在粗糙的尺度上被检测到,此时位置信息未必是最准确的。因此, 通常图像的尺度分析包含两个阶段:首先在粗尺度上进行特征检测,然后再在细尺度上进行精确定位。 如图 4 (摘自:Tony Lindeberg. (1998), "Feature Detection with Automatic Scale Selection")所示。

img_scale_space_localization.jpg

Figure 4: 尺度分析中角点检测实例。第1列为原图像,第2列为不同尺度上检测的结果,第3列是精确定位后的结果。

3 Corner检测

角点(Corner)对应于物体的拐角,如道路的十字路口、丁字路口等。

角点检测方法主要分为两类:基于图像“边缘”的方法和基于图像“灰度”的方法。本文介绍的Harris角点检测和SUSAN角点检测都是后者,即基于图像“灰度”的方法。

3.1 Harris角点检测

Harris角点检测的基本思路如下: 考虑一个局部的小区域(或称小窗口),如果在各个方向上移动这个特定的小窗口,窗口内图像的灰度发生了较大的变化,那么,就认为在窗口内遇到了角点,如图 5 左图所示。 如果这个特定的窗口在图像各个方向上移动时,窗口内图像的灰度没有发生变化,那么窗口内就不存在角点,如图 5 中图所示;如果窗口在某一个方向移动时,窗口内图像的灰度发生了较大的变化,而在另一些方向上没有发生变化,那么,窗口内的图像可能就是一条直线的线段,如图 5 右图所示。

img_harris_window.jpg

Figure 5: 左图(窗口移动与角点);中图(窗口移动与平面);右图(窗口移动与线段)

对于图像 \(I(x,y)\) ,当在点 \((x,y)\) 处平移 \((\Delta x, \Delta y)\) 后的自相似性,可以通过下面自相关函数给出:
\[S_W(\Delta x,\Delta y) = \sum_{x_i \in W} \sum_{y_i \in W} w(x_i, y_i)(\underbrace{I(x_i + \Delta x, y_i + \Delta y)}_{\text{平移后的图像}} - I(x_i, y_i))^2\]
式中, \(W\) 是以点 \((x,y)\) 为中心的窗口(即图 5 中的小窗口), \(w(x_i,y_i)\) 为加权函数,它既可是常数,也可以是高斯加权函数 \(e^{- \frac{(x_i -x)^2 + (y_i -y)^2}{2\sigma^2}}\) (靠近窗口中心的权重比较大,远离窗口中心的权重比较小),如图 6 所示。

img_harris_weight_func.jpg

Figure 6: 左图(常量加权函数),右图(高斯加权函数)

当 \((\Delta x, \Delta y)\) 取不同值时(即相当于在各个方向上移动图 5 所示小窗口),如果 \(S_W(\Delta x,\Delta y)\) 都有较大的变化,则可以认为点 \((x,y)\) 是一个角点。这就是Harris角点检测的基本思路。

3.1.1 两个例子

考虑图 7 所示例子,我们要检测红色框中心点是否为角点。

img_harris_example.png

Figure 7: 实例1:检测红色框中心点是否为角点

先考虑两个特例,假设 \(\Delta y=0\) , \(S_W(\Delta x,\Delta y)\) 的响应如图 8 所示,假设 \(\Delta x=0\) , \(S_W(\Delta x,\Delta y)\) 的响应如图 9 所示。

img_harris_example_y0.gif

Figure 8: 设 \(\Delta y=0\) ,当 \(\Delta x\) 变化时 \(S_W(\Delta x,\Delta y)\) 变化情况(变化较大)

img_harris_example_x0.gif

Figure 9: 设 \(\Delta x=0\) ,当 \(\Delta y\) 变化时 \(S_W(\Delta x,\Delta y)\) 变化情况(变化较小)

对图 7 红色框中心点,考虑 \(\Delta x, \Delta y\) 所有变化情况时,\(S_W(\Delta x,\Delta y)\) 的响应如图 10 所示(峡谷形)。

img_harris_example_xy.png

Figure 10: 7 中当 \(\Delta x, \Delta y\) 变化时,\(S_W(\Delta x,\Delta y)\) 的响应图

从图 10 中明显可看到,当 \(\Delta x, \Delta y\) 取不同值时(相当于从各个方向移动小窗口),\(S_W(\Delta x,\Delta y)\) 不一定都有较大的变化,所以推断图 7 红色框中心点不是角点。

我们再看一个例子。考虑图 11 所示例子,我们要检测红色框中心点是否为角点。

img_harris_example2.png

Figure 11: 实例2:检测红色框中心点是否为角点

考虑 \(\Delta x, \Delta y\) 所有变化情况时,\(S_W(\Delta x,\Delta y)\) 的响应如图 12 所示(锥形)。

img_harris_example2_xy.png

Figure 12: 11 中当 \(\Delta x, \Delta y\) 变化时,\(S_W(\Delta x,\Delta y)\) 的响应图

从图 12 中明显可看到,当 \(\Delta x, \Delta y\) 取不同值时(相当于从各个方向移动小窗口),\(S_W(\Delta x,\Delta y)\) 都有较大的变化,所以推断图 7 红色框中心点是角点。

注:本节介绍的两个例子摘自:https://dsp.stackexchange.com/questions/3336/mathematics-of-harris-corner-point-detection

3.1.2 数学推导

根据二元函数泰勒展开,对图像 \(I(x,y)\) 在平移 \((\Delta x, \Delta y)\) 后的图像进行一阶近似( \(I_x,I_y\) 是图像 \(I(x,y)\) 的偏导数):
\[I(x_i + \Delta x, y_i + \Delta y) \approx I(x_i,y_i)+I_x(x_i,y_i) \Delta x+ I_y(x_i,y_i) \Delta y\]
把上式代入自相关函数 \(S_W(\Delta x,\Delta y)\) 中,从而有:
\[\begin{aligned} S_W(\Delta x,\Delta y) & \approx \sum_{x_i \in W} \sum_{y_i \in W} w(x_i, y_i) \left( I_x(x_i,y_i) \Delta x + I_y(x_i,y_i) \Delta y \right)^2 \\ &= \sum_{x_i \in W} \sum_{y_i \in W} w(x_i, y_i) \left( I_x(x_i,y_i)^2 \Delta x^2 + I_y(x_i,y_i)^2 \Delta y^2 + 2 I_x(x_i,y_i) I_y(x_i,y_i) \Delta x \Delta y \right) \\ &= \sum_{x_i \in W} \sum_{y_i \in W} w(x_i, y_i) \left[ I_x(x_i,y_i)^2 \Delta x + I_x(x_i,y_i) I_y(x_i,y_i) \Delta y, \; I_y(x_i,y_i)^2 \Delta y + I_x(x_i,y_i) I_y(x_i,y_i) \Delta x \right] \cdot \begin{bmatrix} \Delta x \\ \Delta y \\ \end{bmatrix} \\ &= \sum_{x_i \in W} \sum_{y_i \in W} w(x_i, y_i) \left[\Delta x, \Delta y \right] \cdot \begin{bmatrix}I_x(x_i,y_i)^2 & I_x(x_i,y_i)I_y(x_i,y_i) \\ I_x(x_i,y_i)I_y x_i,y_i) & I_y(x_i,y_i)^2\end{bmatrix} \cdot \begin{bmatrix} \Delta x \\ \Delta y \\ \end{bmatrix} \\ & \overset{\text{记为}}{=} \left[ \Delta x, \Delta y \right] M(x,y) \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix} \end{aligned}\]
其中:

\begin{align*} M(x,y) &= \sum_{x_i \in W} \sum_{y_i \in W} w(x_i, y_i) \begin{bmatrix}I_x(x_i, y_i)^2&I_x(x_i, y_i)I_y(x_i, y_i) \\ I_x(x_i, y_i)I_y(x_i, y_i)&I_y(x_i, y_i)^2\end{bmatrix} \\ &= \begin{bmatrix} \sum\limits_{x_i \in W} \sum\limits_{y_i \in W} w(x_i, y_i) I_x(x_i, y_i)^2 & \sum\limits_{x_i \in W} \sum\limits_{y_i \in W} w(x_i, y_i) I_x(x_i, y_i)I_y(x_i, y_i) \\ \sum\limits_{x_i \in W} \sum\limits_{y_i \in W} w(x_i, y_i) I_x(x_i, y_i)I_y(x_i, y_i) & \sum\limits_{x_i \in W} \sum\limits_{y_i \in W} w(x_i, y_i) I_y(x_i, y_i)^2\end{bmatrix} \\ & \overset{\text{记为}}{=} \begin{bmatrix}A&C\\C&B\end{bmatrix} \end{align*}

\(M(x,y)\) 只与图像中当前考查的点 \((x,y)\) 相关,它可以直接从图像中计算出来。通过上面的推导,可知图像 \(I(x,y)\) 在点 \((x,y)\) 处平移 \((\Delta x, \Delta y)\) 后的自相关函数可以近似为二次项函数:
\[\begin{aligned} S_W(\Delta x,\Delta y) & \approx \left[ \Delta x, \Delta y \right] M(x,y) \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix} \\ &= \left[ \Delta x, \Delta y \right] \begin{bmatrix}A&C\\C&B\end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix} \\ &= A \Delta x^2 + 2C \Delta x \Delta y + B \Delta y^2 \end{aligned}\]

由线性代数(可参考《工程数学线性代数(第五版),同济大学版》5.7节 正定二次型)中的知识可知: \(A \Delta x^2 + 2C \Delta x \Delta y + B \Delta y^2 = c\) (其中 \(c\) 为大于零的常数)的图形是以原点为中心的椭圆。椭圆的两个半轴由矩阵 \(\begin{bmatrix}A&C\\C&B\end{bmatrix}\) (它是实对称矩阵,一定有两个特征值)的特征值 \(\lambda_1, \lambda_2\) 决定的,如图 13 所示( \(\lambda_{\text{max}}, \lambda_{\text{min}}\) 分别是两个特征值中较大的值和较小的值)。当把 \(c\) 看作任意常数时,则是一族椭圆。这族椭圆随着常数 \(c \to 0\) 而收缩到原点,如图 14 所示。

img_harris_ellipse.jpg

Figure 13: \(A \Delta x^2 + 2C \Delta x \Delta y + B \Delta y^2=c\) 所表示的椭圆

img_harris_ellipse_c1_c2.jpg

Figure 14: 椭圆 \(A \Delta x^2 + 2C \Delta x \Delta y + B \Delta y^2=c\) 随着常数 \(c \to 0\) 而收缩到原点,这里 \(c_1 < c_2\)

前面介绍过,判断点 \((x,y)\) 是否为角点的办法是从点 \((x,y)\) 开始在各个方向(即向量 \((\Delta x, \Delta y)\) )上移动特定的小窗口 \(W\) ,如果 \(S_W\) (这里就是 \(c\) )都有较大的变化则认为点 \((x,y)\) 就是角点。这个角点判断条件可以换一种说法来表述, 假设 \(c\) 值固定,从各个方向移动小窗口 \(W\) 使 \(A \Delta x^2 + 2C \Delta x \Delta y + B \Delta y^2\) 达到固定响应值 \(c\) 时,如果需要移动的位移很小则可以判断点 \((x,y)\) 是角点。 从图 13 可知,当 \(\lambda_{\text{max}}, \lambda_{\text{min}}\) 都很大时( \(c\) 是固定的),椭圆两个半轴的长度就很小,而椭圆上的点 \((\Delta x, \Delta y)\) 代表的是达到响应 \(c\) 时小窗口 \(W\) 的位移,椭圆很小意味着位移很小。 从而可以得到结论:当矩阵 \(\begin{bmatrix}A&C\\C&B\end{bmatrix}\) 的特征值 \(\lambda_{\text{max}}, \lambda_{\text{min}}\) 都很大时,可以认为 \((x,y)\) 是角点。

3.1.3 矩阵M(x,y)特征值与图像区域类型的关系

前面介绍了如何通过矩阵 \(M(x,y)\) 特征值来判断图像区域是否为角点。还可以通过矩阵 \(M(x,y)\) 特征值来判断是否为平面和直线。总结如下:

矩阵 \(M(x,y)\) 的特征值与图像中不同类型(直线/平面/角点)的区域有下面关系:
(1) 图像中的直线。一个特征值大,另一个特征值小, \(\lambda_1 \gg \lambda_2\) 或 \(\lambda_2 \gg \lambda_1\) 。 自相关函数值在某一方向上大,在其他方向上小。
(2) 图像中的平面。两个特征值都小,且近似相等; 自相关函数数值在各个方向上都小(即 \((\Delta x, \Delta y)\) 取不同值时, \(S_W\) 的变化都很小)。
(3) 图像中的角点。两个特征值都大,且近似相等; 自相关函数值在所有方向上都大(即 \((\Delta x, \Delta y)\) 取不同值时, \(S_W\) 都有较大的变化)。

img_harris_ eigenvalue.jpg

Figure 15: 矩阵M(x,y)特征值与图像区域类型的关系

3.1.3.1 矩阵M的其它名字:"Structure Tensor" or "Second Moment Matrix"

矩阵M:
\[M(x,y) = \begin{bmatrix} \sum\limits_{x_i \in W} \sum\limits_{y_i \in W} w(x_i, y_i) I_x(x_i, y_i)^2 & \sum\limits_{x_i \in W} \sum\limits_{y_i \in W} w(x_i, y_i) I_x(x_i, y_i)I_y(x_i, y_i) \\ \sum\limits_{x_i \in W} \sum\limits_{y_i \in W} w(x_i, y_i) I_x(x_i, y_i)I_y(x_i, y_i) & \sum\limits_{x_i \in W} \sum\limits_{y_i \in W} w(x_i, y_i) I_y(x_i, y_i)^2\end{bmatrix}\]
它还有其它名字:"Structure Tensor" or "Second Moment Matrix"。

参考:https://en.wikipedia.org/wiki/Structure_tensor

3.1.4 Harris角点响应(避免直接求特征值)

求解矩阵特征值的费时的操作,Harris给出了不直接求特征值来判断角点的方法。Harris指出可以通过计算一个角点响应值 \(R\) 来判断角点。 \(R\) 的定义为:

\begin{equation} R= \det M - k(\text{trace} M)^2 \tag{Harris角点响应定义} \end{equation}

上式中, \(\det M\) 是矩阵 \(M=\begin{bmatrix}A&C\\C&B\end{bmatrix}\) 的行列式; \(\text{trace} M\) 是矩阵 \(M\) 的直迹; \(k\) 为经验常数,一般取0.04~0.06。事实上,特征是隐含在 \(\det M\) 和 \(\text{trace} M\) 中,因为(由矩阵相关知识可得):
\[\begin{aligned} \det M &= \lambda_1 \lambda_2 = AB-C^2 \\ \text{trace} M &= \lambda_1 + \lambda_2 = A+B \end{aligned}\]

Harris给出了下面结论:
(1) 当 \(|R|\) 很小时,这时 \(\lambda_1, \lambda_2\) 都很小,相应图像区域为平面;
(2) 当 \(R < 0\) 时,这时 \(\lambda_1 \gg \lambda_2\) 或 \(\lambda_2 \gg \lambda_1\) ,相应图像区域为直线;
(3) 当 \(R\) 比较大时,这时 \(\lambda_1, \lambda_2\) 都很大,且近似相等,相应图像区域为角点。

3.1.5 Harris角点检测算法总结

根据上述讨论,可以将Harris图像角点检测算法归纳如下,共分以下五步:
(1) 计算图像 \(I(x,y)\) 在 \(X\) 和 \(Y\) 两个方向的梯度 \(I_x, I_y\) 。
\[I_x=\frac{\partial I}{\partial x}=I \star (-1\ 0\ 1), I_y =\frac{\partial I}{\partial y}=I \star (-1\ 0\ 1)^T\]
(2) 计算图像两个方向梯度的乘积。
\[I_x^2 = I_x \cdot I_x, \; I_y^2 = I_y \cdot I_y, \; I_{xy} = I_x \cdot I_y\]
(3) 使用高斯函数对 \(I_x^2, I_y^2, I_{xy}\) 进行高斯加权(取 \(\sigma=1\) ),生成矩阵 \(M=\begin{bmatrix}A&C\\C&B\end{bmatrix}\) 的元素 \(A,B,C\) 。
\[A=g(I_x^2)=I_x^2 \star w, \; B=g(I_y^2)=I_y^2 \star w, \; C=g(I_{xy})=I_{xy} \star w\]
(4) 计算每个像素的Harris响应值 \(R\) ,并对小于某一阈值 \(t\) 的 \(R\) 置为零。
\[R=\{R: \det M - k(\text{trace} M)^2 < t\}\]
(5) 在3×3或5×5的邻域内进行“非最大值抑制”,局部最大值点即为图像中的角点。

3.1.6 OpenCV例子:Harris角点检测

OpenCV中实现了Harris角点检测,参见函数 cornerHarris ,下面是它的一个例子:

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import cv2
import numpy as np

filename = 'chessboard.jpg'
img = cv2.imread(filename)
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

gray = np.float32(gray)
dst = cv2.cornerHarris(gray,2,3,0.04)

#result is dilated for marking the corners, not important
dst = cv2.dilate(dst,None)

# Threshold for an optimal value, it may vary depending on the image.
img[dst>0.01*dst.max()]=[0,0,255]

cv2.imwrite('output.jpg', img)

这个例子中,假设输入图像chessboard.jpg如图 16 所示。

img_harris_example_chessboard.jpg

Figure 16: 待检测图像

则检测的输出图像output.jpg如图 17 所示。

img_harris_example_chessboard_output.jpg

Figure 17: Harris角点检测结果(红色标记)

参考:
http://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_feature2d/py_features_harris/py_features_harris.html#harris-corners

3.1.7 Harris角点检测算子的扩展

3.1.7.1 Harris-Laplace角点检测算子

Harris角点检测算子对亮度变化不敏感,且具有旋转不变性(回忆一下前面讨论,从各个方向移动小窗口 \(W\) ,检测自相关函数是否有较大变化)。但原始的Harris角点检测算子不具备尺度不变性。

将Harris角点检测算子和高斯尺度空间表示相结合,可以使Harris角点检测算子具有“尺度不变性”,这种方法也称为“Harris-Laplace检测方法”。

参考:https://en.wikipedia.org/wiki/Corner_detection#The_multi-scale_Harris_operator

3.1.7.2 Harris-Affine角点检测算子

Harris-Affine角点检测算子对原始Harris角点检测算子进行扩展,使其具备“仿射不变性”。

参考:K. Mikolajczyk, K. and C. Schmid (2004). "Scale and affine invariant interest point detectors"

3.2 Shi-Tomasi角点检测(Harris角点检测的变种)

Jianbo Shi和Carlo Tomasi在其论文 Good Features to Track 中提出的Shi-Tomasi算法,开始主要是为了解决跟踪问题,用来衡量两幅图像的相似度,我们可以把它看为Harris角点检测算法的变种。

Shi-Tomasi角点检测和Harris角点检测的不同之处仅在于角点响应的计算公式不同。
Harris角点检测中使用下面公式计算待检测点的响应:
\[R= \det M - k(\text{trace} M)^2 = \lambda_1 \lambda_2 - k (\lambda_1 + \lambda_2)^2\]
Shi-Tomasi角点检测中使用下面公式计算待检测点的响应:
\[R = \min(\lambda_1, \lambda_2)\]

Shi-Tomasi角点检测需要计算矩阵M的两个特征值。

3.3 SUSAN角点检测

SUSAN(Smallest Univalue Segment Assimilatating Nucleus的缩写),即“最小核值相似区”,是由牛津大学的Smith等提出来的。

这里介绍一下SUSAN角点检测的基本思想。

SUSAN检测方法是一种基于窗口模板的检测方法, 在图像每个像素点位置处建立一个圆形窗口,把圆中心像素与圆形模板内其他像素比较,统计出与圆中心像素相似的像素数量,当相似像素数量小于一个阈值,就被认为是要检测的角点。 如图 18 所示。

img_susan.png

Figure 18: SUSAN原理(在a圆形窗口中和它中心像素相似的像素数量最少)

SUSAN角点检测的一个优点就是对局部的噪声不敏感,抗干扰能力强,因为SUSAN算子不依赖图像分割的结果,避免了梯度计算。

参考:https://users.fmrib.ox.ac.uk/~steve/susan/susan/node2.html#SECTION00020000000000000000

4 Blob检测

斑点检测(Blob detection)是计算机视觉的重要内容,是区域检测(Region detection)的一种特例,是许多特征生成、目标识别等方法的重要预处理环节。

斑点通常指与周围有着颜色和灰度区别的区域。如,一棵树是一个斑点,一块草地、一栋房子也可看成一个斑点等等。斑点通常和关键点(key point),兴趣点(intrest point)以及特征点(feature point)表示同一个概念。 由于斑点代表的是一个区域,相比单纯的角点(Corner),它的稳定性要好,抗噪声能力要强。

4.1 利用LoG检测斑点

LoG是高斯拉普拉斯算子(Laplace of Guassian)的缩写,它的定义是:
\[\text{LoG} \triangleq \nabla^2 G(x,y,\sigma) = \frac{\partial^2}{\partial x^2} G(x,y,\sigma) + \frac{\partial^2}{\partial y^2} G(x,y,\sigma)\]

我们知道,高斯拉普拉斯算子 \(\nabla^2 G(x,y,\sigma)\) 和图像卷积后,再求“零交叉点”,即可得到图像的边缘(这就是Marr-Hildreth边缘检测的思路,可参考http://aandds.com/blog/dip.html)。

19 是利用高斯拉普拉斯算子检测一维信号中边缘的例子。

img_LoG_edge.jpg

Figure 19: LoG算子检测一维信号中边缘的例子

下面将介绍如何改造LoG算子,让其可以进行斑点检测。

4.1.1 一维信号斑点检测

一维信号中的斑点可以认为是两个相邻的突跳边缘组成。图 20 是不同大小的一维斑点。

img_1d_blob.jpg

Figure 20: 不同大小的一维斑点

21 是高斯函数( \(\sigma = 1\) )的二阶导数和图 20 中4个斑点信号卷积的结果。

img_1d_LoG_response.jpg

Figure 21: 斑点信号与高斯函数( \(\sigma = 1\) )二阶导数卷积的响应(检测右下子图的极值可以检测出信号斑点)

“零交叉点”对应原信号边缘。此外,我们发现,对于图 21 的右下子图, 检测LoG的“极值点”可以检测出原信号(图 20 右下子图)中的斑点。 检测出的斑点的尺寸(这里为2)正好是高斯函数方差( \(\sigma = 1\) )的2倍。

不过,使用检测LoG极值的方法只能检测出小斑点,而无法检测大斑点(如图 20 中上排两个子图、左下子图中的斑点就不法检测出来)。就算把 \(\sigma\) 增大也无法检测出斑点,如图 22 所示, 我们发现增大 \(\sigma\) 后,LoG响应会衰减,无法从LoG响应中检测出原信号斑点。

img_unnormalized_laplacian_response.jpg

Figure 22: LoG响应随 \(\sigma\) 增大而衰减

实验表明, 使用尺度归一化的LoG(Scale-normalized Laplacian of Gaussian)可以去除方差增大导致的衰减现象。尺度归一化LoG算子(二维情况)的定义为 \(\sigma^2 \nabla^2 G(x,y,\sigma)\) ,即:

\begin{align*} \nabla^2_{norm} G(x,y,\sigma) & \triangleq \sigma^2 \nabla^2 G(x,y,\sigma) \\ & = \sigma^2 \frac{\partial^2}{\partial x^2} G(x,y, \sigma) + \frac{\partial^2}{\partial y^2} G(x,y, \sigma) \\ & = - \frac{1}{\pi \sigma^2} \left(1 - \frac{x^2 + y^2}{2 \sigma^2} \right) e^{- \frac{x^2 + y^2}{2 \sigma^2}} \tag{尺度归一化LoG} \\ \end{align*}

23 是对一维信号的非规范化和尺度归一化LoG响应的对比图。

img_LoG_normLoG.jpg

Figure 23: 非规范化(子图a)和尺度归一化LoG(子图b)响应对比图

4.1.2 尺度归一化LoG检测斑点


前面以一维信号为例介绍了尺度归一化LoG算子检测斑点的过程。对于图像(二维信号)中斑点的检测过程也类似,基本过程为: 首先求尺度归一化LoG算子 \(\nabla^2_{norm} G(x,y,\sigma)\) 与图像的卷积,在响应图像中求极值点,极值点就对应于原图像中的斑点。

\[\nabla^2_{norm} G(x,y,\sigma) = - \frac{1}{\pi \sigma^2} \left(1 - \frac{x^2 + y^2}{2 \sigma^2} \right) e^{- \frac{x^2 + y^2}{2 \sigma^2}}\]

求 \(\nabla^2_{norm} G(x,y,\sigma)\) 的极值点,等价于求取下式:
\[\frac{\partial (\nabla^2_{norm} G(x,y,\sigma))}{\partial \sigma} =0\]
上式可化简为:
\[(x^2 + y^2 - 2 \sigma^2) e^{- \frac{x^2 + y^2}{2 \sigma^2}} = 0\]
即:
\[\sigma = (x^2 + y^2) / \sqrt{2}\]
例如,对于图 24 中的二值化圆形斑点,在尺度 \(\sigma = r / \sqrt{2}\) 时,尺度归一化LoG响应值达到最大。如果图 24 (a)中的圆形斑点黑白反向,则它的尺度归一化LoG响应值就在尺度 \(\sigma = r / \sqrt{2}\) 时达到最小。 尺度归一化LoG响应值到达峰值时的尺度 \(\sigma\) 值,称为特征尺度(Characteristic Scale)。

img_norm_LoG_example.jpg

Figure 24: 图像中的圆形斑点与其尺度归一化LoG响应

25 是实际检测斑点的过程示意图。对于二维图像,计算图像在不同尺度下的尺度归一化LoG响应值(即求 \(\nabla^2_{norm} G(x,y,\sigma)\) 与图像 \(I(x,y)\) 的卷积),然后,检查位置和尺度空间中每个响应值。 如果该点的响应值都大于或都小于其他26个立方空间邻域(图 25 中的26个圆圈)的值,那么该点就是被检测到的图像斑点。

img_maxima_and_minima.jpg

Figure 25: 在位置和尺度空间中搜索峰值点

4.2 利用DoH检测斑点

另一种斑点检测方法是利用图像点Hessian矩阵(二阶微分):
\[H(L) = \begin{bmatrix}L_{xx}&L_{xy}\\L_{xy}&L_{yy}\end{bmatrix}\]
上式中, \(L(x,y,\sigma) = G(x,y,\sigma) \star I(x,y), L_{xx} = \frac{\partial^2 L}{\partial x^2}, L_{yy} = \frac{\partial^2 L}{\partial y^2}, L_{xy} = \frac{\partial^2 L}{\partial x \partial y}\) 。

计算Hessian矩阵的行列式(Determinant of Hessian, DoH):
\[\text{DoH} = \det H(L) = L_{xx}L_{yy} - L^2_{xy}\]
尺度归一化的DoH为:
\[\det H_{norm}(L) \triangleq \sigma^4 ( L_{xx}L_{yy} - L^2_{xy} )\]

当 \(\det H_{norm}(L)\) 取极大值或极小值时的坐标点对应于图像中的斑点,这是使用DoH检测图像斑点的基本依据。 我们知道,使用LoG检测斑点时,斑点对应于 \(\sigma^2 (L_{xx} + L_{yy})\) 的极值点。而DoH则采用 \(\sigma^4 (L_{xx}L_{yy} - L^2_{xy})\) 的极值点来代表斑点。为什么DoH极值点可以代表斑点呢?我们简单地定性(不严谨)分析一下。忽略DoH中的 \(-L^2_{xy}\) ,由于“两数之和” \(L_{xx} + L_{yy}\) 取极值时,“两数之积” \(L_{xx}L_{yy}\) 往往也取极值,所以DoH极值点LoG极值点近似相同。

无论是LoG还是DoH,它们对图像中的斑点进行检测,其步骤都可以分为以下两步:
1)使用不同的 \(\sigma\) 生成 \(\frac{\partial^2 G}{\partial x^2} + \frac{\partial^2 G}{\partial y^2}\) (用于求LoG响应)或 \(\frac{\partial^2 G}{\partial x^2}, \frac{\partial^2 G}{\partial y^2}, \frac{\partial^2 G}{\partial xy}\)(用于求DoH响应)模板,并对图像进行卷积运算;
2)在图像的位置空间与尺度空间中搜索LoG与DoH响应的峰值。

参考:
Lindeberg, 2015 "Image matching using generalized scale-space interest points" (pdf)
https://en.wikipedia.org/wiki/Blob_detection#The_determinant_of_the_Hessian

4.2.1 LoG vs. DoH


与LoG相比, DoH响应仅当两个垂直方向上同时有大的变化时才会有较大的响应,所以,DoH对图像中的细长结构的斑点有较好的抑制作用。 如图 26 (摘自Lindeberg, 2015 "Image matching using generalized scale-space interest points")所示,通过对比第二排右子图(LoG检测结果)和第三排右子图(DoH检测结果)可知,细长斑点在DoH检测中被较好地抑制。

img_LoG_vs_DoH.jpg

Figure 26: 第一排(原图),第二排(LoG响应及其极值点),第三排(DoH响应及其极值点)

5 Blob检测高效算法:SIFT

尺度不变特征转换(Scale-Invariant Feature Transform)是一种机器视觉的算法,用于 检测与描述图像中的局部特征 ,它对旋转、尺度缩放、亮度变化保持不变性,对视角变化、仿射变换、噪声也保持一定程度的稳定性。此算法由 David Lowe 在1999年所发表,2004年完善总结。该算法不仅仅是个特征检测算法, 它将特征点检测和特征失量生成、特征匹配搜索等步骤完整地结合在一起进行优化,达到了接近实时的运算速度。 下面的描述中,不区分特征点、斑点、兴趣点等概念。

SIFT的四个步骤:
(1) Scale-space extrema detection(在尺度空间上查找特征点);
(2) Keypoint localization(特征点精确定位);
(3) Orientation assignment(为每个特征点指定方向参数);
(4) Keypoint descriptor(生成特征点描述子)

参考:
David G. Lowe. (2004), "Distinctive Image Features from Scale-Invariant Keypoints" (pdf)

5.1 SIFT算法检测特征点

SIFT算法的第一步是在尺度空间(Scale-space)上查找特征点。

SIFT算法中,尺度空间在实现时使用高斯金字塔表示,高斯金字塔的构建分为两部分:
(1) 对图像做不同尺度的高斯模糊,这样得到了一组(Octave)尺度空间;
(2) 对图像做降采样(隔点采样),这样得到了多组尺度空间。

5.1.1 初步求原始极值点

4.1.2 中介绍了使用“尺度归一化LoG算子” \(\nabla^2_{norm} G(x,y,\sigma) = \sigma^2 \nabla^2 G(x,y,\sigma)\) 检测斑点的方法。SIFT算法没有直接采用尺度归一化LoG算子 \(\sigma^2 \nabla^2 G(x,y,\sigma)\) (它的计算效率不高),而是利用高斯函数的差分(Difference of Gaussian, DoG)来对 \(\sigma^2 \nabla^2 G(x,y,\sigma)\) 进行近似(这里不推导它们为什么近似相同)。DoG的定义为:
\[\begin{aligned} D(x, y, \sigma) &= (G(x,y,k \sigma) - G(x,y,\sigma)) \star I(x,y) \\ &= L(x,y,k \sigma) - L(x,y,\sigma) \end{aligned}\]
其中, \(k \sigma\) 和 \(\sigma\) 是连续的两个图像的高斯平滑尺度。DoG非常容易计算,两个高斯图像相差即可,如图 27 所示。

img_SIFT_DoG.jpg

Figure 27: DoG的计算

在DoG中检测特征点的过程和节 4.1.2 中描述的在尺度归一化LoG响应中检测斑点过程是相同(即如图 25 所示)。

得到原始极值点的过程可总结为图 28 (摘自http://blog.csdn.net/lhanchao/article/details/52345845 )所示。

img_SIFT_initial_keypoints.png

Figure 28: SIFT得到初步极值点过程

注:SIFT论文中建议采用4个octaves,第个octave中采用5个高斯模糊层。

5.1.2 精确定位极值点

在节 2.3 中介绍过,在尺度空间中寻找到的特征点的位置信息往往不够精确。图 29 是一维情况下采样后极值点偏离真正极值点的例子。

img_SIFT_sample.jpg

Figure 29: 采样后极值点偏离真正极值点

一般,采用“子像元插值(sub pixel interpolation)”的方式来求得更准确的极值点。

首先,我们看一个一维函数插值的例子,如图 30 所示。

img_SIFT_1d_func.jpg

Figure 30: 已知离散一维函数值,求取连续一维函数值极值点的例子

在 \(x_0 = 0\) (它是几个离散点中的极值点)处进行 泰勒展开 ,有:
\[\begin{aligned} f(x_0 + h) & \approx f(x_0) + f'(x_0) h + \frac 1 2 f''(x_0) h^2 \\ & = 6 + 2 h + \frac{-6}{2} h^2 \\ & = 6 + 2 h - 3 h^2 \end{aligned}\]
求 \(f(x_0 + h)\) (看作是 \(h\) 的函数)极值,令 \(f'(x_0 + h) = 2 - 6h = 0\) 得 \(h = \frac{1}{3}\) 。从而,极值点的极值为:
\[f(x_0 + h) = 6 + 2 \times \frac{1}{3} - 3 \times (\frac{1}{3})^2 = 6\frac{1}{3}\]
对于3维函数 \(D(\boldsymbol{X}) = D(x, y, \sigma)\) 在 \(\boldsymbol{X}_0\) 处进行泰勒展开(向量形式)有:
\[D(\boldsymbol{X}_0 + \boldsymbol{h}) \approx D(\boldsymbol{X}_0) + \left. \left(\frac{\partial D}{\partial \boldsymbol{X}} \right)^{\mathsf{T}}\right|_{\boldsymbol{X} = \boldsymbol{X}_0} \boldsymbol{h} + \frac{1}{2} \boldsymbol{h}^{\mathsf{T}} \left. \left(\frac{\partial^2 D}{\partial \boldsymbol{X}^2} \right) \right|_{\boldsymbol{X} = \boldsymbol{X}_0} \boldsymbol{h}\]
其中:
\[\frac {\partial D}{\partial \boldsymbol{X}} = \begin{bmatrix} \frac {\partial D}{\partial x} \\ \frac {\partial D}{\partial y} \\ \frac {\partial D}{\partial \sigma}\end{bmatrix} = \begin{bmatrix} \frac {D(x+1,y,\sigma) - D(x-1,y,\sigma)}{2} \\ \frac {D(x,y+1,\sigma)-D(x,y-1,\sigma)}{2} \\ \frac {D(x,y,\sigma+1) - D(x,y,\sigma-1)}{2} \end{bmatrix}, \quad \frac{\partial^2 D}{\partial \boldsymbol{X}^2} = \begin{bmatrix} D_{xx} & D_{xy} & D_{x\sigma} \\ D_{yx} & D_{yy} & D_{y\sigma} \\ D_{\sigma x} & D_{\sigma y} & D_{\sigma\sigma} \end{bmatrix}\]
\(\frac{\partial^2 D}{\partial \boldsymbol{X}^2}\) 中的元素 \(D_{xy}\) 为(其它元素类似):
\[D_{xy} = \frac { \frac{D(x+1,y+1,s)-D(x-1,y+1,s)}{2} - \frac{D(x+1,y-1,s)-D(x-1,y-1,s)}{2} }{2}\]
求 \(D(\boldsymbol{X}_0 + \boldsymbol{h})\) (看作是 \(\boldsymbol{h}\) 的函数)的极值,令其导数为零,可得极值点:
\[\boldsymbol{h} = - \left(\frac{\partial^2 D}{\partial \boldsymbol{X}^2} \right)^{-1} \frac{\partial D}{\partial \boldsymbol{X}}\]

精确定位极值点这个步骤主要进行下面两个处理:
处理一: 如果 \(\boldsymbol{h}\) (它是相对于插值中心点的偏移量)在任何方向上的偏移大于0.5(两个像素的一半距离)时,说明插值中心点已经偏移到它的邻近点上,这样的点要删除。
处理二: 计算极值点 \(\boldsymbol{X}_0 + \boldsymbol{h}\) 处的响应 \(D(\boldsymbol{X}_0 + \boldsymbol{h})\) ,如果 \(|D(\boldsymbol{X}_0 + \boldsymbol{h})| < 0.03\) (假设图像的灰度值是在0~1.0之间),则认为响应值过小,这样的点容易受噪声干扰而不稳定,所以也要把它删除。

极值点处的极值 \(D(\boldsymbol{X}_0 + \boldsymbol{h})\) 的计算公式为 \(D(\boldsymbol{X}_0 + \boldsymbol{h}) = D(\boldsymbol{X}_0) + \frac{1}{2} \left. \left(\frac{\partial D}{\partial \boldsymbol{X}} \right)^{\mathsf{T}}\right|_{\boldsymbol{X} = \boldsymbol{X}_0} \boldsymbol{h}\) ,其推导过程如下:
把 \(\boldsymbol{h}\) 代入到 \(D(\boldsymbol{X}_0 + \boldsymbol{h})\) 的泰勒展开式中的最后一部分中:
\[\begin{aligned} D(\boldsymbol{X}_0 + \boldsymbol{h}) & = D(\boldsymbol{X}_0) + \left(\frac{\partial D}{\partial \boldsymbol{X}} \right)^{\mathsf{T}} \boldsymbol{h} + \frac{1}{2} \boldsymbol{h}^{\mathsf{T}} \left(\frac{\partial^2 D}{\partial \boldsymbol{X}^2} \right) \boldsymbol{h} \\ & = D(\boldsymbol{X}_0) + \left(\frac{\partial D}{\partial \boldsymbol{X}} \right)^{\mathsf{T}} \boldsymbol{h} + \frac{1}{2} \left( - \left(\frac{\partial^2 D}{\partial \boldsymbol{X}^2} \right)^{-1} \frac{\partial D}{\partial \boldsymbol{X}} \right)^{\mathsf{T}} \left(\frac{\partial^2 D}{\partial \boldsymbol{X}^2} \right) \left( - \left(\frac{\partial^2 D}{\partial \boldsymbol{X}^2} \right)^{-1} \frac{\partial D}{\partial \boldsymbol{X}} \right) \\ & = D(\boldsymbol{X}_0) + \left(\frac{\partial D}{\partial \boldsymbol{X}} \right)^{\mathsf{T}} \boldsymbol{h} + \frac{1}{2} \left(\frac{\partial D}{\partial \boldsymbol{X}} \right)^{\mathsf{T}} \left(\frac{\partial^2 D}{\partial \boldsymbol{X}^2} \right)^{- \mathsf{T}} \frac{\partial^2 D}{\partial \boldsymbol{X}^2} \left(\frac{\partial^2 D}{\partial \boldsymbol{X}^2} \right)^{-1} \frac{\partial D}{\partial \boldsymbol{X}} \\ & = D(\boldsymbol{X}_0) + \left(\frac{\partial D}{\partial \boldsymbol{X}} \right)^{\mathsf{T}} \boldsymbol{h} + \frac{1}{2} \left(\frac{\partial D}{\partial \boldsymbol{X}} \right)^{\mathsf{T}} \left(\frac{\partial^2 D}{\partial \boldsymbol{X}^2} \right)^{-1} \frac{\partial D}{\partial \boldsymbol{X}} \\ & = D(\boldsymbol{X}_0) + \left(\frac{\partial D}{\partial \boldsymbol{X}} \right)^{\mathsf{T}} \boldsymbol{h} + \frac{1}{2} \left(\frac{\partial D}{\partial \boldsymbol{X}} \right)^{\mathsf{T}} (- \boldsymbol{h}) \\ & = D(\boldsymbol{X}_0) + \frac{1}{2} \left(\frac{\partial D}{\partial \boldsymbol{X}} \right)^{\mathsf{T}} \boldsymbol{h} \\ \end{aligned}\]

5.1.3 去除不稳定极值点(删除边缘效应)

4.2.1 中介绍过通过LoG求斑点时,细节结构的斑点(边缘)不会被较好地抑制。边缘是不稳定的,很难定位,应该把它们去掉。SIFT采用DoG来近似LoG,并没有克服这个缺点,所以我们需要想办法来删除位于边缘部分的极值点。

我们通过极值点邻域的曲率来判断的当前极值点是否在边缘上。在跨越边缘的方向有较大的主曲率,在与边缘相切的方向主曲率较小。

主曲率可以通过下面2阶Hessian方阵求出:
\[H(x,y) = \begin{bmatrix} D_{xx} & D_{xy} \\ D_{xy} & D_{yy} \end{bmatrix}\]
某点的主曲率和该点的Hessian矩阵的特征值是成比例的,因此可以通过Hessian矩阵的特征值来确定某点在差分高斯金字塔中的主曲率。

我们可以只关心两个特征值的比例,当比例较小时,说明两个特征值相差不大,从而不同方向的主曲率相关不大。而当两个特征值比例很大时,则说明不同方向的主曲率相关比较大,这样的极值点应该删除。

由于只关心特征值的比例,采用和Harris角点检测类似的思路,我们可以避免直接求解Hessian矩阵特征值。具体过程如下所述。设Hessian矩阵的两个特征值分别为 \(\alpha, \beta\) ,则:
\[\begin{aligned} \text{trace} H &= D_{xx} + D_{yy} = \alpha + \beta \\ \det H &= D_{xx}D_{yy} - (D_{xy})^2 = \alpha \beta \end{aligned}\]
记 \(\gamma\) 为较大特征值为较小特征值的比值,不失一般性假设 \(\alpha\) 比较大,则有 \(\alpha = \gamma \beta\) ,从而:
\[\frac{\text{trace}(H)^2}{\det(H)} = \frac{(\alpha + \beta)^2}{\alpha \beta} = \frac{(\gamma \beta + \beta)^2 }{\gamma \beta^2} = \frac{(\gamma + 1)^2}{\gamma}\]
上式左边可通过Hessian矩阵求解得到,它仅与两个特征值的比例有关,而与具体特征值无关。当两个特征值相等时, \(\frac{(\gamma + 1)^2}{\gamma}\) 的值为最小,随着 \(\gamma\) 增大, \(\frac{(\gamma + 1)^2}{\gamma}\) 的值也增加。所以,要想检测主曲率比较小于某一阈值 \(\gamma\) ,只用检测下式是否成立:
\[\frac{\text{trace}(H)^2}{\det(H)} < \frac{(\gamma + 1)^2}{\gamma}\]
在SIFT论文中,给出的曲率比率阈值 \(\gamma=10\) ,也就是说,只保留主曲率比例小于10的特征点,对于主曲率比例大于10的特征点将被删除。

5.1.4 检测特征点过程演示

31 (摘自 SIFT wikipedia )是检测SIFT算法检测特征点的过程演示。

img_SIFT_keypoints_filtering.jpg

Figure 31: 上图为通过DoG求得的原始特征点;中图为通过子像元插值,去掉低响应点后的特征点;下图为删除边缘效应后的特征点

5.2 SIFT特征描述子

在特征点提取之后,如何使用或选择一种合适的特征描述方法来描述特征点附近局部图像模式,便是图像识别过程中又一重要步骤。图像识别过程往往是两幅图像之间进行相似性比较的过程,为了实现相似性度量,使用或选择一种紧凑而完整的特征描述是十分重要的。特征描述子(Descriptor)是一种图像局部结构特征的定量化数据描述,它能够充分反映特征点附近局部图像的形状和纹理结构特性。

5.2.1 特征点方向分配

为了让特征点具有图像旋转不变性,需要根据检测到的特征点的局部图像结构求得一个方向基准。 我们使用图像梯度的方法求取局部结构的稳定方向。对于已经检测到的特征点,我们知道了该特征点的尺度值 \(\sigma\) ,这个尺度对应的高斯图像为:
\[L(x,y) = G(x,y, \sigma) \star I(x,y)\]
计算以特征点为中心,以 \(3 \times 1.5 \sigma\) 为半径的区域内的所有图像像素点的梯度幅角和幅值(模值)。其计算公式分别为:
\[\begin{aligned} 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) &= \arctan \left(\frac{L(x, y+1) - L(x, y-1)}{L(x+1,y) - L(x-1,y)}) \right) \end{aligned}\]

完成特征点邻域的高斯图像的梯度计算后,使用直方图统计邻域内像素的梯度方向和幅值。梯度方向直方图的横轴是梯度方向角,纵轴是梯度方向角对应的梯度幅值累加值。梯度方向直方图将0-360度的范围分为36个柱(bins),每10度为一个柱。 直方图的峰值代表该特征点处邻域内图像梯度的主方向,也即该特征点的主方向。 如图 32 (摘自:http://aishack.in/tutorials/sift-scale-invariant-feature-transform-keypoint-orientation/ )所示。不过,每个累加到梯度方向直方图的采样点的梯度幅值都要进行权重处理,加权采用圆形高斯加权函数,高斯加权函数的 \(\sigma\) 值为特征点尺度的1.5倍。 通过高斯加权,使特征点附近的梯度幅值有较大的权重,这样可以部分弥补SIFT算法没有考虑仿射不变性而产生的特征点不稳定的问题。

img_SIFT_orientation_histogram.jpg

Figure 32: 梯度方向直方图(直方图的峰值代表该特征点的主方向)

如果在梯度方向直方图中存在一个相当于主峰值能量80%能量的峰值时(如图 32 所示),则将这个方向认为是该特征点的辅方向。一个特征点可能会被指定具有多个方向(一个主方向,一个以上辅方向),这可以增强匹配的鲁棒性,如下图所示。实际编程实现中,就是把该特征点复制成多份特征点,并将方向值分别赋给这些复制后的特征点。通常,离散的梯度方向直方图要进行插值拟合处理,这样可以求得更精确的方向角度值。

在获得了图像特征点主方向 \(\theta\) 后,每个特征点有三个信息:位置、尺度、方向,它们可表示为 \((x,y,\sigma, \theta)\) 。

5.2.2 特征点特征矢量生成

SIFT描述子是对特征点附近“邻域内高斯图像梯度”统计结果的一种表示,它是一个三维的阵列,但通常在编程将它表示成一个矢量(1维向量),这个矢量是通过对三维阵列按一定规律进行排列得到的。

特征描述子与特征点所在的尺度有关,因此,对梯度的求取应在特征点对应的高斯图像上进行。将特征点附近的邻域划分成 \(4 \times 4\) 个子区域,每个子区域的尺寸为 \(3 \sigma\) 个像素。

为了保证特征矢量具有旋转不变形,需要以特征点为中心,将特征点附近邻域内图像的梯度的位置和方向旋转一个方向角 \(\theta\) (即特征点的主方向),如图 33 所示。

img_SIFT_rotate_theta.jpg

Figure 33: 将特征点附近邻域内图像梯度的位置和方向旋转方向角 \(\theta\) (即特征点的主方向)

特征点邻域划分成 \(4 \times 4\) 个子区域,在每个子区域内计算8个方向(每个方向45度,注:在计算特征点主方向时,每个方向是10度,共有36个方向)的“梯度方向直方图”(注:在计算特征点主方向时也用到了梯度方向直方图),绘制每个梯度方向的累加值,形成一个种子点。这样, 每个种子点共有8个方向的梯度强度信息。由于存在 \(4 \times 4\) 个子区域,所以,共有 \(4 \times 4 \times 8 = 128\) 个数据,最终形成128维的SIFT特征矢量。 如图 34 所示,128维的每一维都是某个子区域(共16个子区域)所有像素的某个方向(共8个方向)的梯度累加值。

img_SIFT_128.jpg

Figure 34: SIFT特征矢量由128维矢量组成

特征向量形成后,记为 \(W=(w_1, w_2, \cdots, w_{128})\) ,为了去除光照变化(灰度整体漂移)的影响,需要对它们进行归一化处理,归一化后的向量为 \(L=(l_1, l_2, \cdots, l_{128})\) ,计算公式为:
\[l_i = \frac{w_i}{\sum_{j=1}^{128}w_j}, \quad i = 1,2, \cdots, 128\]

非线性光照,相机饱和度变化可能造成某些方向的梯度值过大,而对方向的影响微弱。因此设置门限值(向量归一化后,一般取门限值为0.2)截断较大的梯度值,即大于0.2的值只取为0.2。然后,再进行一次归一化处理。


Author: cig01

Created: <2016-07-15 Fri 00:00>

Last updated: <2018-02-10 Sat 20:33>

Creator: Emacs 25.3.1 (Org mode 9.1.4)