【計算機視覺】2. 特徵點檢測：Harris, SIFT, SURF, ORB

作者：由 Encoder 發表于攝影時間：2018-05-03

特徵點檢測：Harris， SIFT， SURF & ORB

Harris角點檢測

Def。

［角點（corner point）］

在鄰域內的各個方向上灰度變化值足夠高的點，是影象邊緣曲線上曲率極大值的點。

［基於灰度影象的角點檢測］包括基於梯度的方法（透過計算邊緣的曲率判斷角點），基於模板的方法（考慮畫素鄰域點的灰度變化，將與鄰點亮度對比足夠大的點定義為角點），基於模板梯度組合的方法

［基於二值影象的角點檢測］將二值影象作為單獨的檢測目標，可使用各種基於灰度影象的角點檢測方法

［基於輪廓曲線的角點檢測］透過角點強度或曲線曲率提取角點

角點檢測的基本思想是：使用角點檢測運算元，對影象的每個畫素計算角點響應函式（Corner Response Function ），閾值化角點響應函式，根據實際情況選擇閾值，對閾值化的角點響應函式進行非極大值抑制，並獲取非零點作為角點。透過一個小的滑動視窗在鄰域檢測角點

在任意方向上移動視窗，若視窗內的灰度值都有劇烈的變化，則視窗的中心就是角點。定義角點響應函式

$E(u,v)=\sum_{(x,y)} w(x,y) [I(x+u,y+v)-I(x,y)]^2$

其中w為窗函式（window function），I為影象梯度，那麼E就表示了灰度變化的劇烈程度。1977年，Moravec最先提出瞭如下的角點檢測方法：

對於原始影象，取偏移量（Δx，Δy）為（1，0），（1，1），（0，1），（-1，1），分別計算每一畫素點（xi，yi）的灰度變化

對於每一畫素點（xi，yi），計算角點響應函式R（xi，yi）=min E

設定閾值T，將角點響應函式R（xi，yi）中低於T的值設為0

在視窗範圍內進行非極大值抑制：遍歷角點響應函式，若某個畫素的角點響應函式在視窗內不是最大，該畫素置0

選擇非零點作為角點檢測結果

Moravec角點檢測的缺點

二值的視窗函式導致角點響應函式不夠光滑

只在四個方向上計算灰度值變化，導致角點響應函式在多處都有較大響應

對於每個點只考慮E的最小值，導致演算法對邊緣有很強的反應

1988年，Harris和Plessey對Moravec的方法進行了改進，提出了經典的Harris角點檢測演算法。Harris首先將Moravec演算法中的視窗函式由階躍函式改為二維高斯函式，並透過泰勒展開考察微小移動，也就是說，如果要求E的最大值以明確角點，就可以令$u，v \rightarrow 0$，對E做泰勒展開，得

$E(u,v)=(u,v)M\begin{pmatrix} u\\ v \end{pmatrix}\\ \text{where }M=\sum_{(x,y)}w(x,y)\begin{pmatrix} I_X^2 & I_X I_Y \\ I_X I_Y & I_Y^2 \end{pmatrix}=\begin{pmatrix} \sum_W I_X^2 & \sum_W I_X I_Y \\ \sum_W I_X I_Y & \sum_W I_Y^2 \end{pmatrix}$

記

$M=\begin{pmatrix}A & B\\ B & C\end{pmatrix}$

，則上式可以寫成

的形式，這表示了一個橢圓，自相關矩陣M描述了影象區域性區域的灰度變化趨勢，可以透過橢圓的形狀來判定角點。

Proposition。

對於橢圓

，設其半長軸和半短軸分別為a，b，那麼

$\frac{1}{\sqrt{a}},\frac{1}{\sqrt{b}}$

是矩陣M的特徵值。

Proof。

橢圓上任意一點到原點（也就是橢圓的中心）的距離平方

，條件為

，那麼拉格朗日函式

$L=x^2+y^2+\lambda(Ax^2+2Bxy+Cy^2-1)$

於是

$\begin{cases} F_x$

則

$[(1-A\lambda)x^2-B\lambda xy]+[-B\lambda xy+(1-C\lambda)y^2]=0$

又可以改寫為

$\begin{pmatrix} x & y \end{pmatrix}\begin{pmatrix} 1-A\lambda & -B\lambda \\ -B\lambda & 1-C\lambda \end{pmatrix}\begin{pmatrix} x\\y \end{pmatrix}=0$

由於（x，y）不是零向量，所以

$\begin{vmatrix}1-A\lambda & -B\lambda \\ -B\lambda & 1-C\lambda\end{vmatrix}=0$

，即1/λ是矩陣

的特徵值，所以

$\begin{cases} \frac{1}{\lambda_1}+\frac{1}{\lambda_2}=A+C=\text{tr} M\\ \frac{1}{\lambda_1 \lambda_2}=AC-B^2=\det M \end{cases}$

又

$d^2=x^2+y^2=\lambda(Ax^2+2Bxy+Cy^2)=\lambda\Rightarrow d=\sqrt{\lambda}$

這樣λ1，λ2就分別對應d^2的兩個最值，也就是a^2和b^2，所以，

$\frac{1}{\sqrt{a}},\frac{1}{\sqrt{b}}$

