理解SIFT
<div class="note-container">
<h1 class="title" id="leanote-title">理解SIFT</h1>
<div class="content-container" id="content-container">
<!-- 切换 -->
<div class="tab"><a id="tab-markdown">Markdown</a><a id="tab-html">HTML</a></div>
<textarea id="leanote-content-markdown">## SIFT步骤
找尺度空间极值
可以用DoG找
需要讲理论依据和DoG结合进行推导
s值和$\sigma$值的选取问题
关键点定位
和泰勒展开相关
方向赋值
局部图像描述符
SIFT用于目标识别的步骤
目标图像提取SIFT特征；搭建SIFT特征数据库
用快速kNN将目标特征在数据库中查找
用霍夫变换寻找属于同一单目标的簇(cluster)
用least-square(最小二乘法)验证
SIFT算法步骤
多尺度极值检测(用DoG找)
关键点定位
方向赋值
关键点描述符
尺度空间极值检测
尺度空间定义：
\(
L(x,y,\sigma)=G(x,y,\sigma)*I(x,y)
\)
相关论文已证明，真正具有尺度无关性的算子只有一种：归一化的LoG函数：
\(
\sigma^2\nabla^2G
\)
考虑高斯查分函数DoG:
\(
G(x,y,k\sigma)-G(x,y,\sigma)
\)
对应的尺度空间表达式有：
\(
D(x,y,\sigma)=L(x,y,k\sigma)-L(x,y,\sigma) \\
=(G(x,y,k\sigma)-G(x,y,\sigma))*I(x,y)
\)
考虑热放射等式，有：
\(
\frac{\partial G}{\partial \sigma}
=\sigma\nabla^2G
\)
其证明过程：
$
G(x,y)=\frac{1}{2\pi \sigma2}e{\frac{x2+y2}{2\sigma^2}} \
\frac{\partial G}{\partial x}=-\frac{x}{2\pi\sigma4}e{-\frac{x2+y2}{2\sigma^2}} \
\frac{\partial2G}{\partial2x}=-\frac{\sigma2-x2}{2\pi\sigma6}e{-\frac{x2+y2}{2\sigma^2}} \
\nabla2G=\frac{\partial2G}{\partial2x}+\frac{\partial2G}{\partial^2y} = \frac{x2+y2-2\sigma2}{2\pi\sigma6}e{-\frac{x2+y2}{2\sigma2}} \
\frac{\partial G}{\partial \sigma}=\frac{x2+y2-2\sigma2}{2\pi\sigma5}e{-\frac{x2+y2}{2\sigma2}} \
\therefore \frac{\partial G}{\partial x}=\sigma \nabla^2 G
G(x,y)=\frac{1}{2\pi \sigma2}e{\frac{x2+y2}{2\sigma^2}} \
\frac{\partial G}{\partial x}=-\frac{x}{2\pi\sigma4}e{-\frac{x2+y2}{2\sigma^2}} \
\frac{\partial2G}{\partial2x}=-\frac{\sigma2-x2}{2\pi\sigma6}e{-\frac{x2+y2}{2\sigma^2}} \
\nabla2G=\frac{\partial2G}{\partial2x}+\frac{\partial2G}{\partial^2y} = \frac{x2+y2-2\sigma2}{2\pi\sigma6}e{-\frac{x2+y2}{2\sigma2}} \
\frac{\partial G}{\partial \sigma}=\frac{x2+y2-2\sigma2}{2\pi\sigma5}e{-\frac{x2+y2}{2\sigma2}} \
\therefore \frac{\partial G}{\partial x}=\sigma \nabla^2 G
\(
而考虑\)\frac{\partial G}{\partial \sigma}\(的近似替代公式：
\)
\sigma\nabla^2G=\frac{\partial G}{\partial \sigma} \approx \frac{G(x,y,k\sigma)-G(x,y,\sigma)}{k\sigma-\sigma}
\(
因此有：
\)
G(x,y,k\sigma)-G(x,y,\sigma) \approx (k-1)\sigma2\nabla2G
$
即：用DoG来近似替代LoG算子，得到近似不变的尺度无关特性。相差的倍数是常量没有影响。
当k->1，误差几乎为0.
高斯金字塔的建立
octave：采样率相同、尺寸相同，模糊度逐渐增强的一组相邻图像
level：octave中的每张图像，叫做一个level
每个octae内是持续高斯模糊得到；每个octave的第一个levle通过降采样得到。每个octave内有s次模糊。
\(\sigma\)的选取和\(s\)选取，不能太小，否则不准确；太大又会消耗太多计算资源，得不到实时效果。Lowe的论文实验表明，选取\(\sigma=1.6\)，\(s=3\)是合适的。
不过为了充分利用输入的图像，可以考虑在生成octave前做一次升采样：认为照相机采集的图像本身就是模糊的，其\(\sigma=0.5\)，利用线性插值得到\(\sigma=1.0\)的图像，然后用于生成octave
疑问l升采样后的图像，有没有再次做模糊操作得到octave的第一张？还是说直接作为octave的第一个level？
DoG金字塔
每两个高斯金字塔level之间做差值就得到DoG的一个level。
为了得到所有极值，因为高斯金字塔的首尾两层不会被带入计算，因此在首尾各添加一层，最上面再加一层
这段还是不够清晰，需要修改
局部极值检测
在DoG金字塔中，检测空间中每个像素点周边的26个点。如果是极值，那么继续往金字塔高层进行检测。
极值检测需要考虑检测的策略：检测到的极值点数量会影响到后续的特征检测数量，特征越少则匹配的结果越少，但是往往是正确率较高的匹配；特征越多则匹配的结果越多，但是正确率往往不高。前者更多地考虑了效率，后者更多地考虑了完整性，需要在效率和完整性之间权衡。因而，要选取合适数量的特征数量，在这之前需要选取合适数量的极值点，这就需要考虑采样频率的问题。（实际上，集中在一起的极值点对于小的扰动非常不稳定）。
采样频率的问题，包括尺度上的采样频率、空间域上的采样频率两个方面。
Lowe通过实验得到一些测量结果，从结果的图表中得到\(s\)值和\(\sigma\)值的经验值。Lowe的实验图像“有32张，包括室外场景、人脸、航拍图片、工业图像（这些图像的来源领域对结果几乎没有影响）”。但是本人认为，仅仅从32张图像就得到相关的经验参数，并不很具有说服力，如果有时间，不妨使用更多的图片进行相关测定。
尺度上的采样频率问题，就是考虑图像金字塔中每个octave内的尺度空间数量s。实验中尺度图像由随机的旋转和0.2~0.9之间的随机尺度变换得到，同时也添加了1%的噪声（使用平均分布为每个像素点增加一个随机数）。
实验结果为：随着尺度数量s的增加，同一位置的重复匹配率“先增加后减小”，同一极值点对应的特征描述符在数据库中被正确搜索到的比率，也是“先增加后减小”。折可以被解释为：随着尺度数量s的增加，尽管极值点数量增加，但是过了一定限度后增加的极值点都是低于平均稳定程度的，因而反倒降低了匹配的重复率和正确率。对应于Lowe的实验结果，选取\(s=3\)时，正确匹配率最高。而随着尺度数量s的增加，极值点数量也增加，并且正确匹配的极值点的总数也是增加的；而通常目标识别通常更依赖于正确匹配的极值点数量而不是正确匹配率，因此看起来使用更大的尺度空间数s会得到更好的效果。但是，这会增加计算消耗。因而，Lowe的论文中权衡了正确率和计算消耗，选取\(s=3\)作为每个octave内的尺度数量。也就是说，仅仅检测那些稳定的极值点，同时又能够避免大量的计算消耗，还能够得到很好的正确率。
空间域上的采样频率问题，就是高斯模糊的参数\(\sigma\)，即平滑尺度。考虑到极值可能随机地出现在一起，也需要类似前面一种情况考虑采样频率和检率的问题。实验数据表明，随着\(\sigma\)的增加，同一特征点在变换后的图像中被重复检测到的百分比、极值点对应的特征描述符在数据库中被查找到的正确率，都呈现增加趋势。但是\(\sigma\)增加后计算量也增加，权衡考虑算法效率和正确率，使用\(\sigma=1.6\)。（同样地，这显然也是Lowe的32张图像的测试结果带来的经验数据，个人认为这个数据要根据不同的情况可以适当修改）。
同时也可以考虑图像预先平滑操作，这样能有效地忽略最高空域频率（这句其实并不理解），其做法是：认为相机拍摄得到的原始图像的模糊度为\(\sigma=0.5\)，通过双线性差值得到\(\sigma=1.0\)模糊尺度的图像。这就意味着，在创建尺度空间的第一个octave之前，需要做一次平滑操作（也就是：原图升采样得到\(\sigma=1.0\)-->做第一次预先处理的平滑得到\(\sigma=1.6\)的尺度图像作为第一个octave的第一个level）。
这里发现对于“采样”并不很理解，需要找书仔细看下
确定关键点的精确位置
前面得到的极值，只是离散的图像金字塔中不同尺度下的极值点。真正的尺度空间应当是连续的，而不是离散的，图像金字塔中的极值只是真正极值的近似，还需要进一步的精确检测，得到真正的极值点和极值。
使用泰勒公式对DoG函数进行展开，并求导，解出导数为0的点，得到的解就是“前面极值检测中检测到的极值采样点，与真正的尺度空间中的极值点的偏差（这个偏差是向量(x,y,\(\sigma\))）”。即：
\(
定义向量\vec{x} = (x, y, \sigma)
对D函数，泰勒展开有：D(\vec {x})=D+\frac{\partial D^T}{\partial \vec {x}}\vec {x} + \frac{1}{2}\vec {x}^T\frac{\partial^2 D}{\partial \vec{x}^2}\vec {x}
求解得到：\vec {x_0} = - \frac{\partial^2 D^{-1}}{\partial \vec {x}^2}\frac{\partial D}{\partial \vec {x}}
\)
\( 其中使用到对于向量求导的计算，包括： \\ \frac{\partial \vec {x}}{\partial \vec {x}}=单位矩阵E \\ \frac{\partial \vec{x}^T \vec{x}}{\vec{x}} = 2 \vec{x} \)
将\(\vec{x_0}\)带入\(D(\vec{x})\)的表达式，得到对应的极值：
\(
D(\vec{x}) = D + \frac{1}{2}\frac{\partial D^T}{\partial \vec{x}}\vec{x}
\)
实验中忽略\(D \lt 0.03\)的极值点。
这样，得到的\(\vec{x}\)表示采样点的偏离量。此时，keypoint的表示中包括了：位置(x,y),尺度\(\sigma\)。后面还会包括方向。
消除边界响应
还需要去除边界上的干扰，需要考虑主曲率。 通过计算黑塞矩阵能得到特征值之间的关系，而主曲率和特征值之间是成倍的关系，因此只需要考虑特征值之间的关系。
假设黑塞矩阵算出的特征值中，大的为\(\alpha\)，小的为\(\beta\)，而：
$
H=
\begin{matrix}
D_{xx} & D_{xy} \
D_{xy} & D_{yy}
\end{matrix}
$
\(
Tr(H)=D_{xx}+D_{yy} = \alpha + \beta \\
Det(H)=D_{xx}D_{yy}-(D_{xy})^2 = \alpha \beta
\)
显然，\(\frac{Tr(H)}{Det(H)}\)是增函数，因此只需要保证：
\(
\frac{Tr(H)^2}{Det(H)} < \frac{(r+1)^2}{r}
\)
就能够保证主曲率低于设定值。本paper中，取r=10.
方向赋值
为保证后续算出的SIFT特征向量具备旋转不变形，考虑keypoint的梯度方向角度\(\theta\)，当keypoint及其邻域内的像素都旋转\(\theta\)角度后，再进行特征向量的计算，则算出的特征向量就具备了旋转不变性。
为计算keypoint的梯度方向，需要考虑它的邻域内每个像素的梯度方向，合成一个局部范围内的方向。具体步骤包括两步：
１）计算每个点的梯度
首先计算当前尺度下所有像素点的梯度方向，以备后续使用。包括梯度的幅值\(m\)和角度\(\theta\)：
\(
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}(\frac{L(x,y+1)-L(x,y-1)}{L(x+1,y)-L(x-1,y)})
\)
2）计算keypoint的方向直方图
根据keypoint邻域内的点的梯度，计算出keypoint的方向。
假设当前尺度图像的尺度为\(\sigma\)，则使用keypoint的\(３×1.5\sigma\)为半径的圆形高斯窗口范围内的像素点，将它们的梯度幅值\(m(x,y)\)和对应位置的高斯模板中的权值\(w(x,y)\)对应相乘，得到邻域内的点\((x,y)\)的梯度方向的权重；然后按照这个权重进行圆形邻域内的方向直方图统计，按照每\(10°\)一个柱进行统计，一共有36个柱可供选择，覆盖了\(360°\)范围。
统计好的直方图，大约统计了80个的点的方向，会在直方图中形成至少一个极大值。直方图中的最大值作为keypoint的主方向，其他极大值若其值大于等于最大值的80%，则作为辅方向也被选出来，这样原来由位置\((x,y)\)和尺度\(\sigma\)构成的keypoint的表示，现在增加了一项梯度方向\(gradient\quad direction\)，即原有keypoint+主方向构成一个新的keypoint，原有keypoint+辅方向构成另一个新的keypoint。辅方向带来了新增的keypoint，能增加后续匹配的稳定性。
高斯圆形窗口半径的确定
Lowe的论文中提到，高斯圆形窗口的半径，是当前尺度\(\sigma\)的1.5倍。而实际上在空域操作时，高斯模板通常使用\(3*3\)的规格，意味着覆盖\([-3\sigma,3\sigma]\)的范围，概率上覆盖了99.7%的区域；这里提到高斯窗口权重是1.5倍的尺度值，那么就意味着高斯窗口的半径是\(3*1.5\sigma=4.5\sigma\),也就是高斯窗口的直径（即边长）是\(2*4.5\sigma=9\sigma\)。考虑__尺度__本身的含义，表示度量的单位，那么可以认为在尺度图像\(L\)中，一个\(\sigma\)表示一个像素，那么此时的高斯模板就是使用９×９的大小，并且是圆形的区域了。
使用高斯窗口的原因是，考虑不同距离的点对于keypoint的影响力是不同的。
为什么使用当前尺度\(\sigma\)的１.5倍作为高斯圆形窗口的权重呢？即：为什么用\(\sigma_1=1.5\sigma\)作为高斯模糊的参数，而不是使用\(\sigma\)作为高斯模糊的参数？这是因为考虑到了后续要进行方向旋转，每个keypoint会围绕其算出来的这个方向转动，考虑到平面内旋转对称性，keypoint的一个方向区域在平面内转动所覆盖的区域，可以看做转动范围是\([0°,45°]\)。这样一来方向区域中距离中心点(keypoint)最远的点，即矩形的４个角上的像素点，在最极端情形下会转到keypoint的正上方，那么此时的高斯圆形窗口的半径变成了原来的\(\sqrt 2\)倍。作为一个近似替代，使用1.5替代\(\sqrt 2\)，因此使用了\(1.5\sigma\)作为高斯圆形窗口的权重，这是考虑后续旋转的极端情况后做出的设定。
局部图像描述符
前面一步计算了新的keypoint，使得原有的keypoint都具备了梯度方向，目的在于构造出具备旋转不变性的局部特征向量，即通过旋转keypoint邻域内的每个点，旋转角度为keypoint在上一步算出的角度，使得邻域内的点所具备的性质将与旋转角度无关，个人认为这类似于在计算方差时每个点减去均值的做法。然后，再考虑邻域内每个旋转后的点的梯度方向，这些方向按照一定的规则进行合成，就构成了SIFT特征向量。其具体步骤为：
１）keypoint及其邻域内像素点，按照keypoint的方向角\(\theta\)进行旋转
２）选定合适的keypoint descriptor窗口大小，用高斯圆形窗口作为权重，将keypoint邻域内的子区域分别计算方向直方图。
这次的方向直方图，按照\(45°\)为柱的划分单位，\(360°\)被划分为８个柱；然后，每个\(4*4\)像素区域，按照对应的高斯圆形窗口区域分配的权值，将梯度幅值进行直方图统计。
此时的\(4*4\)窗口区域，只是高斯圆形区域内的一个子区域。
将关键点附近的邻域划分为\(d*d\)个子区域（实验结果表明\(d=4\)效果最好），每个子区域作为一个种子点，每个种子点有８个方向。每个子足浴的大小于关键点方向分配时相同，即每个区域有\(3\sigma\_oct\)个子像素，为每个子区域分配边长为\(3\sigma\_oct\)的举行区域进行采样（子像素实际用边长为\(\sqrt(3\sigma\_oct)\)，并且采样点宜多不宜少）。考虑到世纪计算时，需要采用双线性差值，所需图像窗口边长为\(3\sigma\_oct*(d+1)\)。在考虑到旋转因素，则实际所需图像半径为:
\(
radius = 3\sigma\_oct*\frac{d+1}{2}*\sqrt 2
\)
计算结果四舍五入。因为实验测试表明d=4效果最好，那么半径\(radius\)的最佳效果值约为18（这和Lowe原文中的示意图吻合得比较好：原图是d=2的情形，扩大一倍后d=4得到最好的效果，对应的邻域边长是18）。
也就是说，在Lowe使用的测试数据图片下，每个keypoint对应有４个子区域从属于它，每个子区域用来生成８个方向（每个像素点的梯度幅值乘以高斯圆形模板的权重值，作为直方图权重进行方向统计），那么每个keypoint有\(4*4*8＝128\)个方向，这128个方向（准确说是方向幅值）就构成一个SIFT特征向量。
SIFT用于匹配和识别
总算把SIFT特征向量的生成过程搞清楚了，后续通过为每张图片生成若干个SIFT特征，那么多张图片之间的匹配就通过SIFT特征向量的匹配实现。常用\(nearest\ neighborhood\)这类算法进行匹配，另外需要其他算法进行微调，一般过程如下:
目标图像提取SIFT特征；搭建SIFT特征数据库
用快速kNN将目标特征在数据库中查找
用霍夫变换寻找属于同一单目标的簇(cluster)
用least-square(最小二乘法)验证
实际使用中为了达到实时性，kNN算法会被替换为合适的、优化过的邻近匹配算法。最容易想到的是k-D树，但当向量维数大于10时依然很慢，Lowe等人对其进行优化，修改为按照离搜索距离的大小作为顺序进行搜索，并且当搜索到若干数量(200个)后停止搜索，实验表明此时仅仅丢失５％的正确结果，却能提升两个数量级的时间效率。
个人认为，后续可以考虑将SIFT特征向量表示为二进制特征向量，通过二进制计算能实现最快速度的计算，这需要先想清楚二进制特征的结构和运算结果的整理。
