結合論文[1]和Rob Hess的開源SIFT代碼(發現OpenCV2.3的源碼里也是用的Rob Hess的SIFT代碼)對SIFT算法進行了研究，下面是小結：
在計算機視覺的領域中，圖像匹配是很多問題最重要的一個方面，包括物體和場景識別，通過多幅圖像進行3D重構，立體匹配和運動跟蹤。SIFT特征對於圖像的旋轉和尺度變化具有不變性，對於光照改變和攝像機角度變化具有部分的不變性。SIFT算法生成圖像特征的主要步驟有以下幾個：
（1）尺度空間極值檢測：搜索所有尺度上的圖像位置。通過高斯微分函數來識別潛在的對於尺度和旋轉不變的興趣點。
（2）關鍵點的定位：在每個候選的位置上，通過一個擬合精細的模型來確定位置和尺度。關鍵點的選擇依據於它們的穩定程度。
（3）方向的確定：基於圖像局部的梯度方向，分配給每個關鍵點位置一個或多個方向。所有后面的對圖像數據的操作都相對於關鍵點的方向、尺度和位置進行變換，從而提供對於這些變換的不變性。
（4）關鍵點描述：在每個關鍵點周圍的鄰域內，在選定的尺度上測量圖像局部的梯度。這些梯度被變換成一種表示，這種表示允許比較大的局部形狀的變形和光照變化。
2詳細過程
這里先根據Rob Hess的代碼簡化一個描述SIFT特征的結構體：
struct feature
{
double x; /*!< 在輸入圖像上的x坐標 */
double y; /*!< 在輸入圖像上的y坐標 */
double scl; /*!< 關鍵點的尺度 */
double ori; /*!< 關鍵點的方向 */
int d; /*!< 描述器的長度 */
double descr[FEATURE_MAX_D]; /*!< 描述器 */
int r; /*!< 檢測到關鍵點那一層圖像上關鍵點所在的行 */
int c; /*!< 檢測到關鍵點那一層圖像上關鍵點所在的列 */
int octv; /*!< 檢測到關鍵點那一層圖像所在的組的索引*/
int intvl; /*!< 檢測到關鍵點那一層圖像所在的層的索引 */
double subintvl; /*!< 檢測到關鍵點那一層圖像所在的層的插值 */
double scl_octv; /*!< 檢測到關鍵點那一層圖像所在的組的尺度 */
};
2.1尺度空間表示
圖像的尺度空間表示成一個函數L(x, y, σ)，它是由一個變化尺度的高斯函數G(x, y, σ)與圖像I(x, y)卷積生成的，見Lindeberg.T的論文[2]。即
其中
代表卷積操作，而G(x, y, σ)為：
SIFT算法建議在某一個尺度上的對關鍵點的檢測，可以通過對兩個相鄰高斯尺度空間的圖像相減，得到一個DoG（高斯微分）的響應值圖像D(x, y, σ)。然后，通過對響應值圖像D(x, y, σ)進行局部極大搜索，在位置空間和尺度空間中定位關鍵點。其中：
式中，k為兩相鄰尺度空間倍數的常數。
事實上DoG是對LoG(高斯拉普拉斯)的一個近似。根據熱傳導方程（注：這里t=σ^2）：
對上式進行有限差分運算，得到：
由於常數k-1並不會影響極值點的位置，即關鍵點中心的位置。所以，D(x, y, σ) 是
的近似表示，也即DoG是LoG的近似表示。
接下來是在如何計算DoG。首先是構建圖像的高斯金字塔。將圖像金字塔共分O組，一組稱為一個Octave，每組又分為多層，層間隔數為S，因此有S+3(S+1+2，2代表在上下再各添一層圖像，搜索極值只在中間的S+1層圖像上搜索)層圖像，下一組的第一層圖像由上一組的倒數第三層（如果層索引從0開始，則為第S層）圖像按照隔點采樣得到，這樣做的目的是為了減少卷積運算的工作量。DoG是通過高斯金字塔中的每組上下相鄰兩層的高斯尺度空間圖像相減得到。如下圖：
在高斯金字塔中，σ和o, s的關系如下：
其中，σ0是基准層尺度，o為組Octave的索引，s為組里圖像的層索引。
我們認為初始圖像有一個初始的σ=0.5的高斯模糊（為了防止顯著的混淆效應）。從上式可以看出，（2）式中的k=2^(1/S)。
在最開始建立高斯金字塔時，要預先模糊輸入圖像來作為第0個組的第0層的圖像，這時相當於丟棄了最高的空域的采樣率。因此通常的做法是先將圖像的尺度擴大一倍來生成第0組。我們假定初始的輸入圖像為了抗擊混淆現象，已經對其進行σn=0.5的高斯模糊，如果輸入圖像的尺寸用雙線性插值擴大一倍，那么相當於σn=1.0。因此，實際
上面用到了高斯的半群性質得到，即：
其中t=σ2。
在實際的編程實現中（2）式如下實現：
在D.Lowe的論文和Rob Hess的代碼中，使用的默認參數為：
σn=0.5, σ0=1.6, S=2, O的值通過對圖像的長和寬中較小的一個長度取2的對數再減去2得到（保證最后一組的圖像至少要有4個像素）。
下面是尺度空間表示的一個例子：
2.2尺度空間局部極值檢測（關鍵點的初步探查）
興趣點的初步探查是通過同一組內各DoG相鄰層之間比較完成的。為了尋找尺度空間的極值點，每一個采樣點要和它的相鄰點進行比較，看其是否比它的圖像域和尺度域的相鄰點大或者小。如下圖所示：
上圖中表示×的檢測點和它同尺度的8個相鄰點和上下相鄰尺度對應的9×2個點共26個點進行比較，以確保在尺度空間和二維圖像位置空間都檢測到極值點。搜索過程從每組的第二層開始，也就是從層索引為1的那一層開始，到層索引為S的層結束。
下面是局部極值檢測的結果：
2.3局部極值位置的插值（關鍵點的精確定位）
以上極值點的搜索是在離散空間中進行的，檢測到的極值點並不是真正意義上的極值點。下圖演示了二維函數離散空間得到的極值點與連續空間極值點的差別。利用已知的離散空間點插值得到的連續空間極值點的方法叫做子像素插值（Sub-pixel Interpolation）。
插值的方法是根據泰勒級數展開，具體過程參見[3]，這里只給出結論：
極值點位置的偏移量（位置矢量為）：
極值點的極值為
其中
代表相對於插值中心的偏移量，當它在任一維度上的偏移量大於0.5時（即x或y或σ），意味着插值中心已經偏移到它的鄰近點上，所以必須改變當前關鍵點的位置（加上
,實際編成中對
進行四舍五入取整），同時在新的位置上反復插值直到收斂；也有可能超出所設定的迭代次數或者超出圖像邊界的范圍，此時這樣的點應該刪除。另外，當
的絕對值小於0.03時（假設圖像的灰度值是在0到1.0之間），其響應值過小，這樣的點易受噪聲的干擾而變得不穩定，所以也要刪除。
下面是移除
絕對值較小的關鍵點的結果：
2.4消除邊界響應（關鍵點的進一步篩選）
由於DoG對圖像中的邊緣有比較強的響應值，而一旦特征落在圖形的邊緣上，這些點就是不穩定的點。根據Harris[4]角點可以知道，一個角點在任何方向上平移都應該保證局部窗口內的像素值的劇烈變化，而邊緣上的點沿着邊緣方向移動時局部窗口內的像素值基本沒有什么變化。如下圖所示：
同樣，一個平坦的DoG響應峰值往往在橫跨邊緣的地方有較大的主曲率，而在垂直的方向有較小的主曲率。而主曲率可以通過2×2的Hessain矩陣H求出：
H的特征值與D的主曲率成正比例。可以避免求取具體的特征值，因為我們只關心特征值的比例。令α為較大的特征值，β為較小的特征值。
那么
如果γ為較大特征值與較小特征值之間的比值，α=γβ，這樣便有
上式的結果只與兩個特征的比例有關，而與具體特征值無關。兩個特征值相等時，
最小，隨着γ增大，其值也增加。所以，要想檢查主曲率的比值小於某一閾值γ，只要檢查下式是否成立：
D.Lowe在論文中給出γ=10，對於主曲率比值大於10的特征點將被刪除，否則，這些特征點將被保留。下面是移除邊界響應的結果：
到此為止，特征點的位置已經被確定，假設在D(x,y,σ(o,s))上位置為(xk,yk)處確定了一個特征點，在插值后的偏移Xo=(xo,yo,so)(注：尺度σ的導數在DoG金字塔中是根據同一組的圖像前后差分求得，也就是s+1層圖像減去s-1層圖像再除以2，因此這里用so替代σo)，那么先前給出的描述SIFT特征的結構體將如下被填充：
feat.x = (xk+xo) × 2^o
feat.y = (yk+yo) × 2^o
feat.r = yk;
feat.c = xk;
feat.octv = o;
feat.intvl = s;
feat.subintvl = so;
（注：如果最先開始圖像的尺寸被加倍了，那么feat.x和feat.y要除以2）
2.5關鍵點的尺度計算
這一步計算出關鍵點的尺度，根據2.1節給出的公式(6)，如下：
（注：如果最先開始圖像的尺寸被加倍了，那么feat.scl要除以2）
2.6特征點的方向分配
為了實現圖像旋轉的不變性，需要根據檢測到的特征點的局部圖像結構求得一個方向基准。我們使用圖像梯度的方法求取局部結構的穩定方向。對於已經檢測到的特征點，我們知道該特征點在DoG金字塔中的確切位置和尺度（由feat.r, feat.c, feat.octv, feat.intvl, feat.scl_octv確定），在與尺度相應的高斯圖像L(x,y,σ(feat.octv, feat.intvl))上使用有限差分，計算以特征點為中心，以3×1.5×feat.scl_octv為半徑的區域內圖像梯度的幅角和幅值（模值），如下圖所示，幅角和幅值的計算公式如下：
在完成特征點鄰域的高斯圖像的梯度計算后，使用直方圖統計鄰域內像素的梯度方向和幅值。梯度方向直方圖的橫軸是梯度方向角，縱軸是梯度方向角對應的梯度幅值累加。梯度方向直方圖將0-360度的范圍分為36個柱(bins)，每10度一個柱。直方圖的峰值代表該特征點處鄰域內圖像梯度的主方向，也即該特征點的主方向，如下圖所示：
每個累加到梯度方向直方圖的采樣點的梯度幅值都要進行權重處理，加權采用圓形高斯加權函數，其標准偏差σ為1.5×feat.scl_octv，（上面的局部鄰域半徑通過3σ定理得到）。如下圖所示，由於SIFT算法只考慮了尺度和旋轉不變性，並沒有考慮仿射不變性。通過高斯加權，使特征點附近的梯度幅值有較大的權重，這樣可以部分彌補因沒有仿射不變性而產生的特征點不穩定的問題。
另外，當存在一個相當於主峰值能量80%能量的峰值時，則將這個方向認為是該特征點的輔方向。一個特征點可能會被指定具有多個方向（一個主方向，一個以上輔方向），這可以增強匹配的魯棒性，如下圖所示。實際編程實現中，就是把該特征點復制成多份特征點，並將方向值分別賦給這些復制后的特征點，並且，離散的梯度方向直方圖要進行插值擬合處理，來求得更精確的方向角度值。
下面是特征方向和尺度分配的結果（箭頭指向代表方向，長度代表尺度）：
2.7特征描述子表示（特征點的特征矢量生成）
由SIFT描述子h(x,y,θ)是對特征點附近鄰域高斯圖像梯度統計結果的一種表示，它是一個三維的陣列，但通常在編程實現時只表示成一個矢量（1維向量）。特征描述子與特征點所在的尺度有關，因此，對梯度的求取應在特征點對應的高斯圖像上進行。將特征點附近的鄰域划分成Bp*Bp個子區域，每個子區域的尺寸為m×feat.scl_octv個子像素，其中，m=3，Bp=4。考慮到實際計算時，需要采用雙線性插值，計算的圖像區域為m × feat.scl_octv × (Bp+1)。如果再考慮旋轉因素，那么實際計算的圖像區域應為m × feat.scl_octv × (Bp+1) × √2。如下圖所示：
為了保證特征矢量具有旋轉不變形，需要以特征點為中心，將特征點附近鄰域內(m × feat.scl_octv × (Bp+1) × √2 × m × feat.scl_octv × (Bp+1) × √2)圖像的梯度的位置和方向旋轉一個方向角feat.ori。如下圖所示：
在特征點附近鄰域圖像梯度的位置和角度旋轉后，再以特征點為中心，取一個m × feat.scl_octv × Bp × m × feat.scl_octv × Bp大小的圖像區域，並將它等間隔划分成Bp×Bp個子區域，每個間隔為m × feat.scl_octv像素，在每個子區域內計算8個方向的梯度方向直方圖，繪制每個梯度方向的累加值，形成一個種子點。求與特征點主方向時有所不同，此時，每個子區域的梯度方向直方圖將0-360度的范圍划分為8個方向范圍，每個范圍為45度，這樣，每個種子點共有8個方向的梯度強度信息。由於存在4×4（Bp × Bp）個子區域，所以，共有4×4×8=128個數據，即特征向量的長度feat.d=128。如下圖所示：
同樣，要對特征矢量需要信息高斯加權處理，加權采用標准偏差為m × feat.scl_octv × Bp / 2的標准高斯函數。特征矢量形成后，為了去除光照變化的影響，需要對它們進行歸一化處理。在歸一化處理后，對於特征矢量中值大於0.2的要進行截斷，即大於0.2的值只取0.2，然后，再進行一次歸一化處理，其目的是提高特征的鑒別性。
3疑問
(1)在每次需要確定特征點鄰域半徑大小時，使用的尺度並不是公式(6)確定的尺度，而是去掉所在組序號o計算出的尺度，這是為什么？
(2)前面說到
的絕對值小於0.3的點要被舍棄，在Rob Hess的SIFT代碼中，實際上是
的絕對值小於0.3/S的點被舍棄，為何這個對比度閾值要除以S？
4參考文獻
[1]D.Lowe, "Distinctive Image Features from Scale-Invariant Keypoints", January 5, 2004
[2]T.Lindeberg, "Scale-space theory:A basic tool for analysing structures at different scales"
[3]王永明，王貴錦，“圖像局部不變性特征與描述”，國防工業出版社
[4]C. Harris and M. Stephens (1988). "A combined corner and edge detector". Proceedings of the 4th Alvey Vision Conference. pp. 147–151.
本文出自 “UnderTheHood” 博客，請務必保留此出處http://underthehood.blog.51cto.com/2531780/658350