是矩陣M的特徵值。由此我們還可以得到，橢圓面積

$S=\pi ab=\frac{\pi}{\sqrt{AC-B^2}}$

q。e。d。

得到M的特徵值有什麼用呢？若用奇異值分解的觀點看這個問題。由於Jacobian矩陣

$J=\begin{pmatrix} I_X \\ I_Y\end{pmatrix}$

，故M=JJ^T，這樣M的特徵值開根號後就是J的奇異值，因此M的特徵值就可以體現I

X和I

Y的相對大小。

現在我們需要定義角點響應函式，進一步進行區分，令（一般k=0。04 ~ 0。06）

$R=\text{det}(M)-k \text{tr}^2 (M)$

其中

$\begin{cases} \text{det}⁡(M)=\lambda_1 \lambda_2=\sum_W I_X^2 \cdot \sum_W I_Y^2 - (\sum_W I_X I_Y)^2\\ \text{tr}(M)=\lambda_1+\lambda_2 = \sum_W I_X^2+I_Y^2 \end{cases}$

則定性判斷方法為

［邊緣］

$\lambda_2\geq \lambda_1$

或

$\lambda_2\leq \lambda_1$

［角點］

$\lambda_1,\lambda_2$

都很大，E在各個方向顯著變化

［平滑區域］

$\lambda_1,\lambda_2$

都很小，$E$在各個方向基本不變

定量判斷方法為

［邊緣］ R<<0

［角點］ R>>0

［平滑區域］ |R|很小

一般增大k的值，將減小角點響應值R，降低角點檢測的靈性，減少被檢測角點的數量；減小k值，將增大角點響應值R，增加角點檢測的靈敏性，增加被檢測角點的數量。

下面我們介紹如何利用Harris角點特徵進行特徵點匹配。

對兩幅影象進行特徵匹配的過程是：

建立影象的特徵點資料庫每個特徵點的資料結構，包括：位置座標、尺度、方向、特徵向量，

為新影象的每個特徵點在資料庫中逐個匹配，根據特徵向量的歐氏距離在資料庫中尋找其最近鄰和次近鄰特徵點，若（最近鄰距離/次近鄰距離）大於某一闕值，則特徵匹配成功。

Harris 角點檢測器僅僅能夠檢測出影象中的興趣點，但是沒有給出透過比較影象間的興趣點來尋找匹配角點的方法。我們需要在每個點上加入描述子資訊，並給出一個比較這些描述子的方法。

興趣點描述子是分配給興趣點的一個向量，描述該點附近的影象的表觀資訊。描述子越好，尋找到的對應點越好。我們用對應點或者點的對應來描述相同物體和場景點在不同影象上形成的畫素點。Harris 角點的描述子通常是由周圍影象畫素塊的灰度值，以及用於比較的歸一化互相關矩陣構成的。影象的畫素塊由以該畫素點為中心的周圍矩形部分影象構成。通常，兩個（相同大小）畫素塊 I1（x）和 I2（x）的歸一化互相關矩陣定義為：

$ncc(I_1,I_2 )=\frac{1}{n-1} \sum_x \frac{[I_1 (x)-\mu_1]}{\sigma_1} \frac{[I_2 (x)-\mu_2]}{\sigma_2}$

Harris角點的優點

計算簡單

提取的點特徵均勻且合理

穩定：Harris運算元對影象旋轉、亮度變化、噪聲影響和視點變換不敏感

Harris 運算元的侷限性

對尺度很敏感，不具有尺度不變性

提取的角點精度是畫素級的

需要設計角點匹配演算法

1994年，Shi和Tomasi對Harris角點檢測進行了改進，即定義角點響應函式R=min⁡ （λ1，λ2），該改進是計算適合跟蹤的優質特徵（good features to track），使得特徵分佈更均勻，其檢測精度是亞畫素級別的。

SIFT演算法預備知識：尺度空間與影象金字塔

1999年Lowe提出了SIFT（Scale-invariant feature transform，尺度不變性特徵變換）特徵檢測演算法，並於2003年對其完善總結；2006年，Bay等人對SIFT演算法進行了改進，提升其效率，提出了SURF（Speeded Up Robust Features，加速魯棒性特徵）演算法。

SIFT特徵檢測演算法的特點：

SIFT特徵是影象的區域性特徵，對旋轉、尺度縮放、亮度變化保持不變性，對視角變化、仿射變換、噪聲也保持一定程度的穩定性

資訊量豐富，適用於在海量特徵資料庫中進行匹配

多量性，少數物體也可以產生大量SIFT特徵

高速性，經最佳化的SIFT匹配演算法甚至可以達到實時性

SIFT特徵檢測的步驟：

檢測尺度空間的極值點

精確定位特徵點（Keypoint）

設定特徵點的方向引數

生成特徵點的描述子（128維向量）

如果一張照片的畫素數目被不斷壓縮，或者觀察者距離照片越來越遠，它將逐漸變得模糊，而這個導致照片的呈現內容發生變化的連續的自變數就可以稱為尺度。對物體觀察的尺度不同，物體呈現的方式也不同。對計算機視覺而言，無法預知某種尺度的物體結構是否有意義，因此有必要將所有尺度的結構表示出來。以下關於尺度空間理論的部分，也可以參考：

