1. SIFT算法中一些符號的說明
$I(x,y)$表示原圖像。
$G(x,y,\sigma)$表示高斯濾波器，其中$G(x,y,\sigma) = \frac{1}{2\pi\sigma^2}exp(-(x^2+y^2)/2\sigma^2)$。
$L(x,y,\sigma)$表示由一個高斯濾波器與原圖像卷積而生成的圖像，即$L(x,y,\sigma) = G(x,y,\sigma)\otimes I(x,y)$。一系列的$\sigma_i$，則可以生成一系列的$L(x,y,\sigma_i)$圖像，此時我們把這一系列的$L(x,y,\sigma)$圖像稱為原圖像的一個尺度空間表示。關於尺度空間的知識可以參考：圖像特征提取：尺度空間理論。
$DOG$表示高斯差分（Difference of Gaussians），也可以表示為$D(x,y,\sigma)$，其中$D(x,y,\sigma) = (G(x,y,k\sigma) – G(x,y,\sigma))\otimes I(x,y) = L(x,y,k\sigma) – L(x,y,\sigma)$。
上面特別值得注意的是尺度為$\sigma$的高斯差分圖像由於尺度為$k\sigma$與尺度為$\sigma$的L圖像生成的。$k$為兩相鄰尺度空間倍數的常數。
$O$：高斯金字塔的組數（Octave），其中值得注意的是在實際構建中，第一組的索引可以為0也可以為-1，這個在后面解釋原理。
$S$：高斯金字塔每一組的層數。在實際最開始構建尺度空間圖像，即L圖像的時候，構建了S+3層，一定要把這個S+3與S區分開，為什么是S+3后面分析。
2. 構建高斯差分金字塔
2.1 第一組第一層圖像的生成
很多初涉SIFT的都會被這個問題所困惑，這里要分兩種情況：其一是把第一組的索引定為0；其二是把第一組的索引定為-1。
我們先考慮第一組索引為0的情況，我們知道第一組第一層的圖像是由原圖像與$\sigma_o$（一般設置為1.6）的高斯濾波器卷積生成，那么原圖像是誰呢？是$I(x,y)$嗎？不是！為了圖像反走樣的需要，通常假設輸入圖像是經過高斯平滑處理的，其值為$\sigma_n = 0.5$，即半個像元。意思就是說我們采集到的圖像$I(x,y)$，已經被$\sigma = \sigma_n = 0.5$的高斯濾波器平滑過了。所以我們不能直接對$I(x,y)$用$\sigma_0$的高斯濾波器平滑，而應該用$\sigma = \sqrt{\sigma_0^2 - \sigma_n^2}$的高斯濾波器去平滑$I(x,y)$，即
$$FirstLayer(x,y) = I(x,y)\otimes G(x,y,\sqrt{\sigma_0^2 - \sigma_n^2})$$
其中$FirstLayer(x,y)$表示整個尺度空間為第1組第1層的圖像，$\sigma_o$一般取1.6，$\sigma_n = 0.5$。
現在我們來考慮把第一組的索引定為-1的情況。那么首先第一個問題便是為什么要把索引定為-1。如果索引為0，如上面那種情況所示，整個尺度空間的第1組的第1層圖像已經是由原圖像模糊生成的了，那么也就是說已經丟失了細節信息，那么原圖像我們完全沒有利用上。基於這種考慮，我們先將圖像放大2倍，這樣原圖像的細節就隱藏在了其中。由上面一種情況分析，我們已經知識了I(x,y)看成是已經被$\sigma_n = 0.5$模糊過的圖像，那么將$I(x,y)$放大2倍后得到$I_s(x,y)$，則可以看為是被$2\sigma_n= 1$的高斯核模糊過的圖像。那么由$I_s$生成第1組第1層的圖像用的高斯濾波器的$\sigma = \sqrt{\sigma_0^2 -(2 \sigma_n)^2}$。可以表示為。
$$FirstLayer(x,y) = I_s(x,y)\otimes G(x,y,\sqrt{\sigma_0^2 - (2\sigma_n)^2})$$
其中$FirstLayer(x,y)$表示整個尺度空間為第1組第1層的圖像，$I_s(x,y)$是由$I(x,y)$用雙線性插值放大后的圖像。$\sigma_o$一般取1.6，$\sigma_n = 0.5$。
2.2 尺度空間生成了多少幅圖像
我們知道S是我們最終構建出來的用來尋找特征點的高斯差分圖像，而特征點的尋找需要查找的是空間局部極小值，即在某一層上查找局部極值點的時候需要用到上一層與下一層的高斯差分圖像，所以如果我們需要查找S層的特征點，需要S+2層高斯差分圖像，然后查找其中的第2層到第S+1層。
而每一個高斯差分圖像$G(x,y,\sigma)$都需要兩幅尺度空間的圖像$L(x,y,k\sigma)$與$L(x,y,\sigma)$進行差分生成，這里假設S =3，則我們需要的高斯差分圖像有S+2 = 5張，分別為$G(x,y,\sigma),G(x,y,k\sigma),G(x,y,k^2\sigma),G(x,y,k^3\sigma),G(x,y,k^4\sigma)$。其中的$G(x,y,k\sigma),G(x,y,k^2\sigma),G(x,y,k^3\sigma)$這三張圖像是我們用來查找局部極值點的圖像。那么我們則需要S+3 = 6張尺度空間圖像來生成上面那些高斯差分圖像，它們分別為：$L(x,y,\sigma),L(x,y,k\sigma),L(x,y,k^2\sigma),L(x,y,k^3\sigma),L(x,y,k^4\sigma)，L(x,y,k^5\sigma)$
從上面的分析,我們知道對於尺度空間來說,我們一共需要S+3層圖像來構建出來S+2層高斯差分圖像。所以，如果整個尺度空間一共有O組，每組有S+3層圖像。共O*(S+3)張尺度圖像，如果我們查找OpenCV中的SIFT源碼，則很容易找到如下代碼來說明問題：
pyr.resize(nOctaves*(nOctaveLayers + 3));
上面代碼中的pyr代表了整個尺度空間的圖像，nOctaves為組數，nOctaveLayers即為我們定義的S。
2.3 為什么是倒數第3張
相信你在看很多SIFT算法描述里都這樣寫着，取上一張的倒數第3張圖像隔行采樣后作為下一組的第一張圖像。
答案是為了保證尺度空間的連續性，我們下面來仔細分析。
我們知道對於尺度空間第o組，第s層的圖像，它的尺度為$\sigma = \sigma_o k^{o+s/S}，其中,k =1/2,o\in[0,1,2,\dots,nOctave-1],s\in[0,1,2,\dots,S+2]$。那么我們從第0組開始，看它各層的尺度。
第0組：$\sigma_o \to 2^{1/3}\sigma_0\to 2^{2/3}\sigma_0\to 2^{3/3}\sigma_0\to 2^{4/3}\sigma_0\to 2^{5/3}\sigma_0$
第1組：$2\sigma_o\to 2*2^{1/3}\sigma_0\to 2*2^{2/3}\sigma_0\to 2*2^{3/3}\sigma_0\to 2*2^{4/3}\sigma_0\to 2*2^{5/3}\sigma_0$
我們只分析2組便可以看出，第1組的第0層圖像恰好與第0組的倒數第三幅圖像一致，尺度都為$2\sigma_0$，所以我們不需要再根據原圖來重新卷積生成每組的第0張圖像，只需采用上一層的倒數第3張來降采樣即可。
我們也可以繼續分析，第0組尺度空間得到的高斯差分圖像的尺度為：$\sigma_o\to 2^{1/3}\sigma_0\to 2^{2/3}\sigma_0\to 2^{3/3}\sigma_0\to 2^{4/3}\sigma_0$
而第1組尺度空間得到的高斯差分圖像的尺度為：$2\sigma_o\to 2*2^{1/3}\sigma_0\to 2*2^{2/3}\sigma_0\to 2*2^{3/3}\sigma_0\to 2*2^{4/3}\sigma_0$
如果我們把它們的中間三項取出來拼在一起，則尺度為：$2^{1/3}\sigma_0\to 2^{2/3}\sigma_0\to 2^{3/3}\sigma_0\to 2*2^{1/3}\sigma_0\to 2*2^{2/3}\sigma_0\to 2*2^{3/3}\sigma_0$，正好連續！！這一效果帶來的直接的好處是在尺度空間的極值點確定過程中，我們不會漏掉任何一個尺度上的極值點，而是能夠綜合考慮量化的尺度因子。
2.4 用第i-1層的圖像生成第i層的圖像
值得注意的是，在SITF的源碼里，尺度空間里的每一層的圖像（除了第1層）都是由其前面一層的圖像和一個相對$sigma$的高斯濾波器卷積生成，而不是由原圖和對應尺度的高斯濾波器生成的，這一方面是因為我前面提到的不存在所謂意思上的“原圖”，我們的輸入圖像$I(x,y)$已經是尺度為$\sigma = 0.5$的圖像了。另一方面是由於如果用原圖計算，那么相鄰兩層之間相差的尺度實際上非常小，這樣會造成在做高斯差分圖像的時候，大部分值都趨近於0，以致於后面我們很難檢測到特征點。
基於上面兩點原因（個人認為原因1是最主要的，原因2只是根據實際嘗試后的一個猜想，並無理論依據），所以對於每一組的第$i+1$層的圖像，都是由第$i$層的圖像和一個相對尺度的高斯濾波器卷積生成。
那么相對尺度如何計算呢，我們首先考慮第0組，它們的第$i+1$層圖像與第$i$層圖像之間的相對尺度為$SigmaDiff_i = \sqrt{(\sigma_0k^{i+1})^2 – (\sigma_0 k^i)^2}$，為了保持尺度的連續性，后面的每一組都用這樣一樣相對尺度（SIFT實際代碼里是這樣做的）。這里有一個猜測，比如說尺度為$2\sigma_0$的這一組，第$i$層與第$i+1$層之間的相對尺度計算的結果應該為$\sqrt{(2\sigma_0k^{i+1})^2 – (2\sigma_0 k^i)^2} = 2\cdot SigmaDiff_i$，可是代碼里依然用$SigmaDiff_i$是因為這一層被降維了。
sig[0] = sigma; double k = pow(2., 1. / nOctaveLayers); for (int i = 1; i < nOctaveLayers + 3; i++) { double sig_prev = pow(k, (double)(i - 1))*sigma; double sig_total = sig_prev*k; sig[i] = std::sqrt(sig_total*sig_total - sig_prev*sig_prev); }
3. 特征點的搜索
3.1 搜索策略
斑點的搜索是通過同一組內各DoG相鄰層之間比較完成的。為了尋找尺度空間的極值點，每一個采樣點要和它所有的相鄰點進行比較，看其是否比它的圖像域和尺度域的相鄰點大或小。對於其中的任意一個檢測點都要和它同尺度的8個相鄰點和上下相鄰尺度對應的$9\times2$個點共26個點比較，以確保在尺度空間和二維圖像位置空間都檢測到極值點。也就是，比較是在一個$3\times3$的立方體內進行。
搜索過程從每組的第二層開始，以第二層為當前層，對第二層的DoG圖像中的每個點取一個$3\times3$的立方體，立方體上下層為第一層與第三層。這樣，搜索得到的極值點既有位置坐標（DoG的圖像坐標），又有空間尺度坐標（層坐標）。當第二層搜索完成后，再以第三層作為當前層，其過程與第二層的搜索類似。當S=3時，每組里面要搜索3層。
3.2 子像元插值
上的的極值點的搜索是在離散空間中進行的，檢測到的極值點並不是真正意義上的極值點。下圖顯示了一維信號離散空間得到的極值點與連續空間的極值點之間的差別。利用已知的離散空間點插值到連續空間極值點的方法叫子像元插值。
首先我們來看一個一維函數插值的例子。我們已經$f(x)$上幾個點的函數值$f(-1) = 1,f(0) = 6,f(1) = 5$，求$f(x)$在$[-1,1]$上的最大值。
如果我們只考慮離散的情況，那么只用簡單比較一下，便知最大值為$f(0) = 6$，下面我們用子像元插值法來考慮連續區間的上情況：
利用泰勒級數，可以將$f(x)$在$f(0)$附近展開為：
$$f(x) \approx f(0) + f’(0)x+\frac{f’’(0)}{2}x^2$$
另外我們知道$f(x)$在$x$處的導數寫成離散的形式為$f’(x) = \frac{f(x+1) – f(x)}{2}$，二階導數寫成離散形式為$f’’(x) = f(x+1)+f(x-1)-2f(x)$。
所以，我們可以算出$f(x) \approx 6+2x+\frac{-6}{2}x^2 = 6+2x-3x^2$
求取函數$f(x)$的極大值和極大值所在的位置：
$$f’(x) = 2-6x = 0, \ \ \ \hat{x} = \frac{1}{3}$$
$$f(\hat{x}) = 6+2\times \frac{1}{3} – 3\times(\frac{1}{3})^2 = 6\frac{1}{3}$$
現在回到我們SIFT點檢測中來，我們要考慮的是一個三維問題，假設我們在尺度為$\sigma$的尺度圖像$D(x,y)$上檢測到了一個局部極值點，空間位置為$(x,y,\sigma)$，由上面的分析我們知道，它只是一個離散情況下的極值點，連續情況下，極值點可能落在了$(x,y,\sigma)$的附近，設其偏離了$(x,y,\sigma)$的坐標為$(\Delta x,\Delta y,\Delta \sigma)$。則對$D(\Delta x,\Delta y,\Delta\sigma)$可以表示為在點$(x,y,\sigma)$處的泰勒展開：
$$D(\Delta x,\Delta y,\Delta \sigma) = D(x,y,\sigma)+\begin{bmatrix}
\frac{\partial D}{x} & \frac{\partial D}{y} & \frac{\partial D}{\sigma}
\end{bmatrix}\begin{bmatrix}
\Delta x\\
\Delta y\\
\Delta \sigma
\end{bmatrix}+\frac{1}{2}\begin{bmatrix}
\Delta x &\Delta y & \Delta \sigma
\end{bmatrix}\begin{bmatrix}
\frac{\partial ^2D}{\partial x^2} & \frac{\partial ^2D}{\partial x\partial y} &\frac{\partial ^2D}{\partial x\partial \sigma} \\
\frac{\partial ^2D}{\partial y\partial x}& \frac{\partial ^2D}{\partial y^2} & \frac{\partial ^2D}{\partial y\partial \sigma}\\
\frac{\partial ^2D}{\partial \sigma\partial x}&\frac{\partial ^2D}{\partial \sigma\partial y} & \frac{\partial ^2D}{\partial \sigma^2}
\end{bmatrix}\begin{bmatrix}
\Delta x\\
\Delta y\\
\Delta \sigma
\end{bmatrix}$$
可以將上式寫成矢量形式如下：
$$D(x) = D+\frac{\partial D^T}{\partial x}\Delta x+\frac{1}{2}\Delta x^T\frac{\partial ^2D^T}{\partial x^2}\Delta x$$
令上式的一階導數等於0，可以求得$\Delta x = -\frac{\partial^2D^{-1}}{\partial x^2}\frac{\partial D(x)}{\partial x}$
通過多次迭代(Lowe算法里最多迭代5次)，得到最終候選點的精確位置與尺度$\hat{x}$，將其代入公式求得$D(\hat{x})$，求其絕對值得$|D(\hat{x})|$。如果其絕對值低於閾值的將被刪除。
Vec3f dD((img.at<sift_wt>(r, c + 1) - img.at<sift_wt>(r, c - 1))*deriv_scale, (img.at<sift_wt>(r + 1, c) - img.at<sift_wt>(r - 1, c))*deriv_scale, (next.at<sift_wt>(r, c) - prev.at<sift_wt>(r, c))*deriv_scale); // dD為一階差分矢量Df/Dx float v2 = (float)img.at<sift_wt>(r, c) * 2; float dxx = (img.at<sift_wt>(r, c + 1) + img.at<sift_wt>(r, c - 1) - v2)*second_deriv_scale; float dyy = (img.at<sift_wt>(r + 1, c) + img.at<sift_wt>(r - 1, c) - v2)*second_deriv_scale; float dss = (next.at<sift_wt>(r, c) + prev.at<sift_wt>(r, c) - v2)*second_deriv_scale; float dxy = (img.at<sift_wt>(r + 1, c + 1) - img.at<sift_wt>(r + 1, c - 1) - img.at<sift_wt>(r - 1, c + 1) + img.at<sift_wt>(r - 1, c - 1))*cross_deriv_scale; float dxs = (next.at<sift_wt>(r, c + 1) - next.at<sift_wt>(r, c - 1) - prev.at<sift_wt>(r, c + 1) + prev.at<sift_wt>(r, c - 1))*cross_deriv_scale; float dys = (next.at<sift_wt>(r + 1, c) - next.at<sift_wt>(r - 1, c) - prev.at<sift_wt>(r + 1, c) + prev.at<sift_wt>(r - 1, c))*cross_deriv_scale; Matx33f H(dxx, dxy, dxs, dxy, dyy, dys, dxs, dys, dss); // dD + Hx = 0 --> x = H^-1 * (-dD) Vec3f X = H.solve(dD, DECOMP_LU);
3.3 刪除邊緣效應
為了得到穩定的特征點，只是刪除DoG響應值低的點是不夠的。由於DoG對圖像中的邊緣有比較強的響應值，而一旦特征點落在圖像的邊緣上，這些點就是不穩定的點。一方面圖像邊緣上的點是很難定位的，具有定位歧義性；另一方面這樣的點很容易受到噪聲的干擾而變得不穩定。
一個平坦的DoG響應峰值往往在橫跨邊緣的地方有較大的主曲率，而在垂直邊緣的方向有較小的主曲率。而主曲率可以通過$2\times2$的Hessian矩陣$H$求出：
$$H(x,y) = \begin{bmatrix}D_{xx}(x,y) & D_{xy}(x,y)\\ D_{xy}(x,y) &D_{yy}(x,y) \end{bmatrix}$$
上式中，$D$值可以通過求取鄰近點像素的差分得到。$H$的特征值與$D$的主曲率成正比例。我們可以避免求取具體的特征值，因為我們只關心特征值的比例。令$\alpha = \lambda_{max}$為最大的特征值，$\beta = \lambda_{min}$為最小的特征值，那么，我們通過$H$矩陣直跡計算它們的和，通過$H$矩陣的行列式計算它們的乘積：
$$Tr(H) = D_{xx}+D_{yy} = \alpha +\beta$$
$$Det(H) = D_{xx}D_{yy}-(D_{xy})^2=\alpha \beta$$
如果$\gamma$為最大特征值與最小特征值之間的比例，那么$\alpha = \gamma \beta$，這樣便有
$$\frac{Tr(H)^2}{Det(H)} = \frac{(\alpha+\beta)^2}{\alpha\beta} = \frac{(\gamma+1)^2}{\gamma}$$
上式的結果只與兩個特征值的比例有關，而與具體特征值無關。當兩個特征值相等時，$\frac{(\gamma+1)^2}{\gamma}$的值最小，隨着$\gamma$的增加，$\frac{(\gamma+1)^2}{\gamma}$的值也增加。所以要想檢查主曲率的比例小於某一閾值$\gamma$，只要檢查下式是否成立：
$$\frac{Tr(H)^2}{Det(H)} < \frac{(\gamma+1)^2}{\gamma}$$
Lowe在論文中給出的$\gamma = 10$。也就是說對於主曲率比值大於10的特征點將被刪除。
float t = dD.dot(Matx31f(xc, xr, xi)); //D(\bar{x}) = D + 1/2*dD*\bar{x} contr = img.at<sift_wt>(r, c)*img_scale + t * 0.5f; // 插值得到的極值點的值 if (std::abs(contr) * nOctaveLayers < contrastThreshold) return false; // principal curvatures are computed using the trace and det of Hessian float v2 = img.at<sift_wt>(r, c)*2.f; float dxx = (img.at<sift_wt>(r, c + 1) + img.at<sift_wt>(r, c - 1) - v2)*second_deriv_scale; float dyy = (img.at<sift_wt>(r + 1, c) + img.at<sift_wt>(r - 1, c) - v2)*second_deriv_scale; float dxy = (img.at<sift_wt>(r + 1, c + 1) - img.at<sift_wt>(r + 1, c - 1) - img.at<sift_wt>(r - 1, c + 1) + img.at<sift_wt>(r - 1, c - 1)) * cross_deriv_scale; float tr = dxx + dyy; float det = dxx * dyy - dxy * dxy; if (det <= 0 || tr*tr*edgeThreshold >= (edgeThreshold + 1)*(edgeThreshold + 1)*det) return false;