特徵檢測之Harris角點檢測

特徵點又叫興趣點或者角點。常被用於目標匹配，目標跟蹤，三維重建等應用中。點特徵主要指圖像中的明顯點，如突出的角點、邊緣端點、極值點等等。特徵點具備如下三種特性：旋轉不變性；光照不變性；視角不變性。用於點特徵提取的算子稱爲興趣點提取（檢測）算子。經常使用的有Harris角點檢測；FAST特徵檢測；SIFT特徵檢測；SURF特徵檢測。

預備數學

實對稱矩陣

若是有n階矩陣A，其矩陣的元素都爲實數，且矩陣A的轉置等於其自己，則稱A爲實對稱矩陣。

性質：

• 實對稱矩陣Ａ的不一樣特徵值對應的特徵向量是正交的；

• n階實對稱矩陣Ａ必可類似對角化，且類似對角陣上的元素即爲矩陣自己特徵值。

對角化：

優勢：滿秩的對角方陣，能夠直接看出特徵值，特徵向量等重要特徵。

卷積算子—Sobel算子

主要用於邊緣檢測，分別對水平和垂直方向上的邊緣敏感。對圖像任意一點使用會產生對應的梯度矢量或者其法矢量。對噪聲具備平滑抑制做用，可是獲得的邊緣較粗，且可能出現僞邊緣。該算子包含兩組3x3的矩陣，分別爲橫向及縱向，將之與圖像做平面卷積，便可分別得出橫向及縱向的亮度差分近似值。若是以A表明原始圖像，及分別表明經橫向及縱向邊緣檢測的圖像，其公式以下:

梯度計算公式

梯度方向計算公式

若是等於零，表明圖像該處擁有縱向邊緣，左邊比右邊要暗。

基礎知識

一、角點

使用一個滑動窗口在下面三幅圖中滑動，能夠得出如下結論：

• 左圖表示一個平坦區域，在各方向移動，窗口內像素值均沒有太大變化；

• 中圖表示一個邊緣特徵（Edges），若是沿着水平方向移動(梯度方向)，像素值會發生跳變；若是沿着邊緣移動(平行於邊緣) ，像素值不會發生變化；

• 右圖表示一個角（Corners），無論你把它朝哪一個方向移動，像素值都會發生很大變化。

因此，右圖是一個角點。

二、角點類型 

下圖展現了不一樣角點的類型，能夠發現：若是使用一個滑動窗口以角點爲中心在圖像上滑動，存在朝多個方向上的移動會引發該區域的像素值發生很大變化的現象。

三、圖像梯度 

「像素值發生很大變化」這一現象能夠用圖像梯度進行描述。在圖像局部內，圖像梯度越大表示該局部內像素值變化越大（灰度的變化率越大）。而圖像的梯度在數學上可用微分或者導數來表示。對於數字圖像來講，至關因而二維離散函數求梯度，並使用差分來近似導數：

在實際操做中，對圖像求梯度一般是考慮圖像的每一個像素的某個鄰域內的灰度變化，所以一般對原始圖像中像素某個鄰域設置梯度算子，而後採用小區域模板進行卷積來計算，經常使用的有Prewitt算子、Sobel算子、Robinson算子、Laplace算子等。

Harris角點檢測算法

一、算法思想 

算法的核心是利用局部窗口在圖像上進行移動，判斷灰度是否發生較大的變化。若是窗口內的灰度值（在梯度圖上）都有較大的變化，那麼這個窗口所在區域就存在角點。

這樣就能夠將 Harris 角點檢測算法分爲如下三步：

• 當窗口（局部區域）同時向  （水平）和 （垂直） 兩個方向移動時，計算窗口內部的像素值變化量  ；

• 對於每一個窗口，都計算其對應的一個角點響應函數 ；

• 而後對該函數進行閾值處理，若是 ，表示該窗口對應一個角點特徵。

二、算法步驟 

第一步 — 創建數學模型 

第一步是經過創建數學模型，肯定哪些窗口會引發較大的灰度值變化。讓一個窗口的中心位於灰度圖像的一個位置，這個位置的像素灰度值爲，若是這個窗口分別向  和  方向移動一個小的位移和，到一個新的位置  ，這個位置的像素灰度值就是。

就是窗口移動引發的灰度值的變化值。

設 爲位置 處的窗口函數，表示窗口內各像素的權重，最簡單的就是把窗口內全部像素的權重都設爲1，即一個均值濾波核。 固然，也能夠把  設定爲以窗口中心爲原點的高斯分佈，即一個高斯核。若是窗口中心點像素是角點，那麼窗口移動先後，中心點的灰度值變化很是強烈，因此該點權重係數應該設大一點，表示該點對灰度變化的貢獻較大；而離窗口中心（角點）較遠的點，這些點的灰度變化比較小，因而將權重係數設小一點，表示該點對灰度變化的貢獻較小。 則窗口在各個方向上移動 所形成的像素灰度值的變化量公式以下：

爲了提升計算效率，對上述公式進行簡化，利用泰勒級數展開來獲得這個公式的近似形式：對於二維的泰勒展開式公式爲：

則 爲：

其中和是的微分（偏導），在圖像中就是求和方向的梯度圖：

將

代入，可得：

提出  和  ，獲得最終的近似形式：

其中矩陣爲：

最後是把實對稱矩陣對角化處理後的結果，能夠把R當作旋轉因子，其不影響兩個正交方向的變化份量。