http：//www。

cnblogs。com/ronny/p/388

6013。html

利用影象金字塔的方法，我們可以得到一系列大小不一的影象，由大到小，從下到上構成的塔狀模型。假設高斯金字塔的第l 層影象為G

l，則有：

$G_l(i,j)=\sum^{2}_{m=-2}\sum_{n=-2}^2\omega(m,n)G_{l-1}(2i+m,2j+n)\\ (1\le l \le N,0 \le i \le R_l,0 \le j \le C_l)$

式中，N為高斯金字塔層數；R

l和G

l分別為高斯金字塔第l層的行數和列數；ω（m，n）是一個二維可分離的5× 5視窗函式，表示式為：

$\omega=\frac{1}{256} \begin{pmatrix} 1 & 4 & 6 & 4 & 1\\ 4 & 16 & 24 & 16 & 4\\ 6 & 24 &36 & 24 & 6\\ 4 & 16 & 24 & 16 & 4\\ 1 & 4 & 6 & 4 & 1 \\ \end{pmatrix} =\frac{1}{16}\begin{pmatrix} 1 & 4 & 6 & 4 &1 \end{pmatrix}\times\frac{1}{16}\begin{pmatrix} 1 \\ 4 \\ 6 \\ 4 \\1 \end{pmatrix}$

寫成上面的形式是為了說明，

二維視窗的卷積運算元，可以寫成兩個方向上的1維卷積核（二項核）的乘積。

上面卷積形式的公式實際上完成了2個步驟：1）高斯模糊；2）降維。按上述步驟生成的G

0，G

1，。。。，G

N就構成了影象的高斯金字塔，其中$G

0$為金字塔的底層（與原影象相同），$G

N$為金字塔的頂層。可見高斯金字塔的當前層影象是對其前一層影象先進行高斯低通濾波，然後做隔行和隔列的降取樣（去除偶數行與偶數列）而生成的。

將G

l進行內插（這裡內插用的不是雙線性而是用的與降維時相同的濾波核）得到放大影象

$G_l^\ast$

，使

$G_l^\ast$

的尺寸與

$G_{l-1}$

的尺寸相同，表示為：

$G_l^{\ast}(i,j)=4\sum_{m=-2}^2\sum_{n=-2}^2\omega(m,n)G_l(\frac{i+m}{2},\frac{j+n}{2})$

其中

$0 \le l \le N,0 \le i \le R_l,0 \le j \le G_l$

，上面的係數4，是因為每次能參與加權的項，權值和為4/256，這個與我們用的ω的值有關。式中，

$G_l(\frac{i+m}{2},\frac{j+n}{2})=\begin{cases} G_l(\frac{i+m}{2},\frac{j+n}{2}), & \frac{i+m}{2},\frac{j+n}{2}\in \mathbb{Z} \\ 0, & \text{others} \end{cases}$

令

$\begin{cases} LP_l=G_l-G_{l+1}^{\ast}, & 0 \le l \le N \\ LP_N=G_N, & l=N \end{cases}$

式中，N為拉普拉斯金字塔頂層號，LP

l是拉普拉斯金字塔分解的第$l$層影象。由LP

0，LP

1，\cdots，LP

l，。。。，LP

N構成的金字塔即為拉普拉斯金字塔。它的每一層影象是高斯金字塔本層影象與其高一級的影象經內插放大後圖像的差，此過程相當於帶通濾波，因此拉普拉斯金字塔又稱為帶通金字塔分解。

影象的金字塔化能高效地（計算效率也較高）對影象進行多尺度的表達，但它缺乏堅實的理論基礎，不能分析影象中物體的各種尺度。

我們知道，訊號的尺度空間剛提出是就是透過一系列單引數、寬度遞增的高斯濾波器將原始訊號濾波得到到組低頻訊號。那麼一個很明顯的疑問是：除了高斯濾波之外，其他帶有引數t的低通濾波器是否也可以用來生成一個尺度空間。後來Koenerink、Lindeberg、Florack等人用精確的數學形式透過不同的途徑都證明了

高斯核是實現尺度變換的唯一變換核。

雖然很多研究者從可分性、旋轉不變性、因果性等特性推出高斯濾波器是建立線性尺度空間的最優濾波器。然後在數字影象處理中，需要對核函式進行取樣，離散的高斯函式並不滿足連續高斯函式的的一些優良的性質。所以後來出現了一些非線性的濾波器組來建立尺度空間，如B樣條核函式。

使用高斯濾波器對影象進行尺度空間金塔塔圖的構建，讓這個尺度空間具有下面的性質：

加權平均和有限孔徑效應

訊號在尺度t上的表達可以看成是原訊號在空間上的一系列加權平均，權重就是具有不同尺度引數的高斯核。

訊號在尺度t上的表達也對應於用一個無方向性的孔徑函式（特徵長度為

$\sigma=\sqrt{t}$

）來觀測訊號的結果。這時候訊號中特徵長度小於σ的精細結構會被抑制，理解為一維訊號上小於σ的波動會被平滑掉。

層疊平滑，也叫高斯核族的半群（Semi-Group）性質：兩個高斯核的卷積等同於另外一個不同核引數的高斯核卷積。

