Local Invariant Feature Detectors

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

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

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

图像处理中，仿射变换还可以写为：
\[\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）是仿射变换的一些例子。

由于仿射变换的 \(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

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

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

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

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

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

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

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

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

根据二元函数泰勒展开，对图像 \(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}\]
其中：

\(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 所示。

前面介绍过，判断点 \((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)\) 是角点。

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

矩阵 \(M(x,y)\) 的特征值与图像中不同类型（直线/平面/角点）的区域有下面关系（如图 15 所示）：
(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\) 都有较大的变化）。

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

上式中， \(\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\) 都很大，且近似相等，相应图像区域为角点。

根据上述讨论，可以将 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 的邻域内进行“非最大值抑制”，局部最大值点即为图像中的角点。

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 所示。

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

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

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

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

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

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

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

图 23 是对一维信号的非规范化和尺度归一化 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）。

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

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

节 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 所示。

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

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

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

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

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

首先，我们看一个一维函数插值的例子，如图 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}\]

节 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 的特征点将被删除。

图 31 （摘自 SIFT wikipedia ）是检测 SIFT 算法检测特征点的过程演示。

为了让特征点具有图像旋转不变性，需要根据检测到的特征点的局部图像结构求得一个方向基准。 我们使用图像梯度的方法求取局部结构的稳定方向。对于已经检测到的特征点，我们知道了该特征点的尺度值 \(\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 算法没有考虑仿射不变性而产生的特征点不稳定的问题。

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

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

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

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

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

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

特征向量形成后，记为 \(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。然后，再进行一次归一化处理。