個人博客原文：外鏈網址已屏蔽

這是一次作業，內容是給出兩張圖像，檢測特征點和匹配特征點。要求不能用諸如OpenCV的現成特征點檢測函數。于是就只能造輪子了，寫了這個Matlab版的sift。(其實就是把c語言的opensift翻譯成了matlab

以下是算法流程，其實網上的類似博文已經很多了，只不過我看的過程中也看得不很明白，只能對照著好幾個看，所以干脆自己又寫了一遍。下面的圖均來自于參考資料中，然而參考資料的圖也是來自于參考資料的參考資料中。

1. 構建尺度空間

定理1 對圖像做 σ=σ1 的高斯平滑，再做一次 σ=σ2 的高斯平滑，等效于對原圖像做一次 σ=σ21+σ22??????√ 的高斯平滑。

1.1 構建高斯金字塔

高斯卷積核是實現尺度變換的唯一線性核(Koenderink, 1984; Lindeberg, 1994)。

一幅圖像的尺度空間被定義為對其做可變尺度的高斯卷積：

L(x,y,σ)其中G(x,y,σ)=G(x,y,σ)?I(x,y)=12πσ2e?(x2+y2)/2σ2

對于給定的彩色圖像，轉化為灰度圖像，用不同大小的σ做高斯平滑(按照 3σ 準則，高斯核矩陣的大小設為 (6σ+1)?(6σ+1) ，并保證行和列為奇數)，再此基礎上將圖像降采樣得到不同大小的組(octave)，每組若干圖像(interval)。詳細描述如下：

為了得到更多的特征點，將圖像擴大為原來的兩倍。假設原圖像已有 σ=0.5 的高斯平滑，而我們需要第一個octave的第一張圖像的 σ=1.6 ，按照定理1，我們要對擴大兩倍的圖像做一次高斯平滑，σ=1.62?(0.5×2)2?????????????√ 。

上一個octave的圖像的長度和寬度分別是下一個octave的圖像的兩倍。因此圖像組數(octaves)可由圖像大小決定，將其設為 log2(min(height,width)) ?2 ，這樣將使頂層octave圖像的長度和寬度最小值在8像素左右。

設第m個octave的第n張圖像相對于原始圖像的參數σ為 sigma(m,n)，則sigma(1,1)=σ0=1.6。每個octave有s+1張圖像(即intervals)，這樣得到的高斯差分金字塔(DoG)每個octave將有s張圖像，我們設s為3。為了滿足在不同octave間尺度的連續性，并使 sigma(m,n)= 2?sigma(m?1,n)，按照定理1，則：

sigma(1,n)sigma(m,n)=σ0?kn?1,其中k=21/s=σ0?2m?1?kn?1

如上圖所示，在第一個octave中尺度為k3?σ0的“最后”一張圖像進行下采樣得到第二個octave的第一張圖像，尺度仍為k3?σ0=2?σ0。

但實際上我們需要做出更多不同尺度的高斯平滑圖像，這是因為在后續高斯差分金字塔的極值檢測中，需要前后兩級尺度都存在圖像。如圖中紅框所示，高斯差分金字塔中每個octave有s幅圖像，則需要高斯金字塔中每個octave包含s+3幅圖像。其中第s+1幅圖像用作下一個octave第一幅圖像的降采樣。

具體實現中并未對單幅圖像多次進行高斯平滑，而是由上一幅圖像進行高斯平滑得到下一幅圖像并迭代之，按照定理1計算σ。

1.2 構建高斯差分金字塔

對兩幅高斯金字塔的圖像作差。

1.3 檢測極值點

如上圖，與前后兩幅圖像及自身的共26個鄰域像素點比較灰度值檢測極值。

2. 關鍵點精確定位

檢測到的極值點是離散的，通過三元二次函數擬合來精確確定關鍵點的位置和尺度，達到亞像素精度。以某關鍵點為中心的尺度空間函數 D(x,y,intvl) 的二次泰勒展開式為：

D(X)=D+?D?XTX+12XT?2D?X2X

其中等號右邊第一個D為某關鍵點處的灰度值， X=(x,y,intvl)T 是以此點為中心的偏移量，由于 D(X) 是離散的，其導數用差分法求得。令 D(X) 導數為零，得到精確極值位置的偏移量為：

X^=??2D?X2?1?D?X

若X^在任意一個維度大于0.5，說明極值點精確位置距離另一個點更近，應該將關鍵點定位于更近的那個位置。定位到新點后再進行相同操作，若迭代5次位置仍不收斂，則不認為此點為關鍵點。設定圖像邊緣img_border，若關鍵點落在圖像邊緣區域(以img_border為寬度的矩形外框)也不認為此點為關鍵點。

2.1 去除低反差(low contrast)的點

精確極值點處函數值：

D(X^)=D+12?D?XTX^

若 |D(X^)|<0.04/s ，同樣不認為此點是極值點。在此過程中保存極值點的數據ddata，為特征的構建做準備。計算出σ_octv，即位于一個相同的octave內的尺度，某個octave內第n張圖像的 σ_octv=σ0?kintvl?1 ，此處intvl為精確定位后的intvl。

2.2 消除邊緣響應

高斯差分函數有較強的邊緣響應，對于比較像邊緣的點應該去除掉。這樣的點的特征為在某個方向有較大主曲率，而在垂直的方向主曲率很小。

設r為大主曲率與小主曲率的比值，H為關鍵點處的Hessian矩陣，則有(具體推導可見Lowe的論文)：

Tr(H)2Det(H)=(r+1)2r

若滿足：

Tr(H)2Det(H)

說明此處r較小，認為此關鍵點不位于邊緣，否則則去除該點。

3. 方向指定

根據關鍵點的局部特性為每個關鍵點指定一個方向，可以具備旋轉不變性。關鍵點局部特性在檢測到關鍵點的高斯差分金字塔圖像臨近的高斯金字塔圖像中計算。在關鍵點3σ鄰域窗口計算梯度和方向分布，計算方式如下：

m(x,y)θ(x,y)=[L(x+1,y)?L(x?1,y)]2+[L(x,y+1)?L(x,y?1)]2????????????????????????????????????????????√=tan?1{[L(x,y+1)?L(x,y?1)]/[L(x+1,y)?L(x?1,y)]}

此處的x正方向向右，y正方向向上。其中L為關鍵點在上述精確定位后所處尺度的灰度值，m(x,y)為梯度的幅值，θ(x,y)為關鍵點處梯度方向的弧度(范圍為(?π,π])。將360度的方向劃分為36個區域(bins),第一個區域的范圍是[35π36,37π36),按逆時針方向依次劃分。對m(x,y)按 σ=1.5σ_octv 的高斯分布，在 3σ=3?1.5σ_octv 的鄰域窗口加權計算，得到36個方向的直方圖。然后對直方圖進行兩次平滑處理，即按0.25,0.5,0.25的大小對每3個連續的bin加權兩次。

直方圖最大值的方向代表該關鍵點的主方向，對于其他峰值，若大于或等于主方向值的80%，則再分配一個方向。所以對于一個關鍵點，可能會有多個對應的方向，將帶有方向的關鍵點定義為feature，則一個關鍵點可能對應多個feature。由于第一個octave是雙倍大小的圖像，feature的坐標和尺度應轉換到原始圖像所在的octave處理。最后用拋物插值精確定位feature的方向。

對于x為-1,0,1，y為l,c,r的三個點來說，拋物插值得到極值點的x為：

0.5?l?rl?2c+r

4. 關鍵點描述子

上一步已得到具有主方向的關鍵點，即feature，下一步是對feature的鄰域進行采樣，形成對該局部圖像的描述，然后可用某種度量方法對描述進行匹配。

Lowe提出的sift描述子是一個 4×4×8=128 維的向量。描述子的數學形式可定義為 h(x,y,θ) ，其中的x,y代表 4×4=16 個圖像區域的位置，θ即梯度方向，只能取8個值。h(x,y,θ)的值就是在(x,y)代表的圖像區域計算得到的在θ方向的梯度大小。

4.1 描述子采樣區域

這16個圖像區域中的每一個區域均為 3σ_octv 像素，因此16個區域的半邊長為 4×3σ_octv/2 ，考慮到后續操作需要三線性插值，采樣區域半邊長設為 (4+1)×3σ_octv/2 ，又由于旋轉操作，這個值需要乘以\sqrt{2}2√，得到 radius=(4+1)× 2√×3σ_octv/2 。

如下圖所示，圖中 m=3, Bp=4,σ=σ_octv 。

4.2 旋轉至主方向

為了使描述符具有旋轉不變性，將坐標軸旋轉至關鍵點主方向。設i，j分別為采樣點相對關鍵點的行偏移量和列偏移量，i = -radius:radius，j = -radius:radius，關鍵點左上角i和j均為負數。關鍵點主方向為θ，范圍是(?π,π]。

定理2 在右手平面直角坐標系中，向量(x,y)逆時針旋轉θ，得到的向量(x’,y’)為：

[x′y′]=[cosθsinθ?sinθcosθ][xy]

在左圖以關鍵點為中心建立右手平面直角坐標系o-uv，u的正方向與左圖x方向相同，v的正方向與左圖y方向相反。左圖中x=(1,0),y=(0,-1)，將x,y代入定理2的公式，得到右圖中 x′=(cosθ,sinθ), y′=(sinθ,?cosθ) 。其中θ為左圖坐標系旋轉到右圖坐標系的角度，在上圖中為一負數。設圖像中有一點按左圖的x,y可表示為 j?x+i?y ，按右圖中的x′,y′可表示為 j′?x′+i′?y′ ，則有：

j[10]+i[0?1]=j′[cosθsinθ]+i′[sinθ?cosθ]

解得 j′=j?cosθ?i?sinθ, i′=j?sinθ+i?cosθ 。

得到新的行、列偏移量后，將 3σ_octv 設為單位長度，并將中心轉移至最左上角區域中心，得到新的坐標r_bin和c_bin。對梯度方向弧度值減去主方向弧度，并設 2π8 為一個單位，得到o_bin。

采樣點的梯度幅值按照 σ=0.5?4?3σ_octv (即16個區域邊長的一半)的高斯函數加權：

w=m(a+j,b+i)?e?j′2+i′22σ2

其中a，b為關鍵點在高斯金字塔圖像中的位置坐標。

4.3 三線性插值

上述過程中構造了一個三維的bin空間，如4.1中右圖所示，維度包括r_bin，c_bin和o_bin。注意最上層格子和最底層格子是相連的，因為0度等于360度。所有帶有三維坐標的梯度幅值都將分配到三維格子里。

為了減少一個梯度幅值從一個格子漂移(shift)到另一個格子引起的描述子突變，需要對梯度值做三線性插值。也就是根據三維坐標計算距離周圍格子的距離，按距離的倒數計算權重，將梯度幅值按權重分配到臨近的格子里。

某點在三維bin空間的坐標為(r_bin，c_bin，o_bin)，求出r=?r_bin?, c=?c_bin?, o=?o_bin?, dr=r_bin?r, dc=c_bin?c, do=o_bin?o，它的梯度幅值最多可能分配到周圍的8個格子中。計算公式如下：

=weightedValue[r+i+1][c+j+1][((o+k)mod8)+1]w?dri?(1?dr)1?i?dcj?(1?dc)1?j?dok?(1?do)1?k

其中i，j，k均可取0或1，weightedValue下標加1的目的是使下標從1開始。

為簡化計算，可改為：

=weightedValue[r+i+1][c+j+1][((o+k)mod8)+1]w?(0.5+(dr?0.5)(2i?1))?(0.5+(dc?0.5)(2j?1))?(0.5+(do?0.5)(2k?1))

4.4 生成描述子

將上述直方圖數組按順序排列可轉換為一個128維的向量。

為了減少光照變化的影響，對該向量進行歸一化處理。非線性光照變化仍可能導致梯度幅值的較大變化，然而影響梯度方向的可能性較小。因此對于超過閾值0.2的梯度幅值設為0.2，然后再進行一次歸一化。最后將描述子按照對應高斯金字塔圖像的尺度大小排序。

5. 匹配

描述子向量已經歸一化，所以可直接用向量之間的夾角進行匹配，相當于球面距離。圖像A 的描述子匹配圖像B最近的兩個描述子點積之比小于0.6，則認為匹配成功。

6. 一些廢話

6.1 性能優化

因為用的是Matlab，所以不注重性能。然而又不得不注重性能，因為第一次跑通程序的時候跑了一晚上都沒跑完一半！也就是一張圖片的描述子都沒算完。后來發現是因為在運行次數最多的for循環(描述子計算中的梯度計算)里用到了cell數組。把對這個cell數組的查詢操作提到兩重循環前以后，這個程序好像跑了半個小時左右跑出結果了。然而還是太慢，于是我又用Matlab的計時分析工具分析了程序最耗時的地方:

把cell數組的查詢盡可能減少

充分利用Matlab的向量操作

一些沒用的參數給去掉了(如計算梯度時的三個返回值合并到了一個二維數組)

一個三維數組折疊成了一維的(hist)

程序里用了很多全局變量，是因為我把函數分成了文件而不是放在一個文件，為了節省點內存(以及方便)只能這么做(雖然據說Matlab在不改變變量的情況下函數傳值等于引用，然而我并不清楚究竟是怎樣的)。把for循環換成parfor的時候提示，parfor里似乎不推薦用全局變量，而且實際運行的時候全局變量似乎也會影響性能，于是我把全局變量復制成了局部的再放進parfor里。

我還發現一個奇葩的問題，在運行次數最多的計算梯度的函數里用zeros(1,2)創建一個數組竟然也耗時非常多，改成[0 0]就好了。

經過這些修改后，在開啟parallel pool的情況下運行時間縮短到了7分鐘左右。(然而Lowe的C語言版本只要十幾秒

6.2 運行結果

這次作業老師給的是兩張768x1024的圖片，分別檢測到5288和4798個特征點，最后匹配了906對點。用Lowe的siftDemoV4跑出來的結果是1252對匹配。

這個程序的參數基本都是參照opensift，但最后的匹配用的是Lowe的方案。Lowe的實現畢竟不太一樣，運行的結果和opensift有一些差異。以下是匹配siftDemoV4.zip里的scene.pgm和book.pgm的結果：

Item

siftDemoV4

opensift

sw-sift

scene.pgm

1021

746

766

book.pgm

882

740

741

Matched

98

84

58

sw-sift和opensift的區別主要是在高斯平滑和匹配算法上。opensift的高斯平滑用的是OpenCV的CVSmooth函數，匹配用的是歐式距離(而且把描述子乘以512從double類型轉成了int)。和opensift相比，sw-sift檢測到的特征點數量很接近，但是匹配數量較少，所以可改進的地方主要是匹配算法(然而我不想改了==)。另外，我發現高斯平滑的核矩陣大小對結果有很大影響，根據3σ準則它的寬度應該是 (6σ+1)?(6σ+1) ，然而有人設成 (3σ+1)?(3σ+1) 卻取得了更多特征點，因此調整這個參數再用其它參數限制錯誤數量或許可以得到更好的結果。

6.3 源代碼

代碼發布在github: sw-sift，請注意sift是有專利的。

參考資料

David G. Lowe, “Distinctive image features from scale-invariant keypoints,” International Journal of Computer Vision, 60, 2 (2004), pp. 91-110. [PDF] [CODE]

Rob Hess, OpenSIFT源碼: 外鏈網址已屏蔽

zddhub, SIFT算法詳解: 外鏈網址已屏蔽

Rachel Zhang, SIFT特征提取分析: 外鏈網址已屏蔽

JiePro, SIFT算法：特征描述子: 外鏈網址已屏蔽

總結

以上是生活随笔為你收集整理的matlab实现sift,SIFT算法的Matlab实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。