$g(\mu,\sigma_1)\ast g(\mu,\sigma_2)=g(\mu,\sqrt{\sigma_1^2+\sigma_2^2})$

這個性質的意思就是說不同的高斯核對影象的平滑是連續的。

區域性極值遞性

這個特徵可以從人眼的視覺原理去理解，人在看一件物體時，離得越遠，物體的細節看到的越少，細節特徵是在減少的。對生理學的研究中發現，哺乳動物的視網膜和視覺皮層的感受區域可以很好地用4階以內的高斯微分來建模。

高斯核對影象進行濾波具有壓制區域性細節的性質。

尺度伸縮不變性，對原來的訊號加一個變換函式，對變換後的訊號再進行高斯核的尺度空間生成，新的訊號的極值點等特徵是不變的。

SIFT演算法

令原影象為金字塔的第一層，每次降取樣所得到的新影象為金字塔的一層（每層一張影象），每個金字塔共n層。金字塔的層數根據影象的原始大小和塔頂影象的大小共同決定，其計算公式如下：

$n=\log_2⁡ [\min (M,N)] –t$

其中M，N為原影象的長和寬，t為塔頂影象的最小維數的對數值。接下來，我們透過尺度空間理論模擬影象資料的多尺度特徵，以高斯函式為實現尺度變換的唯一線性核，則二維影象I（x， y）的尺度空間

$L(x,y,\sigma)=I(x,y)*g(x,y,\sigma)\\ \text{where } g(x,y,\sigma)=\frac{1}{2\pi\sigma^2 } \exp⁡ \left(-\frac{x^2+y^2}{2\sigma^2 }\right)$

根據高斯函式的性質，我們知道高斯窗的寬度約為6σ，即在［-3σ，3σ］外的區間外，訊號透過高斯濾波器的響應為0，因此，高斯函式不同的標準差就可以作為不同尺度的量度，事實上對高斯函式進行傅立葉變換後有

$\mathscr{F}[g(x,k\sigma)]=k\mathscr{F}[g(x,\sigma)]\quad \mathscr{F}[g(x,y,k\sigma)]=k^2\mathscr{F}[g(x,y,\sigma)]$

即一維高斯核是線性核，二維高斯核使得

$x\leftarrow kx,y\leftarrow ky$

，所以仍是線性的。這樣就可以定義高斯差分（difference of Gaussian， DoG）運算元

$D(x,y,\sigma,k)=L(x,y,k\sigma)-L(x,y,\sigma) \approx (k-1) \nabla^2 h$

其中

$\nabla^2 h$

表示高斯拉普拉斯運算元（LoG），正則化的高斯-拉普拉斯變換

$\nabla^2_{norm}=\sigma^2\nabla^2g=\sigma^2(\frac{\partial^2g}{\partial x^2}+\frac{\partial^2g}{\partial y^2}) = -\frac{1}{2\pi\sigma^2}[1-\frac{x^2+y^2}{\sigma^2}]\cdot \exp(-\frac{x^2+y^2}{2\sigma^2})$

求

$\frac{\partial(\nabla^2_{norm})}{\partial\sigma} = 0\Rightarrow r^2-2\sigma^2=0$

影象與某一個二維函式進行卷積運算實際就是求取影象與這一函式的相似性。同理，影象與高斯拉普拉斯函式的卷積實際就是求取影象與高斯拉普拉斯函式的相似性。當影象中的斑點尺寸與高斯拉普拉斯函式的形狀趨近一致時，影象的拉普拉斯響應達到最大。

從機率的角度解釋為：假設原影象是一個與位置有關的隨機變數X的密度函式，而LOG為隨機變數Y的密度函式，則隨機變數X+Y的密度分佈函式即為兩個函式的卷積形式。如果想讓X+Y能取到最大值，則X與Y能保持步調一致最好，即X上升時，Y也上升，X最大時，Y也最大。實際上

$\nabla^2[G(x,y)*f(x,y)] = \nabla^2[G(x,y)]*f(x,y)$

每一次取樣返回的結果，就像音樂上的一個八度（octave），為了讓尺度體現其連續性，高斯金字塔在簡單降取樣的基礎上加上了高斯濾波。將影象金字塔每層的一張影象使用不同引數做高斯模糊，使得金字塔的每層含有多張高斯模糊影象，將金字塔每層多張影象合稱為八度（Octave），金字塔每層只有一組影象，組數和金字塔層數相等，每組含有多張（也叫層，Interval）影象。另外，降取樣時，高斯金字塔上一組影象的初始影象（底層影象）是由前一組影象的倒數第三張影象隔點取樣得到的。在下一個八度，高斯運算元的標準差擴大一倍。

在高斯差分尺度空間檢測區域性極大或極小值，檢測點與其同尺度的8個相鄰點、上下相鄰尺度對應的9 × 2個點進行比較，以確保在尺度空間和二維影象空間都檢測到極值點，極值點的位置可以透過對高斯差分運算元求一階導數得到。離散空間的極值點並不是真正的極值點，利用已知的離散空間點插值得到的連續空間極值點的方法叫做亞畫素插值（Sub-pixel Interpolation）。為了提高關鍵點的穩定性，需要對尺度空間DoG函式進行曲線擬合。利用DoG函式在尺度空間的Taylor展開式（擬合函式）為：

