FROM：?http://www.cnblogs.com/ronny/p/4028776.html

1. SIFT算法中一些符號的說明

I(x,y)表示原圖像。

G(x,y,σ)表示高斯濾波器，其中G(x,y,σ)=12πσ2exp(?(x2+y2)/2σ2)。

L(x,y,σ)表示由一個高斯濾波器與原圖像卷積而生成的圖像，即L(x,y,σ)=G(x,y,σ)?I(x,y)。一系列的σi，則可以生成一系列的L(x,y,σi)圖像，此時我們把這一系列的L(x,y,σ)圖像稱為原圖像的一個尺度空間表示。關于尺度空間的知識可以參考：圖像特征提取：尺度空間理論。

DOG表示高斯差分（Difference of Gaussians），也可以表示為D(x,y,σ)，其中D(x,y,σ)=(G(x,y,kσ)–G(x,y,σ))?I(x,y)=L(x,y,kσ)–L(x,y,σ)。

上面特別值得注意的是尺度為σ的高斯差分圖像由于尺度為kσ與尺度為σ的L圖像生成的。k為兩相鄰尺度空間倍數的常數。

O：高斯金字塔的組數（Octave），其中值得注意的是在實際構建中，第一組的索引可以為0也可以為-1，這個在后面解釋原理。

S：高斯金字塔每一組的層數。在實際最開始構建尺度空間圖像，即L圖像的時候，構建了S+3層，一定要把這個S+3與S區分開，為什么是S+3后面分析。

2. 構建高斯差分金字塔

2.1 第一組第一層圖像的生成

很多初涉SIFT的都會被這個問題所困惑，這里要分兩種情況：其一是把第一組的索引定為0；其二是把第一組的索引定為-1。

我們先考慮第一組索引為0的情況，我們知道第一組第一層的圖像是由原圖像與σo（一般設置為1.6）的高斯濾波器卷積生成，那么原圖像是誰呢？是I(x,y)嗎？不是！為了圖像反走樣的需要，通常假設輸入圖像是經過高斯平滑處理的，其值為σn=0.5，即半個像元。意思就是說我們采集到的圖像I(x,y)，已經被σ=σn=0.5的高斯濾波器平滑過了。所以我們不能直接對I(x,y)用σ0的高斯濾波器平滑，而應該用σ=σ20?σ2n???????√的高斯濾波器去平滑I(x,y)，即

FirstLayer(x,y)=I(x,y)?G(x,y,σ20?σ2n???????√)

其中FirstLayer(x,y)表示整個尺度空間為第1組第1層的圖像，σo一般取1.6，σn=0.5。

現在我們來考慮把第一組的索引定為-1的情況。那么首先第一個問題便是為什么要把索引定為-1。如果索引為0，如上面那種情況所示，整個尺度空間的第1組的第1層圖像已經是由原圖像模糊生成的了，那么也就是說已經丟失了細節信息，那么原圖像我們完全沒有利用上。基于這種考慮，我們先將圖像放大2倍，這樣原圖像的細節就隱藏在了其中。由上面一種情況分析，我們已經知識了I(x,y)看成是已經被σn=0.5模糊過的圖像，那么將I(x,y)放大2倍后得到Is(x,y)，則可以看為是被2σn=1的高斯核模糊過的圖像。那么由Is生成第1組第1層的圖像用的高斯濾波器的σ=σ20?(2σn)2??????????√。可以表示為。

FirstLayer(x,y)=Is(x,y)?G(x,y,σ20?(2σn)2??????????√)

其中FirstLayer(x,y)表示整個尺度空間為第1組第1層的圖像，Is(x,y)是由I(x,y)用雙線性插值放大后的圖像。σo一般取1.6，σn=0.5。

2.2 尺度空間生成了多少幅圖像

我們知道S是我們最終構建出來的用來尋找特征點的高斯差分圖像，而特征點的尋找需要查找的是空間局部極小值，即在某一層上查找局部極值點的時候需要用到上一層與下一層的高斯差分圖像，所以如果我們需要查找S層的特征點，需要S+2層高斯差分圖像，然后查找其中的第2層到第S+1層。

而每一個高斯差分圖像G(x,y,σ)都需要兩幅尺度空間的圖像L(x,y,kσ)與L(x,y,σ)進行差分生成，這里假設S =3，則我們需要的高斯差分圖像有S+2 = 5張，分別為G(x,y,σ),G(x,y,kσ),G(x,y,k2σ),G(x,y,k3σ),G(x,y,k4σ)。其中的G(x,y,kσ),G(x,y,k2σ),G(x,y,k3σ)這三張圖像是我們用來查找局部極值點的圖像。那么我們則需要S+3 = 6張尺度空間圖像來生成上面那些高斯差分圖像，它們分別為：L(x,y,σ),L(x,y,kσ),L(x,y,k2σ),L(x,y,k3σ),L(x,y,k4σ)，L(x,y,k5σ)

從上面的分析,我們知道對于尺度空間來說,我們一共需要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層的圖像，它的尺度為σ=σoko+s/S，其中,k=1/2,o∈[0,1,2,…,nOctave?1],s∈[0,1,2,…,S+2]。那么我們從第0組開始，看它各層的尺度。

第0組：σo→21/3σ0→22/3σ0→23/3σ0→24/3σ0→25/3σ0

第1組：2σo→2?21/3σ0→2?22/3σ0→2?23/3σ0→2?24/3σ0→2?25/3σ0

我們只分析2組便可以看出，第1組的第0層圖像恰好與第0組的倒數第三幅圖像一致，尺度都為2σ0，所以我們不需要再根據原圖來重新卷積生成每組的第0張圖像，只需采用上一層的倒數第3張來降采樣即可。