經對角化處理後，將兩個正交方向的變化份量提取出來，就是  和 （特徵值）。這裏利用了線性代數中的實對稱矩陣對角化的相關知識，有興趣的同窗能夠進一步查閱相關資料。

第二步—角點響應函數R 如今咱們已經獲得 的最終形式，別忘了咱們的目的是要找到會引發較大的灰度值變化的那些窗口。灰度值變化的大小則取決於矩陣，爲梯度的協方差矩陣。在實際應用中爲了可以應用更好的編程，因此定義了角點響應函數，經過斷定大小來判斷像素是否爲角點。計算每一個窗口對應的得分（角點響應函數R定義）：

其中

是矩陣的行列式，是矩陣的跡。

 和 是矩陣的特徵值， 是一個經驗常數，在範圍 (0.04, 0.06) 之間。

的值取決於的特徵值，對於角點很大，平坦的區域很小，邊緣的爲負值。

第三步—角點斷定 根據  的值，將這個窗口所在的區域劃分爲平面、邊緣或角點。爲了獲得最優的角點，咱們還可使用非極大值抑制。注意：Harris 檢測器具備旋轉不變性，但不具備尺度不變性，也就是說尺度變化可能會致使角點變爲邊緣。想要尺度不變特性的話，能夠關注SIFT特徵。由於特徵值 和 決定了  的值，因此咱們能夠用特徵值來決定一個窗口是平面、邊緣仍是角點：

• 平面:：該窗口在平坦區域上滑動，窗口內的灰度值基本不會發生變化，因此值很是小，在水平和豎直方向的變化量均較小，即 和 都較小，那麼 和 都較小；

• 邊緣：值爲負數，僅在水平或豎直方向有較大的變化量，即和 只有一個較大，也就是  或 ；

• 角點：[公式] 值很大，在水平、豎直兩個方向上變化均較大的點，即和 都較大，也就是 和 都很大。

以下圖所示：

Harris 角點檢測的結果是帶有這些分數  的灰度圖像，設定一個閾值，分數大於這個閾值的像素就對應角點。

三、算法性質 

Harris角點檢測的性質可總結以下：

1. 閾值決定角點的數量。

Harris角點檢測算子對亮度和對比度的變化不敏感（光照不變性） 在進行Harris角點檢測時，使用了微分算子對圖像進行微分運算，而微分運算對圖像密度的拉昇或收縮和對亮度的擡高或降低不敏感。換言之，對亮度和對比度的仿射變換並不改變Harris響應的極值點出現的位置，可是，因爲閾值的選擇，可能會影響角點檢測的數量。

2. Harris角點檢測算子具備旋轉不變性。

Harris角點檢測算子使用的是角點附近的區域灰度二階矩矩陣。而二階矩矩陣能夠表示成一個橢圓，橢圓的長短軸正是二階矩矩陣特徵值平方根的倒數。當特徵橢圓轉動時，特徵值並不發生變化，因此判斷角點響應值也不發生變化，由此說明Harris角點檢測算子具備旋轉不變性。

3. Harris角點檢測算子不具備尺度不變性。

尺度的變化會將角點變爲邊緣，或者邊緣變爲角點，Harris的理論基礎並不具備尺度不變性。

基於OpenCV的實現

在opencv中有提供實現 Harris 角點檢測的函數 cv2.cornerHarris，咱們直接調用的就能夠，很是方便。函數原型：cv2.cornerHarris(src, blockSize, ksize, k[, dst[, borderType]])。對於每個像素，在(blockSize x blockSize)鄰域內，計算梯度圖的協方差矩陣，而後經過上面第二步中的角點響應函數獲得結果圖。圖像中的角點能夠爲該結果圖的局部最大值。便可以獲得輸出圖中的局部最大值，這些值就對應圖像中的角點。

參數解釋：

• src - 輸入灰度圖像，float32類型

• blockSize - 用於角點檢測的鄰域大小，就是上面提到的窗口的尺寸

• ksize - 用於計算梯度圖的Sobel算子的尺寸

• k - 用於計算角點響應函數的參數k，取值範圍常在0.04~0.06之間

import cv2 as cvfrom matplotlib import pyplot as pltimport numpy as np
# detector parametersblock_size = 3sobel_size = 3k = 0.06image = cv.imread('Scenery.jpg')print(image.shape)height = image.shape[0]width = image.shape[1]channels = image.shape[2]print("width: %s height: %s channels: %s"%(width, height, channels))gray_img = cv.cvtColor(image, cv2.COLOR_BGR2GRAY)
# modify the data type setting to 32-bit floating pointgray_img = np.float32(gray_img)
# detect the corners with appropriate values as input parameterscorners_img = cv.cornerHarris(gray_img, block_size, sobel_size, k)
# result is dilated for marking the corners, not necessarykernel = cv2.getStructuringElement(cv2.MORPH_RECT,(3, 3))dst = cv.dilate(corners_img, kernel)
# Threshold for an optimal value, marking the corners in Green#image[corners_img>0.01*corners_img.max()] = [0,0,255]for r in range(height):for c in range(width):pix=dst[r,c]if pix>0.05*dst.max():cv2.circle(image,(c,r),5,(0,0,255),0)image = cv.cvtColor(image, cv2.COLOR_BGR2RGB)plt.imshow(image)plt.show()