SIFT & SURF
两者都是特征检测的算法，但略有不同。
一、SURF
1. 背景及来源
使用SIFT算法进行关键点检测和描述的执行速度比较慢，需要速度更快的算法。SURF算法由Bay在2006年提出，是SIFT算法的增强版，其计算量小，运算速度快，提取的特征和SIFT几乎无二。
2. 二者对比
特征点检测
SIFT：使用不同尺度的图片与高斯函数进行卷积
SURF：使用不同大小的盒状滤波器和原始图像卷积，易于并行
方向
SIFT：关键点邻接矩形区域内，利用梯度直方图计算
SURF：关键点邻接圆域内，计算x，y方向的haar小波
描述符生成
SIFT：关键点邻域内划分d*d个子区域，每个子区域内计算8个方向的直方图
SURF：关键点邻域内划分d*d个子区域，每个子区域计算采样点的haar小波响应，并记录：∑ d x , ∑ d y , ∑ ∣ d x ∣ , ∑ ∣ d y ∣ \sum dx,\sum dy,\sum |dx|, \sum |dy| ∑ d x , ∑ d y , ∑ ∣ d x ∣ , ∑ ∣ d y ∣
二、SIFT(Scale-invariant feature transform)
尺度不变特征转换，之前提过角点检测，但是当放大图像时，原窗口内的像素集合可能构不成角点，所以我们就提出了SIFT解决此类问题。
1. 概述
来源
作者：David Lowe
专利权属于英属哥伦比亚大学
算法示例
2. 建立高斯差分金字塔
示例
高斯金字塔
高斯差分金字塔
其中，Ocative表示第 n 组，每组的层数不固定，但是上下层之间的像素大小是相差2倍(上面比下面小两倍)。以第一组为例，我们先使用方差为σ \sigma σ 2 σ 2\sigma 2 σ
层级关系
一幅图像的最佳组数：O = l o g 2 ( m i n ( M , N ) ) − 3 O = log_2(min(M,N)) - 3 O = l o g 2 ( min ( M , N )) − 3
每组的最佳层数：S = n + 3 S = n + 3 S = n + 3
各组层的σ \sigma σ
k = 2 1 n , σ 0 = 1. 6 2 − 0. 5 2 = 1.52 k = 2^{\frac{1}{n}}, \sigma_0 = \sqrt{1.6^2 - 0.5^2} = 1.52 k = 2 n 1 , σ 0 = 1. 6 2 − 0. 5 2 = 1.52
其中，将下面一组的倒数第二层作为上面一组的倒数第一层，因为下采样的干系。
3. 极值点的精确定位
阈值化
a b s ( v a l ) = 0.5 ∗ T / n abs(val) = 0.5*T/n ab s ( v a l ) = 0.5 ∗ T / n T是常数，论文中为 0.04。
在高斯差分金字塔中找极值
图示
调整极值点位置
在第二步检测到的极值点X 0 ( x 0 , y 0 , σ 0 ) T X_0(x_0,y_0,\sigma_0)^T X 0 ( x 0 , y 0 , σ 0 ) T f ( [ x y σ ] ) = f ( [ x 0 y 0 σ 0 ] ) + [ ∂ f ∂ x , ∂ f ∂ y , ∂ f ∂ σ ] ( [ x y σ ] − [ x 0 y 0 σ 0 ] ) + 1 2 ( [ x y σ ] − [ x 0 y 0 σ 0 ] ) T [ ∂ 2 f ∂ x ∂ x , ∂ 2 f ∂ x ∂ y , ∂ 2 f ∂ x ∂ σ ∂ 2 f ∂ x ∂ y , ∂ 2 f ∂ y ∂ y , ∂ 2 f ∂ y ∂ σ ∂ 2 f ∂ x ∂ σ , ∂ 2 f ∂ y ∂ σ , ∂ 2 f ∂ σ ∂ σ ] ( [ x y σ ] − [ x 0 y 0 σ 0 ] ) f\left(\left[\begin{array}{l}x \\y \\\sigma\end{array}\right]\right)=f\left(\left[\begin{array}{l}x_{0} \\y_{0} \\\sigma_{0}\end{array}\right]\right)+\left[\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}, \frac{\partial f}{\partial \sigma}\right]\left(\left[\begin{array}{l}x \\y \\\sigma\end{array}\right]-\left[\begin{array}{l}x_{0} \\y_{0} \\\sigma_{0}\end{array}\right]\right)+\frac{1}{2}\left(\left[\begin{array}{l}
x \\
y \\
\sigma
\end{array}\right]-\left[\begin{array}{l}
x_{0} \\
y_{0} \\
\sigma_{0}
\end{array}\right]\right)^{T}\left[\begin{array}{l}
\frac{\partial^{2} f}{\partial x \partial x}, \frac{\partial^{2} f}{\partial x \partial y}, \frac{\partial^{2} f}{\partial x \partial \sigma} \\
\frac{\partial^{2} f}{\partial x \partial y}, \frac{\partial^{2} f}{\partial y \partial y}, \frac{\partial^{2} f}{\partial y \partial \sigma} \\
\frac{\partial^{2} f}{\partial x \partial \sigma}, \frac{\partial^{2} f}{\partial y \partial \sigma}, \frac{\partial^{2} f}{\partial \sigma \partial \sigma}
\end{array}\right]\left(\left[\begin{array}{l}
x \\y \\\sigma\end{array}\right]-\left[\begin{array}{l}x_{0} \\y_{0} \\\sigma_{0}\end{array}\right]\right) f ⎝ ⎛ ⎣ ⎡ x y σ ⎦ ⎤ ⎠ ⎞ = f ⎝ ⎛ ⎣ ⎡ x 0 y 0 σ 0 ⎦ ⎤ ⎠ ⎞ + [ ∂ x ∂ f , ∂ y ∂ f , ∂ σ ∂ f ] ⎝ ⎛ ⎣ ⎡ x y σ ⎦ ⎤ − ⎣ ⎡ x 0 y 0 σ 0 ⎦ ⎤ ⎠ ⎞ + 2 1 ⎝ ⎛ ⎣ ⎡ x y σ ⎦ ⎤ − ⎣ ⎡ x 0 y 0 σ 0 ⎦ ⎤ ⎠ ⎞ T ⎣ ⎡ ∂ x ∂ x ∂ 2 f , ∂ x ∂ y ∂ 2 f , ∂ x ∂ σ ∂ 2 f ∂ x ∂ y ∂ 2 f , ∂ y ∂ y ∂ 2 f , ∂ y ∂ σ ∂ 2 f ∂ x ∂ σ ∂ 2 f , ∂ y ∂ σ ∂ 2 f , ∂ σ ∂ σ ∂ 2 f ⎦ ⎤ ⎝ ⎛ ⎣ ⎡ x y σ ⎦ ⎤ − ⎣ ⎡ x 0 y 0 σ 0 ⎦ ⎤ ⎠ ⎞
矢量形式:f ( X ) = f ( X 0 ) + ∂ f T ∂ X X ^ + 1 2 X ^ τ ∂ 2 f ∂ X 2 X ^ f(X)=f\left(\boldsymbol{X}_{0}\right)+\frac{\partial f^{T}}{\partial \boldsymbol{X}} \widehat{X}+\frac{1}{2} \widehat{X}^{\tau} \frac{\partial^{2} f}{\partial X^{2}} \widehat{X} f ( X ) = f ( X 0 ) + ∂ X ∂ f T X + 2 1 X τ ∂ X 2 ∂ 2 f X
具体过程：
f ( x ) f(x) f ( x ) ∂ f ( X ) ∂ X = ∂ f T ∂ X + 1 2 ( ∂ 2 f ∂ X 2 + ∂ 2 f T ∂ X 2 ) ℓ = ∂ f T ∂ X + ∂ 2 f ∂ X 2 X ^ \frac{\partial f(X)}{\partial \boldsymbol{X}}=\frac{\partial f^{T}}{\partial \boldsymbol{X}}+\frac{1}{2}\left(\frac{\partial^{2} f}{\partial \boldsymbol{X}^{2}}+\frac{\partial^{2} f^{T}}{\partial \boldsymbol{X}^{2}}\right) \ell=\frac{\partial f^{T}}{\partial \boldsymbol{X}}+\frac{\partial^{2} f}{\partial \boldsymbol{X}^{2}} \hat{X} ∂ X ∂ f ( X ) = ∂ X ∂ f T + 2 1 ( ∂ X 2 ∂ 2 f + ∂ X 2 ∂ 2 f T ) ℓ = ∂ X ∂ f T + ∂ X 2 ∂ 2 f X ^ X ^ = x − X 0 \hat{X} = x - X_0 X ^ = x − X 0 令导数为0解得:X ^ = − ∂ 2 f − 1 ∂ X 2 ∂ f ∂ X \hat{X}=-\frac{\partial^{2} f^{-1}}{\partial X^{2}} \frac{\partial f}{\partial X} X ^ = − ∂ X 2 ∂ 2 f − 1 ∂ X ∂ f
代入f ( x ) f(x) f ( x ) f ( X ) : f ( X ) = f ( X 0 ) + ∂ f T ∂ X X ^ + 1 2 ( − ∂ 2 f − 1 ∂ X 2 ∂ f ∂ X ) T ∂ 2 f ∂ X 2 ( − ∂ 2 f − 1 ∂ X 2 ∂ f ∂ X ) = f ( X 0 ) + ∂ f T ∂ X X ^ + 1 2 ∂ f T ∂ X ∂ 2 f − T ∂ X 2 ∂ 2 f ∂ X 2 ∂ 2 f − 1 ∂ X 2 ∂ f ∂ X = f ( X 0 ) + ∂ f T ∂ X X ^ + 1 2 ∂ f T ∂ X ∂ 2 f − 1 ∂ X 2 ∂ f ∂ X = f ( X 0 ) + ∂ f T ∂ X X ^ + 1 2 ∂ f T ∂ X ( − X ^ ) = f ( X 0 ) + 1 2 ∂ f T ∂ X X ^ \begin{array}{l} \boldsymbol{f}(\mathrm{X}): f(X)=f\left(X_{0}\right)+\frac{\partial f^{T}}{\partial X} \widehat{X}+\frac{1}{2}\left(-\frac{\partial^{2} f^{-1}}{\partial X^{2}} \frac{\partial f}{\partial X}\right)^{T} \frac{\partial^{2} f}{\partial X^{2}}\left(-\frac{\partial^{2} f^{-1}}{\partial X^{2}} \frac{\partial f}{\partial X}\right)
\\=f\left(X_{0}\right)+\frac{\partial f^{T}}{\partial X} \widehat{X}+\frac{1}{2} \frac{\partial f^{T}}{\partial X} \frac{\partial^{2} f^{-T}}{\partial X^{2}} \frac{\partial^{2} f}{\partial X^{2}} \frac{\partial^{2} f^{-1}}{\partial X^{2}} \frac{\partial f}{\partial X} \\
=f\left(X_{0}\right)+\frac{\partial f^{T}}{\partial X} \hat{X}+\frac{1}{2} \frac{\partial f^{T}}{\partial X} \frac{\partial^{2} f^{-1}}{\partial X^{2}} \frac{\partial f}{\partial X} \\
=f\left(X_{0}\right)+\frac{\partial f^{T}}{\partial X} \hat{X}+\frac{1}{2} \frac{\partial f^{T}}{\partial X}(-\widehat{X}) \\
=f\left(X_{0}\right)+\frac{1}{2} \frac{\partial f^{T}}{\partial X} \hat{X}
\end{array} f ( X ) : f ( X ) = f ( X 0 ) + ∂ X ∂ f T X + 2 1 ( − ∂ X 2 ∂ 2 f − 1 ∂ X ∂ f ) T ∂ X 2 ∂ 2 f ( − ∂ X 2 ∂ 2 f − 1 ∂ X ∂ f ) = f ( X 0 ) + ∂ X ∂ f T X + 2 1 ∂ X ∂ f T ∂ X 2 ∂ 2 f − T ∂ X 2 ∂ 2 f ∂ X 2 ∂ 2 f − 1 ∂ X ∂ f = f ( X 0 ) + ∂ X ∂ f T X ^ + 2 1 ∂ X ∂ f T ∂ X 2 ∂ 2 f − 1 ∂ X ∂ f = f ( X 0 ) + ∂ X ∂ f T X ^ + 2 1 ∂ X ∂ f T ( − X ) = f ( X 0 ) + 2 1 ∂ X ∂ f T X ^
当然，这也仅仅只是粗略计算，还存在一些细节问题：迭代次数限制、解超出一定范围舍去。
舍去低对比度的点
若∣ f ( X ) ∣ < T n |f(X)|<\frac{T}{n} ∣ f ( X ) ∣ < n T
边缘效应的去除
海森矩阵：H ( x , y ) = [ D x x ( x , y ) D x y ( x , y ) D x y ( x , y ) D y y ( x , y ) ]
\boldsymbol{H}(x, y)=\left[\begin{array}{ll}
D_{x x}(x, y) & D_{x y}(x, y) \\
D_{x y}(x, y) & D_{y y}(x, y)
\end{array}\right] H ( x , y ) = [ D xx ( x , y ) D x y ( x , y ) D x y ( x , y ) D yy ( x , y ) ]
Tr ( H ) = D x x + D y y = α + β Det ( H ) = D x x D y y − ( D x y ) 2 = α β \operatorname{Tr}(\boldsymbol{H})=D_{x x}+D_{y y}=\alpha+\beta\\\operatorname{Det}(\boldsymbol{H})=D_{x x} D_{y y}-\left(D_{x y}\right)^{2}=\alpha \beta Tr ( H ) = D xx + D yy = α + β Det ( H ) = D xx D yy − ( D x y ) 2 = α β α > β 并且 α = γ β \alpha>\beta 并且 \alpha=\gamma \beta α > β 并且 α = γ β γ \gamma γ
若D e t ( H ) < 0 Det(H)<0 De t ( H ) < 0 Tr ( H ) 2 Det ( H ) = ( α + β ) 2 α β = ( γ β + β ) 2 γ β 2 = ( γ + 1 ) 2 γ \frac{\operatorname{Tr}(\boldsymbol{H})^{2}}{\operatorname{Det}(\boldsymbol{H})}=\frac{(\alpha+\beta)^{2}}{\alpha \beta}=\frac{(\gamma \beta+\beta)^{2}}{\gamma \beta^{2}}=\frac{(\gamma+1)^{2}}{\gamma} Det ( H ) Tr ( H ) 2 = α β ( α + β ) 2 = γ β 2 ( γ β + β ) 2 = γ ( γ + 1 ) 2 Tr ( H ) 2 Det ( H ) < ( γ + 1 ) 2 γ \frac{\operatorname{Tr}(\boldsymbol{H})^{2}}{\operatorname{Det}(\boldsymbol{H})}<\frac{(\gamma+1)^{2}}{\gamma} Det ( H ) Tr ( H ) 2 < γ ( γ + 1 ) 2
4. 为关键点赋予方向
思想：统计以特征点为圆心，以该特征点所在的高斯图像的尺度为1.5倍为半径的圆内的所有的像素的梯度方向及其梯度幅值，并作1.5 σ 1.5\sigma 1.5 σ
示意图(在最接近关键点尺度值σ \sigma σ
5. 构建关键点的描述符
关键点的匹配需要用到描述符(KNN)。
旋转不变性
在SIFT中，我们想要让图像在旋转一定角度后仍保持原来特征，所以我们需要将整个特征点描述符的区域正向改为梯度幅度最大方向。
图示：
区域大小
其中，m是3，d是每个小区域的边长是多少(像素单位)
具体过程：
构建一个128维的特征向量
把关键点周围的邻域分为16个子区域(长宽都为4)
在每个子区域内统计8个方向的梯度值
补充：有限差分求导法(利用差值求出导数)
单层差分求导:
图示：
( ∂ f ∂ x ) = f 1 − f 3 2 h ( ∂ f ∂ y ) = f 2 − f 4 2 h ( ∂ 2 f ∂ x 2 ) = f 1 + f 3 − 2 f 0 h 2 ( ∂ 2 f ∂ y 2 ) = f 2 + f 4 − 2 f 0 h 2 ( ∂ 2 f ∂ x ∂ y ) = ( f 8 + f 6 ) − ( f 5 + f 7 ) 4 h 2 \left(\frac{\partial f}{\partial x}\right)=\frac{f_{1}-f_{3}}{2 h} \\
\left(\frac{\partial f}{\partial y}\right)=\frac{f_{2}-f_{4}}{2 h} \\
\left(\frac{\partial^{2} f}{\partial x^{2}}\right)=\frac{f_{1}+f_{3}-2 f_{0}}{h^{2}} \\
\left(\frac{\partial^{2} f}{\partial y^{2}}\right)=\frac{f_{2}+f_{4}-2 f_{0}}{h^{2}} \\
\left(\frac{\partial^{2} f}{\partial x \partial y}\right)=\frac{\left(f_{8}+f_{6}\right)-\left(f_{5}+f_{7}\right)}{4 h^{2}} \\ ( ∂ x ∂ f ) = 2 h f 1 − f 3 ( ∂ y ∂ f ) = 2 h f 2 − f 4 ( ∂ x 2 ∂ 2 f ) = h 2 f 1 + f 3 − 2 f 0 ( ∂ y 2 ∂ 2 f ) = h 2 f 2 + f 4 − 2 f 0 ( ∂ x ∂ y ∂ 2 f ) = 4 h 2 ( f 8 + f 6 ) − ( f 5 + f 7 )
中间层差分求导（x方向）：
图示：
( ∂ f ∂ σ ) = f 2 − f 4 2 h ( ∂ 2 f ∂ σ 2 ) = f 2 + f 4 − 2 f 0 h 2 ( ∂ 2 f ∂ x ∂ σ ) = ( f 8 + f 6 ) − ( f 5 + f 7 ) 4 h 2 \begin{array}{l}
\left(\frac{\partial f}{\partial \sigma}\right)=\frac{f_{2}-f_{4}}{2 h} \\
\left(\frac{\partial^{2} f}{\partial \sigma^{2}}\right)=\frac{f_{2}+f_{4}-2 f_{0}}{h^{2}} \\
\left(\frac{\partial^{2} f}{\partial x \partial \sigma}\right)=\frac{\left(f_{8}+f_{6}\right)-\left(f_{5}+f_{7}\right)}{4 h^{2}}
\end{array} ( ∂ σ ∂ f ) = 2 h f 2 − f 4 ( ∂ σ 2 ∂ 2 f ) = h 2 f 2 + f 4 − 2 f 0 ( ∂ x ∂ σ ∂ 2 f ) = 4 h 2 ( f 8 + f 6 ) − ( f 5 + f 7 )
中间层差分求导（y方向）：
图示：
( ∂ 2 f ∂ y ∂ σ ) = ( f 3 + f 6 ) − ( f 3 + f 7 ) 4 h 2 \left(\frac{\partial^{2} f}{\partial y \partial \sigma}\right)=\frac{\left(f_{3}+f_{6}\right)-\left(f_{3}+f_{7}\right)}{4 h^{2}} ( ∂ y ∂ σ ∂ 2 f ) = 4 h 2 ( f 3 + f 6 ) − ( f 3 + f 7 )
