【小林的OpenCV基礎課 18】兩種經典角點檢測演算法/萌の勝利

來自專欄 小林的CV視覺工坊

記得關注微信公眾號【小林CV工坊】哦，愛你萌

小林想給康娜找個妹妹，不知道新田雛能不能駕馭得了，畢竟一個是一千多歲的龍，一個是超能力女兒(づ￣ 3￣)づ

從這一話開始小林將帶領同學們開始OpenCV圖像特徵基本處理方法的學習。

這一話的結構：

• 基本圖像特徵
• 兩種經典的角點檢測演算法、API與Demo
• 亞像素級角點檢測、API與Demo

三種基本的圖像特徵

• 邊緣，對應下圖黑框部分
• 角點，對應下圖紅框部分
• 團塊，對應下圖藍框部分

下面介紹三種經典的角點檢測的演算法（前方公式高能，非戰鬥人員可直接跳到API和Demo部分）

• Harri演算法
• Shi-Tomasi演算法
• 亞像素級角點檢測

敲黑板了，準備面試CV相關方向的同學和準備深造CV的同學要打牢CV經典演算法的基礎，所以這些同學建議仔細閱讀演算法的原理部分。

Harri演算法

Harri演算法作為一種角點檢測的經典演算法，核心思想還是對像素進行梯度運算，總結角點處梯度的特徵，前面講到了XY角點的梯度特徵，小林就以這個特徵講解Harri演算法。

Harri演算法有兩種常見的表達式，我們先來看第一種。

• u和v：窗口偏移量
• x和y：窗口內像素坐標
• ：窗口函數，內含權重信息，常用的有權重為1和呈二元高斯正太分布的權重，意在突出像素值變化明顯的程度
• 函數 ：像素密度函數，類比與像素值

現在來處理下這個公式，我們對後面平方項的求和式 使用泰勒公式展開：

現在看第二種表達式，OpenCV中對應的函數就使用了第二種表達式運算的：

• ：矩陣行列式，
• ：矩陣的跡，
• ：Harri係數，一般取值為0.04~0.06

那麼為什麼要對矩陣M進行運算呢？從第一種表達式中可以看出，矩陣M包含著圖像的梯度信息。因此我們來看下像素的梯度分布圖：

• 團塊部分（Flat）的梯度隨呈各向異性，但梯度值非常小，因此集中在原點附近
• 邊緣部分（Linear Edge）的梯度方向是確定的，集中在某一坐標軸附近
• 角點部分（Corner）的梯度的方向也是確定的，集中在兩個坐標軸附近

再結合抽象後的圖像：

想必同學們就能理解Harri演算法的大致流程了吧。那麼問題來了，用什麼標準來界定 和 的大小呢？通常我們使用分段的思想來界定，實際使用中使用閾值，這一點小林會在Demo中印證。

小林再來總結下Harri演算法的流程：

• 求兩個方向梯度，並計算出矩陣M
• 對矩陣M計算特徵值、行列式和跡
• 根據特徵值的關係並使用閾值確定圖像特徵

API：

dst=cv2.cornerHarris(src, blockSize, ksize, k[, dst[, borderType]])

• dst：Harri演算法的輸出矩陣（輸出圖像），CV_32FC1類型，與src有同樣的尺寸
• src：輸入圖像，單通道，8位或浮點型
• blockSize：鄰域大小
• ksize：Sobel運算元的孔徑大小
• k：Harri演算法係數，即第二種表達式中的

Harri演算法的Demo：

import cv2import numpy as npimg = cv2.imread(mavic air.jpg)grayImg = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)grayImg = np.float32(grayImg)harriImg = cv2.cornerHarris(grayImg, 2, 3, 0.04)# 膨脹運算kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (3, 3))harriImgD = cv2.dilate(harriImg, kernel)# 設置閾值 實際是對Harri的輸出矩陣進行過濾 而輸出矩陣又可以用圖像形象的表現img[harriImgD > 0.05*harriImgD.max()] = [0, 0, 255]cv2.imshow(Harri, harriImgD)cv2.imshow(Output, img)cv2.waitKey()cv2.destroyAllWindows()

Shi-Tomasi演算法

Shi-Tomasi演算法是J. Shi和C. Tomasi兩位大神在1994年對Harri演算法改進後得出的一種演算法，並發表在了一篇叫Good Features To Track的論文中。

Shi-Tomasi演算法相比Harri演算法的不同點在於最後一步，即通過特徵值確定圖像特徵。前文寫到，Harri演算法要判斷特徵值的相對關係，但畢竟是相對的，用起來很糾結。而Shi-Tomasi演算法就專治各種糾結，我給一個閾值（下圖中的 ），然後直接用閾值來界定特徵值的關係，進而簡單又粗暴的確定圖像特徵。下圖中灰色為角點，橙色為邊緣，綠色為團塊。

數學表達式為：

API，沒錯，就是用這兩個大神的論文名稱作的函數名：

這個函數可以直接提供角點的坐標，而不像Harri的函數需要進一步操作才能獲取角點坐標。

corners=cv2.goodFeaturesToTrack(image, maxCorners, qualityLevel, minDistance[, corners[, mask[, blockSize[, useHarrisDetector[, k]]]]])