$D(X)=D+\frac{\partial D^T}{\partial X} X+\frac{1}{2} X^T \frac{\partial^2 D}{\partial X^2} X \\ \text{where } X=[x,y,\sigma,k]^T$

求導並讓方程等於零，可以得到極值點的偏移量為：

$\hat X= -\frac{\partial^2 D^{-1}}{\partial X^2 } \frac{\partial D}{\partial X}$

對於這樣得到的極值點，首先要捨去高斯差分運算元絕對值很小（一般指|D（\hat X）|<0。03）的點，透過主曲率分析去除邊緣響應過大的極值點，即計算差分影象D的Hessian矩陣：

$H= \begin{pmatrix} D_{XX} & D_{XY} \\ D_{XY} & D_{YY} \end{pmatrix}$

保留滿足如下條件的極值點（一般r=10）

$\frac{\text{tr}^2 H}{\text{det} H} < \frac{(r+1)^2}{r}$

現在對L（x，y，σ）進行差分，得到鄰域梯度方向直方圖（這裡常用Roberts運算元計算），即

$\begin{gather} m(x,y)=\sqrt{(L(x+1,y)-L(x-1,y))^2+(L(x,y+1)-L(x,y-1))^2}\\ \theta(x,y)=\arctan \frac{L(x,y+1)-L(x,y-1)}{L(x+1,y)-L(x-1,y)} \end{gather}$

進而確定特徵的方向引數：

梯度方向直方圖的範圍：［0， 360］；

特徵點的主方向：直方圖的峰值；

特徵點的輔方向：直方圖的次峰值，能量超過主峰值能量的80\%；

當影象發生旋轉時，特徵點的主方向相對於畫素的梯度方向不變；將多幅待匹配的影象都旋轉到令特徵點方向為0的位置再匹配，使特徵具有旋轉不變性。

現在我們將座標軸方向旋轉為特徵點的方向，以特徵點為中心取視窗，透過高斯加權增強特徵點附近畫素梯度方向資訊的貢獻，即在4 × 4的小塊上計算梯度方向直方圖（取8個方向），計算梯度方向累加值，形成種子點，構成4× 4 × 8= 128維特徵向量。SIFT特徵對旋轉、尺度縮放、亮度變化保持不變性，對視角變化、仿射變換、噪聲也保持一定程度的穩定性。

對描述子的進一步處理

［去除光照影響］將特徵向量的長度歸一化

［減弱扭曲影響］仿射變換

關於SIFT演算法的幾個問題

本小節主要參考：

http：//www。

cnblogs。com/ronny/p/402

8776。html#3962356

第一個問題：第一組第一層影象的生成

這裡要分兩種情況：其一是把第一組的索引定為0；其二是把第一組的索引定為-1；我們先考慮第一組索引為0的情況，我們知道第一組第一層的影象是由原影象與

$\sigma_0$

（一般設定為1。6）的高斯濾波器卷積生成，為了影象反走樣的需要，通常假設輸入影象是經過高斯平滑處理的，其值為

$\sigma_n=0.5$

，即半個像元。

也就是說我們採集到的影象I（x，y），已經被

$\sigma=\sigma_n=0.5$

的高斯濾波器平滑過了。所以我們不能直接對I（x，y）用

$\sigma_0$

的高斯濾波器平滑，而應該用

$\sigma=\sigma_{20}−\sigma_{2n}$

的高斯濾波器去平滑I（x，y），即

$\text{FirstLayer}(x,y) = I_s(x,y)* G(x,y,\sqrt{\sigma_0^2 - (2\sigma_n)^2})$

其中FirstLayer（x，y）表示整個尺度空間為第1組第1層的影象。

現在我們來考慮把第一組的索引定為-1的情況。那麼首先第一個問題便是為什麼要把索引定為-1；如果索引為0，整個尺度空間的第1組的第1層影象已經是由原影象模糊生成的了，那麼也就是說已經丟失了細節資訊，那麼原影象我們完全沒有利用上。基於這種考慮，我們先將影象放大2倍，這樣原影象的細節就隱藏在了其中。

由上面一種情況分析，我們已經知識了I（x，y）看成是已經被

$\sigma_n=0.5$

模糊過的影象，那麼將I（x，y）放大2倍後得到Is（x，y）

，則可以看為是被

$2\sigma_n=1$

的高斯核模糊過的影象。那麼由Is

生成第1組第1層的影象用的高斯濾波器的

$\sigma = \sqrt{\sigma_0^2 -(2 \sigma_n)^2}$

，即

