sift算法解析

　　SIFT綜述

　　尺度不變特征轉(zhuǎn)換（Scale-invariant feature transform或SIFT）是一種電腦視覺的算法用來偵測與描述影像中的局部性特征，它在空間尺度中尋找極值點，并提取出其位置、尺度、旋轉(zhuǎn)不變量，此算法由 David Lowe在1999年所發(fā)表，2004年完善總結(jié)。

　　其應(yīng)用范圍包含物體辨識、機器人地圖感知與導(dǎo)航、影像縫合、3D模型建立、手勢辨識、影像追蹤和動作比對。

　　此算法有其專利，專利擁有者為英屬哥倫比亞大學(xué)。

　　局部影像特征的描述與偵測可以幫助辨識物體，SIFT 特征是基于物體上的一些局部外觀的興趣點而與影像的大小和旋轉(zhuǎn)無關(guān)。對于光線、噪聲、些微視角改變的容忍度也相當(dāng)高?；谶@些特性，它們是高度顯著而且相對容易擷取，在母數(shù)龐大的特征數(shù)據(jù)庫中，很容易辨識物體而且鮮有誤認(rèn)。使用 SIFT特征描述對于部分物體遮蔽的偵測率也相當(dāng)高，甚至只需要3個以上的SIFT物體特征就足以計算出位置與方位。在現(xiàn)今的電腦硬件速度下和小型的特征數(shù)據(jù)庫條件下，辨識速度可接近即時運算。SIFT特征的信息量大，適合在海量數(shù)據(jù)庫中快速準(zhǔn)確匹配。

　　SIFT算法的特點有：

　　1、SIFT特征是圖像的局部特征，其對旋轉(zhuǎn)、尺度縮放、亮度變化保持不變性，對視角變化、仿射變換、噪聲也保持一定程度的穩(wěn)定性；

　　2、獨特性（Distinctiveness）好，信息量豐富，適用于在海量特征數(shù)據(jù)庫中進行快速、準(zhǔn)確的匹配；

　　3、 多量性，即使少數(shù)的幾個物體也可以產(chǎn)生大量的SIFT特征向量；

　　4、 高速性，經(jīng)優(yōu)化的SIFT匹配算法甚至可以達到實時的要求；

　　5、可擴展性，可以很方便的與其他形式的特征向量進行聯(lián)合。

　　SIFT算法可以解決的問題：

　　目標(biāo)的自身狀態(tài)、場景所處的環(huán)境和成像器材的成像特性等因素影響圖像配準(zhǔn)/目標(biāo)識別跟蹤的性能。而SIFT算法在一定程度上可解決：

　　1、目標(biāo)的旋轉(zhuǎn)、縮放、平移（RST）

　　2、圖像仿射/投影變換（視點viewpoint）

　　3、光照影響（illumination）

　　4、目標(biāo)遮擋（occlusion）

　　5、雜物場景（clutter）

　　6、噪聲

　　SIFT算法的實質(zhì)是在不同的尺度空間上查找關(guān)鍵點（特征點），并計算出關(guān)鍵點的方向。SIFT所查找到的關(guān)鍵點是一些十分突出，不會因光照，仿射變換和噪音等因素而變化的點，如角點、邊緣點、暗區(qū)的亮點及亮區(qū)的暗點等。

　　SIFT的缺點

　　SIFT在圖像的不變特征提取方面擁有無與倫比的優(yōu)勢，但并不完美，仍然存在：

　　1、實時性不高。

　　2、有時特征點較少。

　　3、對邊緣光滑的目標(biāo)無法準(zhǔn)確提取特征點。

　　sift算法matlab代碼的實現(xiàn)

　　1、構(gòu)建尺度空間

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

　　1.1、構(gòu)建高斯金字塔

　　高斯卷積核是實現(xiàn)尺度變換的唯一線性核（Koenderink， 1984; Lindeberg， 1994）。

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

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

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

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

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

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

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

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

　　1.2 構(gòu)建高斯差分金字塔

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

　　1.3 檢測極值點

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

　　2、關(guān)鍵點精確定位

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

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

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

　　2.1 去除低反差（low contrast）的點精確極值點處函數(shù)值：

　　同樣不認(rèn)為此點是極值點。在此過程中保存極值點的數(shù)據(jù)ddata，為特征的構(gòu)建做準(zhǔn)備。計算出σ_octv，即位于一個相同的octave內(nèi)的尺度，某個octave內(nèi)第n張圖像的 σ_octv=σ0?kintvl?1 ，此處intvl為精確定位后的intvl。

　　2.2 消除邊緣響應(yīng)

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

　　設(shè)r為大主曲率與小主曲率的比值，H為關(guān)鍵點處的Hessian矩陣，則有（具體推導(dǎo)可見Lowe的論文）：

　　其中rt為一閾值，設(shè)為10說明此處r較小，認(rèn)為此關(guān)鍵點不位于邊緣，否則則去除該點。

　　3、方向指定

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

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

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

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

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

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

　　4.1、描述子采樣區(qū)域

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

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

　　4.2、旋轉(zhuǎn)至主方向

　　為了使描述符具有旋轉(zhuǎn)不變性，將坐標(biāo)軸旋轉(zhuǎn)至關(guān)鍵點主方向。設(shè)i，j分別為采樣點相對關(guān)鍵點的行偏移量和列偏移量，i = -radius:radius，j = -radius:radius，關(guān)鍵點左上角i和j均為負數(shù)。關(guān)鍵點主方向為θ，范圍是（?π，π］。

　　定理2 在右手平面直角坐標(biāo)系中，向量（x，y）逆時針旋轉(zhuǎn)θ，得到的向量（x’，y’）為：

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

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

　　采樣點的梯度幅值按照 σ=0.5?4?3σ_octv （即16個區(qū)域邊長的一半）的高斯函數(shù)加權(quán)：

　　其中a，b為關(guān)鍵點在高斯金字塔圖像中的位置坐標(biāo)。

　　4.3、三線性插值

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

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

　　某點在三維bin空間的坐標(biāo)為（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個格子中。計算公式如下：

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

　　為簡化計算，可改為：

　　4.4、生成描述子

　　將上述直方圖數(shù)組按順序排列可轉(zhuǎn)換為一個128維的向量。

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

　　5、匹配

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

? ? ? 6、總結(jié)

　　6.1、性能優(yōu)化

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

　　1）把cell數(shù)組的查詢盡可能減少

　　2）充分利用Matlab的向量操作

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

　　4）一個三維數(shù)組折疊成了一維的（hist）

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

　　我還發(fā)現(xiàn)一個奇葩的問題，在運行次數(shù)最多的計算梯度的函數(shù)里用zeros（1，2）創(chuàng)建一個數(shù)組竟然也耗時非常多，改成［0 0］就好了。

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

　　6.2 運行結(jié)果

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

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

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