• corners：輸出角點向量
• image：輸入圖像，單通道，8位或32位浮點型
• maxCorners：輸出角點的數量
• qualityLevel：演算法可接受的最小特徵值，實際用於角點檢測的最小特徵值是該參數與與圖像最大特徵值的乘積，所以該參數通常不會大於1，常用值為0.10或0.01。檢測完角點後還要過濾掉距離較近的點
• minDistance：角點之間的最小距離，小於該距離的角點將被過濾掉
• mask：掩膜操作，用於對ROI參數進行處理
• blockSize：默認值為3，鄰域大小
• useHarrisDetector：是否使用Harri角點檢測
• k：默認值0.04，用於設置權重係數

import numpy as npimport cv2img = cv2.imread(mavic air.jpg)grayImg = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)# 輸出最多50個角點 qualityLevel為0.01 角點間最小距離為50cornerPoints = cv2.goodFeaturesToTrack(grayImg, 50, 0.01, 50)cornerPoints = np.int0(cornerPoints)for i in cornerPoints: x, y = i.ravel() cv2.circle(img, (x, y), 3, [0, 0, 255], -1)cv2.imshow(Shi-Tomasi, img)cv2.waitKey()

亞像素角點檢測

上述講到的兩種角點檢測演算法能較好的找到角點，但其實結果並不精確。當我們需要進行精確的角點計算時，就要用到亞像素角點檢測（顧名思義，都亞像素了，能不精確一點嗎）。亞像素角點檢測在攝像機標定、跟蹤並重建攝像機的軌跡或進行三維重建時有著重要作用。

那亞像素到底「亞」在哪兒了？通常我們計算出的坐標都是正整數，這就意味著我們是在對像素進行操作（注意，像素是圖像處理的基本單位），而亞像素計算出來的坐標是實數，形象地說，就是把像素細分成若干像素（就像原子物理學裡的夸克理論，夸克比原子還小），但由於最終還是要回歸到像素級（整數級）操作，所以還要做數據類型的轉換。

亞像素角點檢測演算法的思想是對正交向量的點乘進行迭代。我們來看下圖：

我們假設點 是起始角點，並定義向量 。

• 左圖中，點 在團塊區域中，梯度為0。
• 右圖中，點 在邊緣上，所以點 的梯度方向垂直與邊緣，向量 與點 的梯度向量的點乘積為0。
• 當點 固定時，若有足夠多的向量 滿足向量 與點 的梯度向量的點乘積為0，將每個向量看成一個方程，對方程組求解，則點 就是亞像素角點的精確位置。

API：

corners=cv2.cornerSubPix(image, corners, winSize, zeroZone, criteria)

• corners：輸出角點向量，提供精確坐標
• image：輸入圖像
• winSize：搜索窗口的一半尺寸，例如，若該參數為 ，則搜索窗口尺寸為
• zeroZone：死區的一半尺寸，死區為不對搜索窗口中央區域作求和運算的區域，用於避免矩陣的奇異性， 表示沒有死區
• criteria：求角點位置的迭代過程的結束條件，可以是最大的迭代次數，亦可是迭代結束所期望達到的精度，亦可是二者結合

criteria的定義格式

(type,maxCount,epsilon)

• type：終止迭代的方式，見下圖，python中對應為cv2.TERM_CRITERIA_EPS、cv2.TERM_CRITERIA_ITER和cv2.TERM_CRITERIA_COUNT
• maxCount：最大迭代次數
• epsilon：期望迭代後達到的精度

Demo：

這個Demo改編自官方指導，源代碼中使用Harri演算法得到角點後進行了連通域處理，目的是確定Harri角點周圍的微小團塊區域，為得到精確的亞像素角點位置做準備。

小林在其中添加了使用Shi-Tomasi演算法後直接進行亞像素角點的尋找，運行Demo可以看出，此方法相比官方指導的方法，精度略低。

import cv2import numpy as npfilename = mavic air.jpgimg = cv2.imread(filename)gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)# 使用Harri演算法和連通域找到亞像素點gray32 = np.float32(gray)dst = cv2.cornerHarris(gray, 2, 3, 0.04)dst = cv2.dilate(dst, None)ret, dst = cv2.threshold(dst, 0.01*dst.max(), 255, 0)dst = np.uint8(dst)cv2.imshow(dst, dst)# 尋找連通域矩心ret, labels, stats, centroids = cv2.connectedComponentsWithStats(dst)# 確定亞像素角點檢測的迭代條件criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 100, 0.001)corners = cv2.cornerSubPix(gray32, np.float32(centroids), (5, 5), (-1, -1), criteria)# 使用Shi-Tomasi演算法找到角點cornerPoints = cv2.goodFeaturesToTrack(gray, 50, 0.01, 50)cornersS = cv2.cornerSubPix(gray, cornerPoints, (5, 5), (-1, -1), criteria)# 橫向堆疊矩心和角點的坐標矩陣res = np.hstack((centroids, corners))res = np.int0(res)print(res)# 藍色點表示Harri演算法後找到的角點# 綠色點表示Harri演算法後和連通域找到的亞像素角點# 紅色點表示Shi-Tomasi演算法後找到的亞像素角點for i in cornersS: x, y = i.ravel() cv2.circle(img, (x, y), 4, [0, 0, 255])for i in range(0, res[:, 3].size): cv2.circle(img, (res[i, 0], res[i, 1]), 4, [255, 0, 0]) cv2.circle(img, (res[i, 2], res[i, 3]), 3, [0, 255, 0])cv2.imshow(Image, img)cv2.waitKey()cv2.destroyAllWindows()

這一話同學們看起來可能會有些頭疼，那就讓Menhera醬為你加油吧

這一話的Demo已經同步到Github了，對應Class 5 Feature Detection and Description文件夾下的C5 Harri.py、C5 ShiTomasi.py和C5 Subpixel.py。