$\text{FirstLayer}(x,y) = I_s(x,y)\otimes G(x,y,\sqrt{\sigma_0^2 - (2\sigma_n)^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,\sigma ),G(x,y,k\sigma ),G(x,y,k^2\sigma ),G(x,y,k^3\sigma ),G(x,y,k^4\sigma )$

其中的

$G(x,y,k\sigma ),G(x,y,k^2\sigma ),G(x,y,k^3\sigma )$

這三張影象是我們用來查詢區域性極值點的影象。那麼我們則需要S+3 = 6張尺度空間影象來生成上面那些高斯差分影象，它們分別為：

$L(x,y,\sigma ),L(x,y,k\sigma ),L(x,y,k^2\sigma ),L(x,y,k^3\sigma ),L(x,y,k^4\sigma )，L(x,y,k^5\sigma )$

從上面的分析，我們知道對於尺度空間來說，我們一共需要S+3層影象來構建出來S+2層高斯差分影象。所以，如果整個尺度空間一共有n組，每組有S+3層影象，就共有n（S+3）張尺度影象。

第三個問題：為什麼取上一張的倒數第3張影象隔行取樣後作為下一組的第一張影象？

對於尺度空間第o組，第s層的影象，它的尺度為

$\sigma = \sigma_o k^{o+s/S},k =1/2,o\in[0,1,2,\dots,\text{nOctave}-1],s\in[0,1,2,\dots,S+2]$

那麼我們從第0組開始，看它各層的尺度。

第0組

$\sigma_0 \to 2^{1/3}\sigma_0\to 2^{2/3}\sigma_0\to 2^{3/3}\sigma_0\to 2^{4/3}\sigma_0\to 2^{5/3}\sigma_0$

第1組

$2\sigma_0\to 2*2^{1/3}\sigma_0\to 2*2^{2/3}\sigma_0\to 2*2^{3/3}\sigma_0\to 2*2^{4/3}\sigma_0\to 2*2^{5/3}\sigma_0$

……

可以看出，第1組的第0層影象恰好與第0組的倒數第三幅影象一致，尺度都為

$2\sigma_0$

，所以我們不需要再根據原圖來重新卷積生成每組的第0張影象，只需採用上一層的倒數第3張來降取樣即可。

我們也可以繼續分析，第0組尺度空間得到的高斯差分影象的尺度為：

$\sigma_o\to 2^{1/3}\sigma_0\to 2^{2/3}\sigma_0\to 2^{3/3}\sigma_0\to 2^{4/3}\sigma_0$

而第1組尺度空間得到的高斯差分影象的尺度為

$2\sigma_0\to 2*2^{1/3}\sigma_0\to 2*2^{2/3}\sigma_0\to 2*2^{3/3}\sigma_0\to 2*2^{4/3}\sigma_0$

：如果我們把它們的中間三項取出來拼在一起，則尺度為：

$2^{1/3}\sigma_0\to 2^{2/3}\sigma_0\to 2^{3/3}\sigma_0\to 2*2^{1/3}\sigma_0\to 2*2^{2/3}\sigma_0\to 2*2^{3/3}\sigma_0$

正好連續。這一效果帶來的直接的好處是

在尺度空間的極值點確定過程中，我們不會漏掉任何一個尺度上的極值點，而是能夠綜合考慮量化的尺度因子。

第四個問題：如何用第i-1層的影象生成第i層的影象？

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

基於上面兩點原因，所以

對於每一組的第i+1層的影象，都是由第i層的影象和一個相對尺度的高斯濾波器卷積生成。

我們首先考慮第0組，它們的第i+1層影象與第i層影象之間的相對尺度為

$\text{SigmaDiff}_i=\sqrt{(\sigma_0k^{i+1})^2 – (\sigma_0 k^i)^2}$

為了保持尺度的連續性，後面的每一組都用這樣一樣相對尺度（SIFT實際程式碼裡是這樣做的）。這裡有一個猜測，比如說尺度為2σ0

的這一組，第i層與第i+1層之間的相對尺度計算的結果應該為

$\sqrt{(2\sigma_0k^{i+1})^2 – (2\sigma_0 k^i)^2} = 2\cdot \text{SigmaDiff}_i$

可是，依然可以用

$\text{SigmaDiff}{}_i$

，是因為這一層被降維了。

第五個問題：離散空間得到的極值點與連續空間的極值點之間的差別

為了尋找尺度空間的極值點，每一個取樣點要和它所有的相鄰點進行比較，看其是否比它的影象域和尺度域的相鄰點大或小。對於其中的任意一個檢測點都要和它同尺度的8個相鄰點和上下相鄰尺度對應的9× 2個點共26個點比較，以確保在尺度空間和二維影象位置空間都檢測到極值點。也就是，比較是在一個3× 3× 3的立方體內進行。

搜尋過程從每組的第二層開始，以第二層為當前層，對第二層的DoG影象中的每個點取一個3× 3× 3的立方體，立方體上下層為第一層與第三層。這樣，搜尋得到的極值點既有位置座標（DoG的影象座標），又有空間尺度座標（層座標）。當第二層搜尋完成後，再以第三層作為當前層，其過程與第二層的搜尋類似。當S=3時，每組裡面要搜尋3層。

極值點的搜尋是在離散空間中進行的，檢測到的極值點並不是真正意義上的極值點。利用已知的離散空間點插值到連續空間極值點的方法，就是子像元插值方法。

SURF演算法

SURF特徵（Speeded Up Robust Features，加速魯棒性特徵）是對SIFT特徵的進一步最佳化，基於Hessian矩陣構造金字塔尺度空間，利用箱式濾波器（box filter）簡化二維高斯濾波，不需要再進行降取樣；透過Haar小波特徵設定特徵點主方向，這樣構建的特徵描述子就是64維的；surf構造的金字塔影象與sift有很大不同，就是因為這些不同才加快了其檢測的速度。Sift採用的是DOG影象，而surf採用的是Hessian矩陣行列式近似值影象，也寫作DOH運算元。

$H= \begin{pmatrix} L_{XX} & L_{XY} \\ L_{XY} & L_{YY} \end{pmatrix}\Rightarrow \text{det}H=L_{XX} L_{YY}-0.9 L_{XY}^2$

由於求Hessian時要先高斯平滑，然後求二階導數，這在離散的畫素點是用模板卷積形成的，這2種操作合在一起用一個模板代替就可以了。相比於sift演算法的高斯金字塔構造過程，surf演算法速度有所提高。在sift演算法中，每一組（octave）的影象大小是不一樣的，下一組是上一組影象的降取樣（1/4大小）；在每一組裡面的幾幅影象中，他們的大小是一樣的，不同的是他們採用的尺度$\sigma$不同。而且在模糊的過程中，他們的高斯模板大小總是不變的，只是尺度$\sigma$改變。對於surf演算法，影象的大小總是不變的，改變的只是高斯模糊模板的尺寸，當然，尺度也是在改變的，但不需要降取樣過程，節省時間。

為了保證旋轉不變性，在SURF中，不統計其梯度直方圖，而是統計特徵點領域內的Harr小波特徵。即以特徵點為中心，計算半徑為6s（S為特徵點所在的尺度值）的鄰域內，統計60度扇形內所有點在x（水平）和y（垂直）方向的Haar小波響應總和（Haar小波邊長取4s），並給這些響應值賦高斯權重係數，使得靠近特徵點的響應貢獻大，而遠離特徵點的響應貢獻小，然後60度範圍內的響應相加以形成新的向量，遍歷整個圓形區域，選擇最長向量的方向為該特徵點的主方向。在特徵點周圍取正方形框，方框的方向為特徵點主方向，把方框分為16個子區域，在每個子區域統計水平方向和垂直方向的haar小波特徵，在每個子區域計算haar小波特徵的水平方向值之和，水平方向絕對值之和、垂直方向值之和、垂直方向絕對值之和，構成16× 4=64維特徵向量

在完成採集後，還需要建立影象的特徵點資料庫，每個特徵點的資料結構包括：位置座標、尺度、方向、特徵向量（128或64維）；為新影象的每個特徵點在資料庫中逐個匹配，即根據特徵向量的歐氏距離在資料庫中尋找其最近鄰和次近鄰特徵點，若最近鄰距離或次近鄰距離大於某一闕值，則特徵匹配成功。

綜上所述，可知SURF採用Henssian矩陣獲取影象區域性最值還是十分穩定的，但是在求主方向階段太過於依賴區域性區域畫素的梯度方向，有可能使得找到的主方向不準確，後面的特徵向量提取以及匹配都嚴重依賴於主方向，即使不大偏差角度也可以造成後面特徵匹配的放大誤差，從而匹配不成功；另外影象金字塔的層取得不足夠緊密也會使得尺度有誤差，後面的特徵向量提取同樣依賴相應的尺度，在這個問題上我們只能採用折中解決方法：取適量的層然後進行插值。

SIFT特徵與SURF特徵的比較：

［構建影象金字塔］ SIFT特徵利用不同尺寸的影象與高斯差分濾波器卷積；SURF特徵利用原圖片與不同尺寸的方框濾波器卷積。

［特徵描述子］ SIFT特徵有4×4×8=128維描述子，SURF特徵有4×4×4=64維描述子

［特徵點檢測方法］ SIFT特徵先進行非極大抑制，再去除低對比度的點，再透過Hessian矩陣去除邊緣響應過大的點；SURF特徵先利用Hessian矩陣確定候選點，然後進行非極大抑制

［特徵點主方向］ SIFT特徵在正方形區域內統計梯度幅值的直方圖，直方圖最大值對應主方向，可以有多個主方向；SURF特徵在圓形區域內計算各個扇形範圍內x、y方向的haar小波響應，模最大的扇形方向作為主方向

ORB特徵

本節主要參考了：

https：//

blog。csdn。net/zouzoupao

pao229/article/details/52625678

ORB特徵是將FAST特徵點的檢測方法與BRIEF特徵描述子結合起來，並在它們原來的基礎上做了改進與最佳化。ORB演算法的速度大約是SIFT的100倍，是SURF的10倍。

ORB（Oriented FAST and Rotated BRIEF）是一種快速特徵點提取和描述的演算法。這個演算法是由Ethan Rublee， Vincent Rabaud， Kurt Konolige以及Gary R。Bradski在2011年一篇名為“ORB：An Efficient Alternative to SIFT or SURF”的文章中提出。ORB演算法分為兩部分，分別是特徵點提取和特徵點描述。特徵提取是由FAST（Features from Accelerated Segment Test）演算法發展來的，特徵點描述是根據BRIEF（Binary Robust IndependentElementary Features）特徵描述演算法改進的。

ORB使用ID3演算法訓練一個決策樹，將特徵點圓周上的16個畫素輸入決策樹中，以此來篩選出最優的FAST特徵點。

接著，非極大值抑制去除區域性較密集特徵點。使用非極大值抑制演算法去除臨近位置多個特徵點的問題。為每一個特徵點計算出其響應大小。計算方式是特徵點P和其周圍16個特徵點偏差的絕對值和。在比較臨近的特徵點中，保留響應值較大的特徵點，刪除其餘的特徵點。

現在我們可以建立金字塔，來實現特徵點的多尺度不變性。設定一個比例因子scaleFactor（opencv預設為1。2）和金字塔的層數nlevels（opencv預設為8）。將原影象按比例因子縮小成nlevels幅影象。縮放後的影象為：I’= I/scaleFactork（k=1，2，…， nlevels）。nlevels幅不同比例的影象提取特徵點總和作為這幅影象的oFAST特徵點。

ORB演算法提出使用矩（moment）法來確定FAST特徵點的方向。也就是說透過矩來計算特徵點以r為半徑範圍內的質心，特徵點座標到質心形成一個向量作為該特徵點的方向。矩定義如下：

$m_{pq}=\sum x^p y^q I(x,y)$

其中，I（x，y）為影象灰度表示式。該矩的質心為

$(m_{10}/m_{00},m_{01}/m_{00})$

假設角點座標為O，則向量的角度即為該特徵點的方向。計算公式

$\theta=\arctan (m_{10}/m_{01})$

rBRIEF特徵描述是在BRIEF特徵描述的基礎上加入旋轉因子改進的。BRIEF演算法計算出來的是一個二進位制串的特徵描述符。它是在一個特徵點的鄰域內，選擇n對畫素點pi、qi（i=1，2，…，n）。然後比較每個點對的灰度值的大小。如果I（pi）> I（qi），則生成二進位制串中的1，否則為0。所有的點對都進行比較，則生成長度為n的二進位制串。一般n取128、256或512，opencv預設為256。另外，值得注意的是為了增加特徵描述符的抗噪性，演算法首先需要對影象進行高斯平滑處理。

在旋轉不是非常厲害的影象裡，用BRIEF生成的描述子的匹配質量非常高，作者測試的大多數情況中都超越了SURF。但在旋轉大於30°後，BRIEF的匹配率快速降到0左右。BRIEF的耗時非常短，在相同情形下計算512個特徵點的描述子時，SURF耗時335ms，BRIEF僅8。18ms；匹配SURF描述子需28。3ms，BRIEF僅需2。19ms。在要求不太高的情形下，BRIEF描述子更容易做到實時。

改進BRIEF演算法—rBRIEF（Rotation-AwareBrief）

（1）steered BRIEF（旋轉不變性改進）

在使用oFast演算法計算出的特徵點中包括了特徵點的方向角度。假設原始的BRIEF演算法在特徵點SxS（一般S取31）鄰域內選取n對點集。經過旋轉角度θ旋轉，得到新的點對，在新的點集位置上比較點對的大小形成二進位制串的描述符。這裡需要注意的是，在使用oFast演算法是在不同的尺度上提取的特徵點。因此，在使用BRIEF特徵描述時，要將影象轉換到相應的尺度影象上，然後在尺度影象上的特徵點處取SxS鄰域，然後選擇點對並旋轉，得到二進位制串描述符。

（2）rBRIEF-改進特徵點描述子的相關性

使用steeredBRIEF方法得到的特徵描述子具有旋轉不變性，但是卻在另外一個性質上不如原始的BRIEF演算法。是什麼性質呢，是描述符的可區分性，或者說是相關性。這個性質對特徵匹配的好壞影響非常大。描述子是特徵點性質的描述。描述子表達了特徵點不同於其他特徵點的區別。我們計算的描述子要儘量的表達特徵點的獨特性。如果不同特徵點的描述子的可區分性比較差，匹配時不容易找到對應的匹配點，引起誤匹配。

為了解決描述子的可區分性和相關性的問題，ORB使用統計學習的方法來重新選擇點對集合。

首先建立300k個特徵點測試集。對於測試集中的每個點，考慮其31x31鄰域。這裡不同於原始BRIEF演算法的地方是，這裡在對影象進行高斯平滑之後，使用鄰域中的某個點的5x5鄰域灰度平均值來代替某個點對的值，進而比較點對的大小。這樣特徵值更加具備抗噪性。另外可以使用積分影象加快求取5x5鄰域灰度平均值的速度。

從上面可知，在31x31的鄰域內共有（31-5+1）x（31-5+1）=729個這樣的子視窗，那麼取點對的方法共有M=265356種，我們就要在這Ｍ種方法中選取256種取法，選擇的原則是這256種取法之間的相關性最小。怎麼選取呢？

在300k特徵點的每個31x31鄰域內按M種方法取點對，比較點對大小，形成一個300kxM的二進位制矩陣Q。矩陣的每一列代表300k個點按某種取法得到的二進位制數。

對Q矩陣的每一列求取平均值，按照平均值到0。5的距離大小重新對Q矩陣的列向量排序，形成矩陣T。

將T的第一列向量放到R中。

取T的下一列向量和R中的所有列向量計算相關性，如果相關係數小於設定的閾值，則將T中的該列向量移至R中。

按照上一步的方式不斷進行操作，直到R中的向量數量為256。

這就是rBRIEF演算法。

標簽：影象高斯特徵角點尺度

上一篇:8種掛麵吃法，有滋有味，只要10分鐘，一個禮拜不重樣

下一篇：實用！一文讀懂開關插座如何佈局和選購