SIFT(Scale-Invariant Feature Transform)是一種具有尺度不變性和光照不變性的特征描述子，也同時是一套特征提取的理論，首次由D. G. Lowe於2004年以《Distinctive Image Features from Scale-Invariant Keypoints[J]》發表於IJCV中。開源算法庫OpenCV中進行了實現、擴展和使用。
本文主要依據原始論文和網絡上相關專業分析，對SIFT特征提取的算法流程進行簡單分析。由於涉及到的知識概念較多，本人能力經驗有限，如果有錯誤的地方麻煩幫助指出，大家一起維護文檔的正確性；有一些內容和圖片來自網絡，文章最后會給出原文地址，侵刪。
一、本文結構
依據原始論文，SIFT算法流程分為4個部分：
1、DoG尺度空間構造（Scale-space extrema detection）
2、關鍵點搜索與定位（Keypoint localization）
3、方向賦值（Orientation assignment）
4、關鍵點描述（Keypoint descriptor）
此外再添加第5個部分，用於說明一些需要特殊說明的重要內容。類似游戲攻略的主線和分線。
以上內容分別對應二至六章，第七章用來放置一些參考文獻。
二、DoG尺度空間構造（Scale-space extrema detection）
2.1"尺度"簡介
圖像處理中的"尺度"可以理解為圖像的模糊程度，類似眼睛近視的度數。尺度越大細節越少，SIFT特征希望提取所有尺度上的信息，也就是無論圖像是否經過放大縮小都能夠提取特征。這種思考，是和人的生理特征類似的，比如我們即使是在模糊的情況下仍然能夠識別物體的種類，生理上人體對物體的識別和分類和其尺度（模糊程度）沒有直接關系。
高斯函數在圖像處理中得到了廣泛的運用，源於其本身具有的優良特性（高斯函數的優良特性放在6.1節）；在圖像處理中對圖像構建尺度空間的方法,是使用不同sigma值的高斯核對圖像進行卷積操作(注意，
sigma就是“尺度”）。
是高斯核，其數學表述為：
整個高斯卷積表示為：
其中
是原圖像，*是卷積符號，
對應尺度下的尺度圖像，這里簡單了解一下，我們看一下在OpenCV中實現不同sigma的高斯卷積的代碼和結果
是原圖像，*是卷積符號，對應尺度下的尺度圖像，這里簡單了解一下，我們看一下在OpenCV中實現不同sigma的高斯卷積的代碼和結果
Mat matSrc
= imread(
"e:/template/lena.jpg");
pyrDown(matSrc,matSrc);
imshow( "原始圖",matSrc);
Mat matDst;
GaussianBlur(matSrc,matDst,Size( 0, 0), 1);
imshow( "sigma是1",matDst);
GaussianBlur(matSrc,matDst,Size( 0, 0), 3);
imshow( "sigma是3",matDst);
GaussianBlur(matSrc,matDst,Size( 0, 0), 5);
imshow( "sigma是5",matDst);
pyrDown(matSrc,matSrc);
imshow( "原始圖",matSrc);
Mat matDst;
GaussianBlur(matSrc,matDst,Size( 0, 0), 1);
imshow( "sigma是1",matDst);
GaussianBlur(matSrc,matDst,Size( 0, 0), 3);
imshow( "sigma是3",matDst);
GaussianBlur(matSrc,matDst,Size( 0, 0), 5);
imshow( "sigma是5",matDst);
其中左上角為原圖，右上角為sigma=1尺度下的圖像，下行左右分別是sigma=3,sigma=5的圖像，尺度越大，圖像越模糊，這是直觀印象。但是即使
sigma=5，你也能夠勉強認出這是一個人像；如果你的頭腦中有lena這個模型，見過lena的圖片，你就能夠很快識別出這是lena的圖片變模糊的結果。這就是基本上一個人腦的識別流程。“尺度”就是這些。
2.2"圖像金字塔"簡介
圖像金字塔也是圖像處理中的常用的概念，一般通過對原始圖像逐步進行下（降）采樣（每隔2個像素抽取一個像素）得到x和y都為原始圖像1/2的縮小圖像。
對於二維圖像，一個傳統的金字塔中，每一層圖像由上一層分辨率的長、寬各一半，也就是四分之一的像素組成；注意金字塔上一個上小下大的東西，它的層數是由下向上來數的，底層叫做第0層，然后往上依次為1、2、3層。
2.3“高斯金字塔”
現在我們已經知道，金字塔就能夠用來解決大小的變化；而不同尺度的gauss核用來解決模糊程度的變化。聰明的計算機工程師想到把這兩者結合起來。
因為有高斯和金字塔嘛，直來直往的就叫做“高斯金字塔”吧。它的模型應該是看上去這樣的：大體看上去是一個金字塔，但是每層都是“復合”結構，是有若干張不同尺度的處理結果疊加起來的。
按照論文里面的叫法，金字塔每層多張圖像合稱為一組（Octave/
'ɔktiv/
），每組有多張（也叫層Interval/
'intəvəl/
）圖像。
高斯金字塔在多分辨率金字塔簡單
降采樣基礎上加了高斯濾波
，但是具體來做的時候有細節在里面：
令每一層金字塔(octave)中尺度變化范圍為
,每層金字塔測量s個尺度。
那么第t層金字塔的尺度范圍就是
,第
層金字塔的第一幅圖像由上層塔中
尺度下的下（降）采樣的到，下采樣比例為0.5。
降采樣時，金字塔上邊一組圖像的第一張圖像（最底層的一張）是由前一組（金字塔下面一組）圖像的倒數第三張下（降）采樣得到。這里有些小混亂，聰明如你一定會提出一些問題，我們把這個叫做“尺度變化連續性”，並放在6.2中結合OpenCV實現的代碼進行講解。
,每層金字塔測量s個尺度。
那么第t層金字塔的尺度范圍就是
,第
層金字塔的第一幅圖像由上層塔中
尺度下的下（降）采樣的到，下采樣比例為0.5。
降采樣時，金字塔上邊一組圖像的第一張圖像（最底層的一張）是由前一組（金字塔下面一組）圖像的倒數第三張下（降）采樣得到。這里有些小混亂，聰明如你一定會提出一些問題，我們把這個叫做“尺度變化連續性”，並放在6.2中結合OpenCV實現的代碼進行講解。
2.4“高斯差分金字塔”
所謂圖像差分，簡單的說就是兩幅圖像不一樣的地方；差分運算就是對兩個圖像做一個相減的運算。之所以這樣做，是因為拍攝連續動作圖像序列，或者是同一場景下不同時刻的圖像，其絕大部分都是相同點。直接相減留下來的就可能特征。
在圖像增強中有一種非常簡單的特征提取方法，就是用原始圖像減去經過高斯模糊后的圖像，能夠得到特征信息；從信號理論這塊來說明，模糊后的圖像保存的是低頻信息，原始圖像減去
低頻信息留下高頻信息。
前面已經有了高斯金字塔，下面可以進一步計算高斯差分金字塔。再次強調：差分就是相減
它有一個比較有名的名字
Difference of Gaussian，簡稱
dog（和狗這種生物沒有什么關系
）。如果對於前面的尺度和金字塔的概念比較清楚，那么一定能夠很快接受下面這幅DOG圖：
或者縮小一些，原始論文中的這張圖看的清楚一些，但是實際上是一個意思：
原文中特別支持這里DoG近似等於尺度歸一化的高斯拉普拉斯算子
,而尺度歸一化的高斯拉普拉斯算子相較於其他角點檢測算子，如梯度，Hessian或Harris焦點特征能夠更穩定的圖像特征。這個算子光是看名字就很厲害，不但有高斯還有
拉普拉斯最后還歸一化了。但這塊內容和算法主線關系不大，我們放到6.3中進行說明吧。
,而尺度歸一化的高斯拉普拉斯算子相較於其他角點檢測算子，如梯度，Hessian或Harris焦點特征能夠更穩定的圖像特征。這個算子光是看名字就很厲害，不但有高斯還有
拉普拉斯最后還歸一化了。但這塊內容和算法主線關系不大，我們放到6.3中進行說明吧。
到目前位置，第一節DoG尺度空間構造結束了，希望你能夠知道什么事DoG尺度空間。是不是內容很多？這里還只是原理主要流程簡析，我們還有幾個比較重要的知識點要在第六章中具體說明；在變成算法的過程中，圖像處理工程師們還做出了更多工作，比如這些層取多少比較合適？有沒有簡化算法？這些都是可以拿來說事的，因為算法沒有效率就和沒有這個理論差不多，這一點的觀點還是要建立。
二、關鍵點搜索與定位（Keypoint localization）
2.1初略尋找
像“掃雷”游戲一樣，除了邊緣點外，每個點有8個相聯的位置（最多有8個地雷）；那么在已經建立起來的DoG空間中，除了邊緣的點，每個點有9*3-1 = 26 個領域。在這個26個領域里面尋找最大最小的點，認為可能是關鍵點。如下圖
標注的點要和周圍26個綠色的點比較。
標注的點要和周圍26個綠色的點比較。
在做卷積運算的時候，專門要考慮邊緣問題，也有很多種解決方法。對於這里存在邊緣問題，如何來處理了？論文中采用了比較直接的方法，就是邊緣的都不要。反過來，我們也需要知道對於
某一層來說，如果我們想獲得s個極值，那么一共需要s+2個dog值（最前面和最后面多兩個），而由於1個dog至少需要2個gauss尺度才能夠計算出來，所以結論是如果想獲得s個極值，那么需要s+3個gauss（這也是為什么“高斯金字塔”在降采樣的時候，選用第3個層來做抽值運算）。下面這幅圖能夠幫助你理解這個概念。
從27個點中找一個極值點，這是比價簡單的運算；所以這里的特征點初略的尋找是比較簡單的，我們很快能從DoG空間中粗略的找到了關鍵點。
2.2細化尋找
下一步是細化，要去除干擾就要分別處理這兩種情況。
2.2.1、較小的極值
對關鍵點進行3D二次函數的擬合(可以類比最優化里牛頓法的做法)，然后對函數求極值。
將每個關鍵點的一段函數
進行泰勒二次展開:
注意，泰勒公式是最為基本的連續函數量化（離散化）的相似公式。
其中x是相對於關鍵點的距離。
和
分別表示函數在關鍵點的取值和一階導數，二階導數。這個展開就是在關鍵點周圍擬合出的3D二次曲線。學導數的時候學過，當導數為0的時候原函數取極值。那么這里二次泰勒的導數為（這種怪異的函數如何求導？放在6.6節,讀正文的時候你就先接受這個結果就可以了，但應該知道這里的這些數學問題）
此時的極值為
，
當這個極值過小的時候，被去除。
原始論文中剔除
的點。
的點。
2.2.2、邊緣噪聲
Hessian矩陣的特征值和特征值所在特征向量方向上的曲率成正比（因為我做過基於hessian方面的增強，所以對於這個方面有一些多的了解，放在6.4中，我們可以多畫幾張圖）。我們希望使用曲率的比值剔除邊緣點，就可以使用Hessian的特征值的比值替代曲率的比值。
計算特征值；以及計算特征值的比例
Hessian矩陣定義如下：
可以由二階差分計算得到，那么接下來就可以計算特征值的比值了。
先等一等，假設兩個特征值分別為
,看下面兩個式子：
其中
分別是矩陣的跡和矩陣的行列式，這樣發現特征值的比值可以通過這兩個兩計算得到，而不需要特征值分解，這樣簡單多了。
分別是矩陣的跡和矩陣的行列式，這樣發現特征值的比值可以通過這兩個兩計算得到，而不需要特征值分解，這樣簡單多了。
假設
為較大特征值，且
,則
為較大特征值，且
,則
當
是
最小，所以當
越大時，對應的
越大。所以我們要將
的點剔除就相當於將
的點剔除。論文中。
但是，原論文里面說這個比率理論來自另一片論文，這里計算“跡”和“行列式”就一定比計算矩陣特征值本身更為簡單嗎？我們放在6.5中進一步研究。這里我們需要知道進一步調用並計算hessian特征值，我們得到了特征點的精細提取。那么第二部分結束了。
這個地方，需要看到，實際上這里出現的兩處“精細篩選”的算法是有相互的關系的：它們一個是去掉了特征點中不是非常“鋒利”的特征點、一個是去掉了屬於“邊緣”的特征點。目的是統一的，都是挑選出更為合適的特征點。
三、方向賦值
在DoG空間中已經找到並細化出若干關鍵點，現在這些特征點還是存儲在類似Mat這樣結構里面的（實際上現在已經有很多可以叫做“特征點”的點了）。但原論文認為還不夠，還需要考慮“旋轉適應性”，就是同樣一個特征，旋轉后了還能夠被認為是同樣的特征。為了解決這個問題，首先就是需要表示“方向”。
在圖像處理中，最直觀的是采用“梯度”表示方向，比如sobel算子等，這里也是類似的方法。
計算圖像局部方向的基本方法，那就是定義啦，理解這個公式應該沒有問題
（注意原論文里這里有個1.5倍尺度重新運算的部分，和主線關系不大）得到這個角度以后，
在完成關鍵點鄰域內高斯圖像梯度計算后，使用直方圖統計鄰域內像素對應的梯度方向和幅值。
梯度方向直方圖的橫軸是梯度方向角，縱軸是剃度方向角對應的梯度幅值累加值。梯度方向直方圖將0°~360°的范圍分為36個柱，每10°為一個柱。下圖是從高斯圖像上求取梯度，再由梯度得到梯度方向直方圖的例圖（注意這里的直方圖是8個柱的簡化版本）。
既然有了直方圖，那么就有高低了（兩個一樣高的這種細節就不要考慮了）。直方圖峰值代表該關鍵點處鄰域內圖像梯度的主方向，也就是該關鍵點的主方向。在梯度方向直方圖中，當存在另一個相當於主峰值 80%能量的峰值時，則將這個方向認為是該關鍵點的輔方向。
所以一個關鍵點可能檢測得到多個方向，這可以增強匹配的魯棒性。Lowe的論文指出大概有15%關鍵點具有多方向，但這些點對匹配的穩定性至為關鍵。
獲得圖像關鍵點主方向后，每個關鍵點有三個信息(x,y,σ,θ)：位置（x,y)、尺度、方向。由此我們可以確定一個SIFT特征區域。通常使用一個帶箭頭的圓或直接使用箭頭表示SIFT區域的三個值：中心表示特征點位置，半徑表示關鍵點尺度,箭頭表示主方向。具有多個方向的關鍵點可以復制成多份，然后將方向值分別賦給復制后的關鍵點。如下圖(這圖認真看看要懂）
那么方向也得到了。
第四部分：關鍵點描述
為找到的關鍵點即SIFT特征點賦了值，包含位置、尺度和方向的信息。接下來的步驟是關鍵點描述，即用一組向量將這個關鍵點描述出來。圖像處理必須考慮相關性，也就是你描述一個點，還必須考慮到它周圍的情況。所以這個描述子不但包括關鍵點，也包括關鍵點周圍對其有貢獻的像素點。
特征描述子與關鍵點所在尺度有關，因此對梯度的求取應在特征點對應的高斯圖像上進行。將關鍵點附近划分成d×d個子區域，每個子區域尺寸為mσ個像元（d=4，m=3，σ為尺特征點的尺度值）。考慮到實際計算時需要雙線性插值，故計算的圖像區域為mσ(d+1)，再考慮旋轉，則實際計算的圖像區域為
在子區域內計算8個方向的梯度直方圖，繪制每個方向梯度方向的累加值，形成一個種子點。
與求主方向不同的是，此時，每個子區域梯度方向直方圖將0°~360°划分為8個方向區間，每個區間為45°。即每個種子點有8個方向區間的梯度強度信息。由於存在d×d，即4×4個子區域，所以最終共有4×4×8=128個數據，形成128維SIFT特征矢量。簡單地看來，就是將第3節的部分擴展為一個更大的矩陣。128維用來表示一個特征點，那是相當大的數據量，這也就是為什么sift計算緩慢的原因吧。
好了，這樣我們已經找到並且表示出來了sift特征點，下面就是比較、存儲等后續操作，相比較原理來說較為次要。后期我們結合OpenCV代碼具體分析的時候進行細化。
6、輔助資料
其實我這篇博客專門將這些輔助的資料挑出來單獨講，一方面是前面的那些系統的資料，網絡上已經有很多寫的非常好的，我也參考了許多，其實相同的地方很多；另一方面我在“輔助資料”這塊融入了許多自己項目實現的一些經驗和整理，希望它能夠成為獨立主要內容，但是又能夠幫助理解主要內容的一部分。
6.1、Gauss的優良特性
什么是Gauss？如果在一維空間里面，基本上能夠理解Gauss就是“正態分布",它天生具有許多優良特性，而且並不是僅僅在圖像處理這塊才體現出來的。
如果整理一下：
1. Lindeberg等人已證明高斯卷積核是實現尺度變換的唯一變換核，並且是唯一的線性核 。
2.高斯函數具有旋轉對稱性，用其對圖像進行平滑運算時在各個方向上的平滑程度相同，因此后續邊緣檢測等操作中不會偏袒某一方向上的圖像的細節。
3.隨着離高斯模板中心點越遠，權值越小，這使得高斯濾波器比起普通的平滑濾波器更能更好地保留圖像細節。如果距離越遠的點權值越重的話，那圖像就會失真。
4.高斯函數有個sigma參數可以非常容易地調節高斯函數對高低頻信息的保留程度。sigma參數越大，高斯函數的圖譜就越低矮平緩，表現在頻譜上就是頻帶越寬，平滑程度高；反之，sigma參數越小，平滑程度越低。而且二維高斯函數的取值半徑（即卷積核大小）越大，平滑程度越高。
舉其中第二點來具體說明，我之前在實現Frangi算法的時候了解到，如果你想求整個圖像的Hessian矩陣，然后再做個Gauss平滑什么的，就可以換成
先”通過先對高斯函數進行偏導操作，然后利用這個核對原圖像進行卷積求解（
圖像的二階偏導數可以通過使用多尺度高斯核的二階偏導數對原圖像的卷積獲得），因為：
我們知道，OpenCV對卷積操作能夠通過FFT進行加速計算，這樣的話，本來是要算整個圖像的偏導的（想想800*600個像素的偏導），現在只要求一個5*5的偏導就可以，速度上就很適合，而且代碼實現也能夠方便一些。
6.2 尺度變化連續性
通過理解這張圖，能夠把問題理解的比較清楚。
假設s=3，也就是每個塔里有3層，則k=21/s=21/3，
那么按照上圖可得Gauss Space和DoG space 分別有3個（s個）和2個（s-1個）分量，在DoG space中，1st-octave兩項分別是
σ,k
σ; 2nd-octave
兩項分別是2
σ,2k
σ;由於無法比較極值，我們必須在高斯空間繼續添加高斯模糊項，使得形成σ,kσ,k2σ,k3σ,k4σ這樣就可以選擇DoG space中的中間三項kσ,k2σ,k3σ（只有左右都有才能有極值），那么下一octave中（由上一層降采樣獲得）所得三項即為2kσ,2k2σ,2k3σ，其首項2kσ=24/3。剛好與上一octave末項k3σ=23/3尺度變化連續起來，所以每次要在Gaussian space添加3項，每組（塔）共S+3層圖像，相應的DoG金字塔有S+2層圖像。雖然我也不能給理解在不同分辨率下保持尺度的連續有何意義，但是至少具備了某種程度的完備性吧。
”尺度越大，圖像越模糊“，對於金字塔圖像來說，分辨率越小，那么包含的信息越少，通過模糊而初略的提取方法提取主要特征。
6.3、
尺度歸一化的高斯拉普拉斯算子
首先理解什么是“高斯拉普拉斯算子”，它的名字也很有名“LoG,Laplacian of Gaussian"
假設對圖像
,使用尺度為
的高斯平滑
然后再使用Laplace算子檢測邊緣
其中*是卷積符號，這是因為
定義高斯拉普拉斯算子
相當於對高斯模板先進行Laplace再與圖像進行平滑。重點在這里，前面是對高斯模板先進行Hessian求導，再和圖像平滑
將LoG展開，看其具體形式(就是二次偏導求和）
所以
定量計算的話，就可以直接認為LoG模板可以為
然后我們進一步理解
先來看看歸一化之后的LoG算子(該過程將高斯函數帶入化簡可得)，歸一的意思就是定下單位來的意思。
6.4Hessian矩陣的特征值和特征值所在特征向量方向上的曲率成正比
我當時最主要的跨越就是這幅圖：
那么應該能夠曲率的概念。Hessian矩陣的特征值分別對應了原始圖像的點垂直於線條方向和平行線條的方向。
Frangi的截圖
這個結果就是Hessian的結果，看左邊，當兩個特征值一個比較大，一個比較小的時候是直線（發分明亮和黑暗）
6.5計算“跡”和“行列式”就一定比計算矩陣特征值本身更為簡單嗎
我找到了1988這篇
里面關於使用矩陣的跡和行列式來簡化計算的描述是這樣的，也是只給了一個解雇。
6.6
二階泰勒級數相關知識
一階泰勒函數表述方式
這個公式很好理解呀，但是原文中是這樣表述
“
將每個關鍵點的一段函數進行泰勒二次展開:
其中x是相對於關鍵點的距離。和分別表示函數在關鍵點的取值和一階導數，二階導數。求導得到極值為，此時
，
當這個極值過小的時候，被去除。
原始論文中剔除
的點。
的點。
”那么這里面的這許多符號是從哪里來的？
實際上
對於高維函數的泰勒展開式為：
這個方程式就是比較像的了。迭代近視的方法叫做牛頓公式。雖然我們在圖像處理中需要解決的還是二維的問題，但是既然已經出現了高緯統一的表示方式，下面我們來推到其：
f(X),X為n維列向量，表示含有n個變量。
那么它每一項的，比如第k項可以認為是：
函數為二維時的展開公式：f(x)=f(xk)+f′(xk)(x−xk)+12f′′(xk)(x−xk)2
相應的高維展開式，符號定義：
f′(Xk)=gk=∇f(Xk)
f′′(Xk)=Gk=∇2f(Xk)
僅僅是定義而已，問題不大，那么
這也是很直接。
可以看出∇f(X)和X−Xk都是列向量，無法相乘取內積。所以這里要將gk轉置。后面的（X−Xk）2也是一樣的道理。
所以對於高維函數的泰勒展開式為：
所以對於高維函數的泰勒展開式為：
現在就能夠看到轉置的用途了。函數表示出來以后，為了求極值，那么就是要求導：分為3個部分
f(X)1=f(Xk)
f(X)2=f′(Xk)(X−Xk)
f(X)3=12f′′(Xk)(X−Xk)2
（1） 對
f(X)1求導，
f(Xk)不含有變量
X所以此項導數為0。
（2） 對 f(X)2求導。
gTk不含有變量 X所以為常數，不需要考慮，這時要求導的部分變成了 y=(X−Xk)。
為簡便，我設X為3維的列向量
（2） 對 f(X)2求導。
gTk不含有變量 X所以為常數，不需要考慮，這時要求導的部分變成了 y=(X−Xk)。
為簡便，我設X為3維的列向量
y對
X求導：
y為矩陣
X為列向量。屬於矩陣對列向量求導的那一類。
根據求導公式：
根據求導公式：
求導，屬於元素對列向量求導的類型。
同理可得：
也就是說：
所以 f′(X)2=f′(Xk)=gk
（3） 對f(X)3求導。
與上一步類似。
f′(X)3=f′′(Xk)=Gk
與上一步類似。
f′(X)3=f′′(Xk)=Gk
所以：
f′(X)=gk+Gk(X−Xk)=0
若Gk為非奇異（則矩陣可逆），解這個方程，記其解為Xk+1即得牛頓法的迭代公式：
Xk+1=Xk−G−1kgk
Xk+1=Xk−G−1kgk
第七章 小結和參考文獻
http://blog.csdn.net/u011310341/article/details/49496229
http://blog.csdn.net/abcjennifer/article/details/7372880
http://blog.csdn.net/abcjennifer/article/details/7365882
http://en.wikipedia.org/wiki/Scale-invariant_feature_transform#David_Lowe.27s_method
http://blog.sciencenet.cn/blog-613779-475881.html
http://www.cnblogs.com/linyunzju/archive/2011/06/14/2080950.html
http://www.cnblogs.com/linyunzju/archive/2011/06/14/2080951.html
http://blog.csdn.net/ijuliet/article/details/4640624
http://www.cnblogs.com/cfantaisie/archive/2011/06/14/2080917.htmS