我們也可以繼續分析，第0組尺度空間得到的高斯差分圖像的尺度為：σo→21/3σ0→22/3σ0→23/3σ0→24/3σ0

而第1組尺度空間得到的高斯差分圖像的尺度為：2σo→2?21/3σ0→2?22/3σ0→2?23/3σ0→2?24/3σ0

如果我們把它們的中間三項取出來拼在一起，則尺度為：21/3σ0→22/3σ0→23/3σ0→2?21/3σ0→2?22/3σ0→2?23/3σ0，正好連續！！這一效果帶來的直接的好處是在尺度空間的極值點確定過程中，我們不會漏掉任何一個尺度上的極值點，而是能夠綜合考慮量化的尺度因子。

2.4 用第i-1層的圖像生成第i層的圖像

值得注意的是，在SITF的源碼里，尺度空間里的每一層的圖像（除了第1層）都是由其前面一層的圖像和一個相對sigma的高斯濾波器卷積生成，而不是由原圖和對應尺度的高斯濾波器生成的，這一方面是因為我前面提到的不存在所謂意思上的“原圖”，我們的輸入圖像I(x,y)已經是尺度為σ=0.5的圖像了。另一方面是由于如果用原圖計算，那么相鄰兩層之間相差的尺度實際上非常小，這樣會造成在做高斯差分圖像的時候，大部分值都趨近于0，以致于后面我們很難檢測到特征點。

基于上面兩點原因（個人認為原因1是最主要的，原因2只是根據實際嘗試后的一個猜想，并無理論依據），所以對于每一組的第i+1層的圖像，都是由第i層的圖像和一個相對尺度的高斯濾波器卷積生成。

那么相對尺度如何計算呢，我們首先考慮第0組，它們的第i+1層圖像與第i層圖像之間的相對尺度為SigmaDiffi=(σ0ki+1)2–(σ0ki)2???????????????√，為了保持尺度的連續性，后面的每一組都用這樣一樣相對尺度（SIFT實際代碼里是這樣做的）。這里有一個猜測，比如說尺度為2σ0的這一組，第i層與第i+1層之間的相對尺度計算的結果應該為(2σ0ki+1)2–(2σ0ki)2????????????????√=2?SigmaDiffi，可是代碼里依然用SigmaDiffi是因為這一層被降維了。

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×2個點共26個點比較，以確保在尺度空間和二維圖像位置空間都檢測到極值點。也就是，比較是在一個3×3的立方體內進行。

搜索過程從每組的第二層開始，以第二層為當前層，對第二層的DoG圖像中的每個點取一個3×3的立方體，立方體上下層為第一層與第三層。這樣，搜索得到的極值點既有位置坐標（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)≈f(0)+f′(0)x+f′′(0)2x2

另外我們知道f(x)在x處的導數寫成離散的形式為f′(x)=f(x+1)–f(x)2，二階導數寫成離散形式為f′′(x)=f(x+1)+f(x?1)?2f(x)。

所以，我們可以算出f(x)≈6+2x+?62x2=6+2x?3x2

求取函數f(x)的極大值和極大值所在的位置：

f′(x)=2?6x=0,???x^=13

f(x^)=6+2×13–3×(13)2=613

現在回到我們SIFT點檢測中來，我們要考慮的是一個三維問題，假設我們在尺度為σ的尺度圖像D(x,y)上檢測到了一個局部極值點，空間位置為(x,y,σ)，由上面的分析我們知道，它只是一個離散情況下的極值點，連續情況下，極值點可能落在了(x,y,σ)的附近，設其偏離了(x,y,σ)的坐標為(Δx,Δy,Δσ)。則對D(Δx,Δy,Δσ)可以表示為在點(x,y,σ)處的泰勒展開：

D(Δx,Δy,Δσ)=D(x,y,σ)+[?Dx?Dy?Dσ]???ΔxΔyΔσ???+12[ΔxΔyΔσ]????????????2D?x2?2D?y?x?2D?σ?x?2D?x?y?2D?y2?2D?σ?y?2D?x?σ?2D?y?σ?2D?σ2??????????????ΔxΔyΔσ???

可以將上式寫成矢量形式如下：

D(x)=D+?DT?xΔx+12ΔxT?2DT?x2Δx

令上式的一階導數等于0，可以求得Δx=??2D?1?x2?D(x)?x

通過多次迭代(Lowe算法里最多迭代5次)，得到最終候選點的精確位置與尺度x^，將其代入公式求得D(x^)，求其絕對值得|D(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×2的Hessian矩陣H求出：

H(x,y)=[Dxx(x,y)Dxy(x,y)Dxy(x,y)Dyy(x,y)]

上式中，D值可以通過求取鄰近點像素的差分得到。H的特征值與D的主曲率成正比例。我們可以避免求取具體的特征值，因為我們只關心特征值的比例。令α=λmax為最大的特征值，β=λmin為最小的特征值，那么，我們通過H矩陣直跡計算它們的和，通過H矩陣的行列式計算它們的乘積：

Tr(H)=Dxx+Dyy=α+β

Det(H)=DxxDyy?(Dxy)2=αβ

如果γ為最大特征值與最小特征值之間的比例，那么α=γβ，這樣便有

Tr(H)2Det(H)=(α+β)2αβ=(γ+1)2γ

上式的結果只與兩個特征值的比例有關，而與具體特征值無關。當兩個特征值相等時，(γ+1)2γ的值最小，隨著γ的增加，(γ+1)2γ的值也增加。所以要想檢查主曲率的比例小于某一閾值γ，只要檢查下式是否成立：

Tr(H)2Det(H)<(γ+1)2γ

Lowe在論文中給出的γ=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;

總結

以上是生活随笔為你收集整理的SIFT定位算法关键步骤的说明的全部內容，希望文章能夠幫你解決所遇到的問題。

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