OpenCV可以檢測圖像的主要特征，然后提取這些特征、使其成為圖像描述符，這類似於人的眼睛和大腦。這些圖像特征可作為圖像搜索的數據庫。此外，人們可以利用這些關鍵點將圖像拼接起來，組成一個更大的圖像，比如將許多圖像放在一塊，然后形成一個360度全景圖像。
這里我們將學習使用OpenCV來檢測圖像特征，並利用這些特征進行圖像匹配和搜索。我們會選取一些圖像，並通過單應性，檢測這些圖像是否在另一張圖像中。
一 特征檢測算法
有許多用於特征檢測和提取的算法，我們將會對其中大部分進行介紹。OpenCV最常使用的特征檢測和提取算法有：
- Harris：該算法用於檢測角點；
- SIFT：該算法用於檢測斑點；
- SURF：該算法用於檢測角點；
- FAST：該算法用於檢測角點；
- BRIEF：該算法用於檢測斑點；
- ORB：該算法代表帶方向的FAST算法與具有旋轉不變性的BRIEF算法；
通過以下方法進行特征匹配：
- 暴力(Brute-Force)匹配法；
- 基於FLANN匹配法；
可以采用單應性進行空間驗證。
二 特征定義
那么，究竟什么是特征呢？為什么一副圖像的某個特定區域可以作為一個特征，而其他區域不能呢？粗略的講，特征就是有意義的圖像區域，該區域具有獨特特征和易於識別性。因此角點及高密度區域都是很好的特征，而大量重復的模式或低密度區域(例如圖像中的藍色天空)則不是很好的特征。邊緣可以將圖像分為兩個區域，因此也可以看做好的特征。斑點是與周圍有很大差別的像素區域，也是有意義的特征。
大多數特征檢測算法都會涉及圖像的角點、邊和斑點的識別，也有一些涉及脊向的概念，可以認為脊向是細長物體的對稱軸，例如識別圖像中的一條路。角點和邊都好理解，那什么是斑點呢？斑點通常是指與周圍有着顏色和灰度差別的區域。在實際地圖中，往往存在着大量這樣的斑點，如一顆樹是一個斑點，一塊草地是一個斑點，一棟房子也可以是一個斑點。由於斑點代表的是一個區域，相比單純的角點，它的穩定性要好，抗噪聲能力要強，所以它在圖像配准上扮演了很重要的角色。
由於某些算法在識別和提取某類型特征的時候有較好的效果，所以知道輸入圖像是什么很重要，這樣做有利於選擇最合適的OpenCV工具包。
三 Harris檢測角點特征
在這之前其實我們已經接觸過角點檢測了，在相機標定的時候，我們就利用到了角點檢測。不過那時候沒有深入的去研究。在這里，我們將會深入原理去學習角點檢測。
下面我們從使用cornerHarris()函數講起。
cornerHarris(src, blockSize, ksize,k[,dst[,borderType]]);
- 參數詳解:
- image:輸入的單通道8位或者浮點圖像；
- blockSize：就是掃描時候窗口的大小。
- ksize：cornerHarris函數會使用Sobel算子，該參數定義了Sobel算子的中孔。簡單來說，該函數定義了角點檢測的敏感度，其值必須介於3~31之間的奇數。
- k：harris 計算響應公式中的$k$值，一般取0.04~0.06；
- borderType：像素插值方法；
函數 cornerHarris 對輸入圖像進行 Harris 邊界檢測。輸出是一幅浮點值圖像，大小與輸入圖像大小相同，浮點值越高，表明越可能是特征角點(我們可以對圖像進行閾值化)。
# -*- coding: utf-8 -*- """ Created on Mon Aug 20 20:17:34 2018 @author: lenovo """ ''' Harris角點檢測 ''' import cv2 import numpy as np img = cv2.imread('./image/cali.bmp') img = cv2.resize(img,dsize=(600,400)) gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY) gray = np.float32(gray) #角點檢測 第三個參數為角點檢測的敏感度，其值必須介於3~31之間的奇數 dst = cv2.cornerHarris(gray,3,23,0.04) print(dst.shape) #(400, 600) img[dst>0.01*dst.max()] = [0,0,255] cv2.imshow('',img) cv2.waitKey(0) cv2.destroyAllWindows()
運行結果如下：
如果我們把第三個參數改為3:，可以看到：
四 Harris檢測原理
上面我們已經通過實例演示了Harris檢測的效果，相信你對Harris角點檢測已經有了初步的認識。這里我將帶你深入了解Harris角點檢測的原理。
我們先來看一幅圖片，了解一下什么是角點？
上圖中E,F中的角我們通常稱作角點（corner points），他們具有以下特征：
- 輪廓之間的交點；
- 對於同一場景，即使視角發生變化，通常具備穩定性質的特征；
- 該點附近區域的像素點無論在梯度方向上還是其梯度幅值上有着較大變化；
harris特征角最早在paper A Combined Corner and Edge Detector中被Chris Harris & Mike Stephens提出。
Harris角點檢測的基本思想：算法基本思想是使用一個固定窗口在圖像上進行任意方向上的滑動，比較滑動前與滑動后兩種情況，窗口中的像素灰度變化程度，如果存在任意方向上的滑動，都有着較大灰度變化，那么我們可以認為該窗口中存在角點。
1、灰度變化描述
當窗口發生$[u,v]$移動時，那么滑動前與滑動后對應的窗口中的像素點灰度變化描述如下：
$$E(u,v)=\sum\limits_{(x,y)€W}w(x,y)[I(x+u,y+v)-I(x,y)]^2$$
參數解釋：
- $[u,v]$是窗口$W$的偏移量；
- $(x,y)$是窗口$W$所對應的像素坐標位置，窗口有多大，就有多少個位置；
- $I(x,y)$是像素坐標位置$(x,y)$的圖像灰度值；
- $I(x+u,y+v)$是像素坐標位置$(x+u,y+v)$的圖像灰度值；
- $w(x,y)$是窗口函數，最簡單情形就是窗口$W$內的所有像素所對應的$w$權重系數均為1.但有時候，我們會將$w(x,y)$函數設置為以窗口$W$中心為原點的二元正太分布。如果窗口$W$中心點是角點時，移動前與移動后，該點在灰度變化貢獻最大；而離窗口$W$中心(角點)較遠的點，這些點的灰度變化幾近平緩，這些點的權重系數，可以設定小值，以示該點對灰度變化貢獻較小，那么我們自然而然想到使用二元高斯函數來表示窗口函數；
我們的窗口函數通常有如下兩種形式：
根據上述表達式，當窗口在平潭區域上移動，可以想象得到，灰度不會發生什么變換。$E(u,v)=0$；如果窗口處在紋理比較豐富的區域上滑動，那么灰度變化會很大。算法最終思想就是計算灰度發生較大變化時所對應的位置，當然這個較大是指任意方向上的滑動，並非單指某個方向。
2、$E(u,v)$化簡
首先需要了解泰勒公式，任何一個函數表達式，均可有泰勒公式進行展開，以逼近原函數，我們可以對下面函數進行一階展開(如果對泰勒公式忘記了，可以翻翻本科所學的高等數學)。
$$f(x+u,y+v)≈f(x,y)+uf_x(x,y)+vf_y(x,y)$$
那么
$$\sum\limits_{(x,y)€W}w(x,y)[I(x+u,y+v)-I(x,y)]^2$$
$$≈\sum\limits_{(x,y)€W}w(x,y)[I(x,y)+uI_x+vI_y-I(x,y)]^2$$
$$=\sum\limits_{(x,y)€W}w(x,y)[u^2I_x^2+2uvI_xI_y+v^2I_y^2]$$
$$=\sum\limits_{(x,y)€W}w(x,y)\begin{bmatrix}u & v\end{bmatrix}\begin{bmatrix} I_x^2 & I_xIy \\ I_xI_y & I_y^2\end{bmatrix}\begin{bmatrix}u \\ v\end{bmatrix}$$
$$=\begin{bmatrix}u & v\end{bmatrix}(\sum\limits_{(x,y)€W}w(x,y)\begin{bmatrix} I_x^2 & I_xIy \\ I_xI_y & I_y^2\end{bmatrix})\begin{bmatrix}u \\ v\end{bmatrix}$$
所以$E(u,v)$表達式可以更新為：
$$E(u,v)=\begin{bmatrix}u & v\end{bmatrix}M\begin{bmatrix}u \\ v\end{bmatrix}$$
其中：$M=\sum\limits_{(x,y)€W}w(x,y)\begin{bmatrix} I_x^2 & I_xIy \\ I_xI_y & I_y^2\end{bmatrix}$，$I_x$,$I_y$分別為窗口內像素點$(x,y)$在$x$方向上和$y$方向上的梯度值。
3、矩陣$M$的關鍵性
難道我們是直接求上述的$E(u,v)$值來判斷角點嗎？Harris角點檢測並沒有這樣做，而是通過對窗口內的每個像素的$x$方向上的梯度與$y$方向上的梯度進行統計分析。這里以$I_x$和$I_y$為坐標軸，因此每個像素的梯度坐標可以表示成$(I_x,I_y)$。針對平坦區域，邊緣區域以及角點區域三種情形進行分析：
下圖是對這三種情況窗口中的對應像素的梯度分布進行繪制：
不知道大家有沒有注意到這三種區域的特點:
- 平坦區域上的每個像素點所對應的$(I_x,I_y)$坐標分布在原點附近，其實也很好理解，針對平坦區域的像素點，他們的梯度方向雖然各異，但是其幅值都不是很大，所以均聚集在原點附近；
- 邊緣區域有一坐標軸分布較散，至於是哪一個坐標上的數據分布較散不能一概而論，這要視邊緣在圖像上的具體位置而定，如果邊緣是水平或者垂直方向，那么$I_y$軸方向或者$I_x$方向上的數據分布就比較散；
- 角點區域的$x$、$y$方向上的梯度分布都比較散。
我們是不是可以根據這些特征來判斷哪些區域存在角點呢？
上面我們已經計算出了$M$矩陣:
$M=\sum\limits_{(x,y)€W}w(x,y)\begin{bmatrix} I_x^2 & I_xIy \\ I_xI_y & I_y^2\end{bmatrix}=\begin{bmatrix} A & C \\ C & B \end{bmatrix}$
我們可以將$E(u,v)$近似為二項函數：
$$E(u,v)=Au^2+2Cuv+Bv^2$$
其中
$$A=\sum\limits_{(x,y)€W}w(x,y)*I_x^2$$
$$B=\sum\limits_{(x,y)€W}w(x,y)*I_y^2$$
$$C=\sum\limits_{(x,y)€W}w(x,y)*I_xI_y$$
二次項函數本質上就是一個橢圓函數。橢圓的長和寬是由$M$的特征值$λ_1,λ_2$決定的(橢圓的長短軸正是矩陣$M$特征值平方根的倒數)，橢圓的方向是由$M$的特征向量決定的，橢圓方程為：
$$\begin{bmatrix} u & v \end{bmatrix}M\begin{bmatrix} u \\ v \end{bmatrix}=1$$
如果使用橢圓進行數據集表示，則繪制圖示如下：
雖然我們利用$E(u,v)$來描述角點的基本思想，然而最終我們僅僅使用的是矩陣$M$。讓我們看看矩陣$M$形式，是不是跟協方差矩陣形式很像，像歸像，但是還是有些不同，哪兒不同？一般協方差矩陣對應維的隨機變量需要減去該維隨機變量的均值：
把$I_x$,$I_y$看成兩個字段，假設窗口內有m個像素點，也就是等價於有m個樣本，我們先計算每個字段的均值：
$$\bar{I_x}=\sum\limits_{i=1}^{m}I_{xi}$$
$$\bar{I_y}=\sum\limits_{i=1}^{m}I_{yi}$$
我們仍然使用$(I_{xi},I_{yi})$表示樣本$(I_{xi},I_{yi})$去均值后的值，則由這m個樣本組成的矩陣：
$$X=\begin{bmatrix} I_{x1} & I_{x2} & ... & I_{xm} \\ I_{y1} & I_{y2} & ... & I_{ym}\end{bmatrix}$$
則對應協方差矩陣可以寫成：$$C=\frac{1}{m}XX^T=\frac{1}{m}\sum\limits_{i=1}^{m}\begin{bmatrix} I_x^2 & I_xIy \\ I_xI_y & I_y^2 \end{bmatrix}$$
但矩陣$M$中並沒有這樣做，所以在矩陣$M$里，我們先進行各維的均值化處理，那么各維所對應的隨機變量的均值為0，協方差矩陣就大大簡化了，簡化的最終結果就是矩陣$M$，是否明白了(注意為了簡化運算，我們先假設$M$矩陣中的權重系數$w(x,y)=1$，並且省略掉了求均值)我們的目的是分析數據的主要成分，相信了解PCA原理的，應該都了解均值化的作用。$$M=\sum\limits_{(x,y)€W}\begin{bmatrix} I_x^2 & I_xIy \\ I_xI_y & I_y^2\end{bmatrix}$$
如果我們對協方差矩陣$M$進行對角化，很明顯，其對角線就是各個字段的方差，這點大家應該明白吧？不明白的話可以復習下PCA原理。
- 如果兩個字段$(I_x,I_y)$所對應的特征值都比較大，說明什么？ 像素點的梯度分布比較散，梯度變化程度比較大，符合角點在窗口區域的特點；
- 如果是平坦區域，那么像素點的梯度所構成的點集比較集中在原點附近，因為窗口區域內的像素點的梯度幅值非常小，此時矩陣$M$的對角化的兩個特征值比較小；
- 如果是邊緣區域，在計算像素點的$x$、$y$方向上的梯度時，邊緣上的像素點的某個方向的梯度幅值變化比較明顯，另一個方向上的梯度幅值變化較弱，其余部分的點都還是集中原點附近，這樣$M$對角化后的兩個特征值理論應該是一個比較大，一個比較小，當然對於邊緣這種情況，可能是呈45°的邊緣，致使計算出的特征值並不是都特別的大，總之跟含有角點的窗口的分布情況還是不同的。
注：$M$為協方差矩陣，需要大家自己去理解下，窗口中的像素集構成一個矩陣（$X€R^{2×m}$，假設這里有m個像素點），使用該矩陣乘以該矩陣的轉置，即是協方差矩陣。
因此可以得出下列結論：
- 特征值都比較大時，即窗口中含有角點；
- 特征值一個較大，一個較小，窗口中含有邊緣；
- 特征值都比較小，窗口處在平坦區域；
4、如何度量角點響應
通常用下面表達式進行度量，對每一個窗口計算得到一個分數$R$，根據$R$的大小來判定窗口內是否存在harris特征角。分數$R$根據下面公式計算得到：
$$R=det(M)-k(trace(M))^2$$
$$det(M)=λ_1λ_2$$
$$trace(M)=λ_1+λ_2$$
這里$λ_1,λ_2$是矩陣$M$的2個特征值，$k$是一個指定值，這是一個經驗參數，需要實驗確定它的合適大小，通常它的值在0.04和0.06之間,它的存在只是調節函數的形狀而已。
$R$取決於$M$的特征值，對於角點$|R|$很大，平坦的區域$|R|$很小，邊緣的$R$為負值;
但是為什么會使用這樣的表達式呢？一下子是不是感覺很難理解？其實也不難理解，函數表達式一旦出來，我們就可以繪制它的圖像，而這個函數圖形正好滿足上面幾個區域的特征。 通過繪制函數圖像，直觀上更能理解。繪制的$R$函數圖像如下：
所以說難點不在於理解這個函數表達式，而在於如何創造出這個函數表達式。Harris也許對很多函數模型非常了解，對於創造出這樣的一個函數表達式，易如反掌，當然在我們看來感覺是很了不起的，那是因為我們見過的函數模型太少。
最后設定R的閾值，進行角點判斷。當然其中還有些后處理步驟就不再說了，比如說角點的極大值抑制等。
非極大值抑制的原理，在一個窗口內，如果有很多角點則用值最大的那個角點，其他的角點都刪除。
注意：有些朋友在問是如何通過矩陣判斷角點的？ 其實上面，我們已經推導出$E(u,v)$的表達式，大家看看這個表達式有什么特征，其中矩陣M是實對稱矩陣，那么E表達式其實就是二次型，對於二次型想必大家會有印象，$u$,$v$代表窗口滑動方向以及滑動量，$E$代表灰度變化，通過矩陣M進行特征值求解，而特征值所對應的特征向量即為灰度變化方向。如果兩個特征值較大，則表示有兩個方向灰度變化較快。所以可以直接通過求解$M$的特征值進行角點判斷.
五 自己實現Harris角點檢測
Harris角點檢測可以分為5個步驟：
1、計算圖像$I(x,y)$在$x$和$y$兩個方向的梯度$I_x$,$I_y$;
$$I_x=\frac{∂I}{∂x}=I×\begin{bmatrix} -1 & 0 & 1 \end{bmatrix}$$
$$I_y=\frac{∂I}{∂y}=I×\begin{bmatrix} -1 & 0 & 1 \end{bmatrix}^T$$
2、計算圖像兩個方向梯度的乘積;
$$I_x^2=I_x*I_x$$
$$I_y^2=I_y*I_y$$
$$I_xI_y=I_x*I_y$$
3、使用高斯函數對$I_x^2$、$I_y^2$、$I_xI_y$進行高斯加權(取σ=2,ksize=3)，計算中心點為$(x,y)$的窗口$W$對應的矩陣$M$；
$$A=\sum\limits_{(x,y)€W}g(I_x^2)=\sum\limits_{(x,y)€W}I_x^2*w(x,y)$$
$$B=\sum\limits_{(x,y)€W}g(I_y^2)=\sum\limits_{(x,y)€W}I_y^2*w(x,y)$$
$$C=\sum\limits_{(x,y)€W}g(I_xI_y)=\sum\limits_{(x,y)€W}I_xI_y*w(x,y)$$
其中$M=\begin{bmatrix} A& C \\C & B\end{bmatrix}=\sum\limits_{(x,y)€W}w(x,y)\begin{bmatrix} I_x^2 & I_xIy \\ I_xI_y & I_y^2\end{bmatrix}$.
4、計算每個像素點$(x,y)處的$Harris響應值$R$；
$$R=det(M)-k(trace(M))^2$$
5、過濾大於某一閾值$t$的$R$值；
$$R=\{R:det(M)-k(trace(M))^2 > t \}$$
如果需要在3×3或者5×5的鄰域進行非最大值抑制，則局部最大值點即為圖像中的角點。
以上所有公式*表示逐像素點乘。
# -*- coding: utf-8 -*- """ Created on Mon Aug 20 20:17:34 2018 @author: lenovo """ ''' Harris角點檢測 ''' import cv2 import numpy as np def test(): ''' 調用系統庫函數進行測試 ''' img = cv2.imread('./image/cali.bmp') img = cv2.resize(img,dsize=(600,400)) gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY) gray = np.float32(gray) #角點檢測 第三個參數為角點檢測的敏感度，其值必須介於3~31之間的奇數 dst = cv2.cornerHarris(gray,3,3,0.04) print(dst.shape) #(400, 600) img[dst>0.01*dst.max()] = [0,0,255] cv2.imshow('',img) cv2.waitKey(0) cv2.destroyAllWindows() def harris_detect(img,ksize=3): ''' 自己實現角點檢測 params: img:灰度圖片 ksize：Sobel算子窗口大小 return： corner：與源圖像一樣大小，角點處像素值設置為255 ''' k = 0.04 #響應函數k threshold = 0.01 #設定閾值 WITH_NMS = False #是否非極大值抑制 #1、使用Sobel計算像素點x,y方向的梯度 h,w = img.shape[:2] #Sobel函數求完導數后會有負值，還有會大於255的值。而原圖像是uint8，即8位無符號數，所以Sobel建立的圖像位數不夠，會有截斷。因此要使用16位有符號的數據類型，即cv2.CV_16S。 grad = np.zeros((h,w,2),dtype=np.float32) grad[:,:,0] = cv2.Sobel(img,cv2.CV_16S,1,0,ksize=3) grad[:,:,1] = cv2.Sobel(img,cv2.CV_16S,0,1,ksize=3) #2、計算Ix^2,Iy^2,Ix*Iy m = np.zeros((h,w,3),dtype=np.float32) m[:,:,0] = grad[:,:,0]**2 m[:,:,1] = grad[:,:,1]**2 m[:,:,2] = grad[:,:,0]*grad[:,:,1] #3、利用高斯函數對Ix^2,Iy^2,Ix*Iy進行濾波 m[:,:,0] = cv2.GaussianBlur(m[:,:,0],ksize=(ksize,ksize),sigmaX=2) m[:,:,1] = cv2.GaussianBlur(m[:,:,1],ksize=(ksize,ksize),sigmaX=2) m[:,:,2] = cv2.GaussianBlur(m[:,:,2],ksize=(ksize,ksize),sigmaX=2) m = [np.array([[m[i,j,0],m[i,j,2]],[m[i,j,2],m[i,j,1]]]) for i in range(h) for j in range(w)] #4、計算局部特征結果矩陣M的特征值和響應函數R(i,j)=det(M)-k(trace(M))^2 0.04<=k<=0.06 D,T = list(map(np.linalg.det,m)),list(map(np.trace,m)) R = np.array([d-k*t**2 for d,t in zip(D,T)]) #5、將計算出響應函數的值R進行非極大值抑制，濾除一些不是角點的點，同時要滿足大於設定的閾值 #獲取最大的R值 R_max = np.max(R) #print(R_max) #print(np.min(R)) R = R.reshape(h,w) corner = np.zeros_like(R,dtype=np.uint8) for i in range(h): for j in range(w): if WITH_NMS: #除了進行進行閾值檢測 還對3x3鄰域內非極大值進行抑制(導致角點很小，會看不清) if R[i,j] > R_max*threshold and R[i,j] == np.max(R[max(0,i-1):min(i+2,h-1),max(0,j-1):min(j+2,w-1)]): corner[i,j] = 255 else: #只進行閾值檢測 if R[i,j] > R_max*threshold : corner[i,j] = 255 return corner if __name__=='__main__': img = cv2.imread('./image/cali.bmp') img = cv2.resize(img,dsize=(600,400)) #轉換為灰度圖像 gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY) dst = harris_detect(gray) print(dst.shape) #(400, 600) img[dst>0.01*dst.max()] = [0,0,255] cv2.imshow('',img) cv2.waitKey(0) cv2.destroyAllWindows()
運行效果如下：
六 Harries角點的性質
1、參數$k$對角點檢測的影響
假設已經得到了矩陣$M$的特征值$λ_1≥λ_2≥0$，令$λ_2=αλ_1$，$0≤α≤1$。由特征值與矩陣$M$的直跡和行列式的關系可得：
$$detM=\prod\limits_{i}λ_i$$ $$trace(M)=\sum\limits_{i}λ_i$$
從而可以得到角點的響應：
$$R=λ_1λ_2-k(λ_1+λ_2)^2=λ_1^2(α-k(1+α)^2)$$
假設$R≥0$，則有：
$$0≤k≤\frac{α}{(1+α)^2}≤0.25$$
對於較小的$α$值，$R≈λ_1^2(α-k),α<k$。
由此，可以得出這樣的結論，增大$k$的值，降低角點檢測的靈敏度，減少被檢測角點的數量；減少$k$值，增加角點檢測的靈敏度，增加被檢測角點的數量。
2、Harris角點檢測算子對亮度和對比度的變化不靈敏
這是因為在進行Harris角點檢測時，使用了微分算子對圖像進行微分運算，而微分運算對圖像密度的拉升或收縮和對亮度的抬高或下降不敏感。換言之，對亮度和對比度的仿射變換並不改變Harris響應的極值點出現的位置，但是，由於閾值的選擇，可能會影響角點檢測的數量。
左圖表示亮度變化，右圖表示對比度變化。
3、Harris角點檢測算子具有旋轉不變性
Harris角點檢測算子使用的是角點附近的區域灰度二階矩矩陣。而二階矩矩陣可以表示成一個橢圓，橢圓的長短軸正是二階矩矩陣特征值平方根的倒數。當特征橢圓轉動時，特征值並不發生變化，所以判斷角點響應值$R$也不發生變化，由此說明Harris角點檢測算子具有旋轉不變性。
4.、Harris角點檢測算子不具有尺度不變性
如下圖所示，當圖像被縮小時，在檢測窗口尺寸不變的前提下，在窗口內所包含圖像的內容是完全不同的。左側的圖像可能被檢測為邊緣或曲線，而右側的圖像則可能被檢測為一個角點。
參考文獻：