Harris角點

時間 2021-03-13

標籤 html git github 算法 app 函數學習編碼 url es5 欄目 HTML 简体版

原文原文鏈接

1. 不一樣類型的角點

在現實世界中，角點對應於物體的拐角，道路的十字路口、丁字路口等。從圖像分析的角度來定義角點能夠有如下兩種定義：html

角點能夠是兩個邊緣的角點；
角點是鄰域內具備兩個主方向的特徵點；

前者每每須要對圖像邊緣進行編碼，這在很大程度上依賴於圖像的分割與邊緣提取，具備至關大的難度和計算量，且一旦待檢測目標局部發生變化，極可能致使操做的失敗。早期主要有Rosenfeld和Freeman等人的方法，後期有CSS等方法。git

基於圖像灰度的方法經過計算點的曲率及梯度來檢測角點，避免了第一類方法存在的缺陷，此類方法主要有Moravec算子、Forstner算子、Harris算子、SUSAN算子等。github

這篇文章主要介紹的Harris角點檢測的算法原理，比較著名的角點檢測方法還有jianbo Shi和Carlo Tomasi提出的Shi-Tomasi算法，這個算法開始主要是爲了解決跟蹤問題，用來衡量兩幅圖像的類似度，咱們也能夠把它看爲Harris算法的改進。OpenCV中已經對它進行了實現，接口函數名爲GoodFeaturesToTrack()。另外還有一個著名的角點檢測算子即SUSAN算子，SUSAN是Smallest Univalue Segment Assimilating Nucleus（最小核值類似區）的縮寫。SUSAN使用一個圓形模板和一個圓的中心點，經過圓中心點像素與模板圓內其餘像素值的比較，統計出與圓中心像素近似的像元數量，當這樣的像元數量小於某一個閾值時，就被認爲是要檢測的角點。我以爲能夠把SUSAN算子看爲Harris算法的一個簡化。這個算法原理很是簡單，算法效率也高，因此在OpenCV中，它的接口函數名稱爲：FAST() 。算法

2. Harris角點

2.1 基本原理

人眼對角點的識別一般是在一個局部的小區域或小窗口完成的。若是在各個方向上移動這個特徵的小窗口，窗口內區域的灰度發生了較大的變化，那麼就認爲在窗口內遇到了角點。若是這個特定的窗口在圖像各個方向上移動時，窗口內圖像的灰度沒有發生變化，那麼窗口內就不存在角點；若是窗口在某一個方向移動時，窗口內圖像的灰度發生了較大的變化，而在另外一些方向上沒有發生變化，那麼，窗口內的圖像可能就是一條直線的線段。app

對於圖像$I(x,y)$，當在點$(x,y)$處平移$(\Delta x,\Delta y)$後的自類似性，能夠經過自相關函數給出：函數

$$c(x,y;\Delta x,\Delta y) = \sum_{(u,v)\in W(x,y)}w(u,v)(I(u,v) – I(u+\Delta x,v+\Delta y))^2$$學習

其中，$W(x,y)$是以點$(x,y)$爲中心的窗口，$w(u,v)$爲加權函數，它既但是常數，也能夠是高斯加權函數。編碼

根據泰勒展開，對圖像$I(x,y)$在平移$(\Delta x,\Delta y)$後進行一階近似：url

$$I(u+\Delta x,v+\Delta y) = I(u,v)+I_x(u,v)\Delta x+I_y(u,v)\Delta y+O(\Delta x^2,\Delta y^2)\approx I(u,v)+I_x(u,v)\Delta x+I_y(u,v)\Delta y$$es5

其中，$I_x,I_y$是圖像$I(x,y)$的偏導數，這樣的話，自相關函數則能夠簡化爲：

$$c(x,y;\Delta x,\Delta y)\approx \sum_w (I_x(u,v)\Delta x+I_y(u,v)\Delta y)^2=[\Delta x,\Delta y]M(x,y)\begin{bmatrix}\Delta x \\ \Delta y\end{bmatrix}$$

其中

$$M(x,y)=\sum_w \begin{bmatrix}I_x(x,y)^2&I_x(x,y)I_y(x,y) \\ I_x(x,y)I_y(x,y)&I_y(x,y)^2\end{bmatrix} = \begin{bmatrix}\sum_w I_x(x,y)^2&\sum_w I_x(x,y)I_y(x,y) \\\sum_w I_x(x,y)I_y(x,y)&\sum_w I_y(x,y)^2\end{bmatrix}=\begin{bmatrix}A&C\\C&B\end{bmatrix} $$

也就是說圖像$I(x,y)$在點$(x,y)$處平移$(\Delta x,\Delta y)$後的自相關函數能夠近似爲二項函數：

$$c(x,y;\Delta x,\Delta y)\approx A\Delta x^2+2C\Delta x\Delta y+B\Delta y^2$$

其中

$$A=\sum_w I_x^2, B=\sum_w I_y^2,C=\sum_w I_x I_y$$

二次項函數本質上就是一個橢圓函數。橢圓的扁率和尺寸是由$M(x,y)$的特徵值$\lambda_一、\lambda_2$決定的，橢賀的方向是由$M(x,y)$的特徵矢量決定的，以下圖所示，橢圓方程爲：

$$[\Delta x,\Delta y]M(x,y)\begin{bmatrix}\Delta x \\ \Delta y\end{bmatrix} = 1$$

橢圓函數特徵值與圖像中的角點、直線（邊緣）和平面之間的關係以下圖所示。共可分爲三種狀況：

圖像中的直線。一個特徵值大，另外一個特徵值小，$\lambda_1\gg \lambda_2$或$\lambda_2\gg \lambda_1$。自相關函數值在某一方向上大，在其餘方向上小。
圖像中的平面。兩個特徵值都小，且近似相等；自相關函數數值在各個方向上都小。
圖像中的角點。兩個特徵值都大，且近似相等，自相關函數在全部方向都增大。

根據二次項函數特徵值的計算公式，咱們能夠求$M(x,y)$矩陣的特徵值。可是Harris給出的角點差異方法並不須要計算具體的特徵值，而是計算一個角點響應值$R$來判斷角點。$R$的計算公式爲：

$$R=det \boldsymbol{M} - \alpha(trace\boldsymbol{M})^2$$

式中，$det\boldsymbol{M}$爲矩陣$\boldsymbol{M}=\begin{bmatrix}A&B\\B&C\end{bmatrix}$的行列式；$trace\boldsymbol{M}$爲矩陣$\boldsymbol{M}$的直跡；$\alpha$爲常常常數，取值範圍爲0.04~0.06。事實上，特徵是隱含在$det\boldsymbol{M}$和$trace\boldsymbol{M}$中，由於，

$$det\boldsymbol{M} = \lambda_1\lambda_2=AC-B^2$$

$$trace\boldsymbol{M}=\lambda_2+\lambda_2=A+C$$

2.2 Harris角點算法實現

根據上述討論，能夠將Harris圖像角點檢測算法概括以下，共分如下五步：

1. 計算圖像$I(x,y)$在$X$和$Y$兩個方向的梯度$I_x、I_y$。

$$I_x=\frac{\partial I}{\partial x}=I\otimes(-1\ 0\ 1)，I_y =\frac{\partial I}{\partial x}=I\otimes(-1\ 0\ 1)^T $$

2. 計算圖像兩個方向梯度的乘積。

$$I_x^2=I_x\cdot I_y，I_y^2=I_y\cdot I_y，I_{xy}=I_x\cdot I_y$$

3. 使用高斯函數對$I_x^二、I_y^2和I_{xy}$進行高斯加權（取$\sigma=1$），生成矩陣$M$的元素$A、B$和$C$。

$$A=g(I_x^2)=I_x^2\otimes w，C=g(I_y^2)=I_y^2\otimes w，B=g(I_{x,y})=I_{xy}\otimes w$$

4. 計算每一個像素的Harris響應值$R$，並對小於某一閾值$t$的$R$置爲零。

$$R=\{R:det \boldsymbol{M} - \alpha(trace\boldsymbol{M})^2<t\}$$

5. 在$3\times3$或$5\times5$的鄰域內進行非最大值抑制，局部最大值點即爲圖像中的角點。

Harris角點檢測的C++實現代碼：https://github.com/RonnyYoung/ImageFeatures/blob/master/source/harris.cpp

2.3 Harris角點的性質

1. 參數$\alpha$對角點檢測的影響

假設已經獲得了矩陣$\boldsymbol{M}$的特徵值$\lambda_1\ge\lambda_2\ge0$，令$\lambda_2=k\lambda_1,0\le k\le 1$。由特徵值與矩陣$\boldsymbol{M}$的直跡和行列式的關係可得：

$$det\boldsymbol{M}=\prod_i\lambda_i \ \ \ \ \ \ trace\boldsymbol{M}=\sum_i\lambda_i$$

從而能夠獲得角點的響應

$$R=\lambda_2\lambda_2=\alpha(\lambda_2+\lambda_2)^2=\lambda^2(k-\alpha(1+k)^2)$$

假設$R\ge0$，則有：

$$0\le \alpha \le\frac{k}{(1+k)^2}\le0.25$$

對於較小的$k$值，$R\approx\lambda^2(k-\alpha),\alpha<k$。

由此，能夠得出這樣的結論：增大$\alpha$的值，將減少角點響應值$R$，下降角點檢測的靈性，減小被檢測角點的數量；減少$\alpha$值，將增大角點響應值$R$，增長角點檢測的靈敏性，增長被檢測角點的數量。

2. Harris角點檢測算子對亮度和對比度的變化不敏感

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

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

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

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

以下圖所示，當右圖被縮小時，在檢測窗口尺寸不變的前提下，在窗口內所包含圖像的內容是徹底不一樣的。左側的圖像可能被檢測爲邊緣或曲線，而右側的圖像則可能被檢測爲一個角點。

2.4 Harris的OpenCV接口

OpenCV的Hairrs角點檢測的函數爲cornerHairrs()，可是它的輸出是一幅浮點值圖像，浮點值越高，代表越多是特徵角點，咱們須要對圖像進行閾值化。

C++: void cornerHarris(InputArray src, OutputArray dst, int blockSize, int apertureSize, double k, int borderType = BORDER_DEFAULT);

src – 輸入的單通道8-bit或浮點圖像。
dst – 存儲着Harris角點響應的圖像矩陣，大小與輸入圖像大小相同，是一個浮點型矩陣。
blockSize – 鄰域大小。
apertureSize – 擴展的微分算子大。
k – 響應公式中的，參數$\alpha$。
boderType – 邊界處理的類型。

int main()
{
    Mat image = imread("../buliding.png");
    Mat gray;
    cvtColor(image, gray, CV_BGR2GRAY);

    Mat cornerStrength;
    cornerHarris(gray, cornerStrength, 3, 3, 0.01);
    threshold(cornerStrength, cornerStrength, 0.001, 255, THRESH_BINARY);
    return 0;
}

從上面上間一幅圖像咱們能夠看到，有不少角點都是粘連在一塊兒的，咱們下面經過加入非極大值抑制來進一步去除一些粘在一塊兒的角點。

非極大值抑制原理是，在一個窗口內，若是有多個角點則用值最大的那個角點，其餘的角點都刪除，窗口大小這裏咱們用3*3，程序中經過圖像的膨脹運算來達到檢測極大值的目的，由於默認參數的膨脹運算就是用窗口內的最大值替代當前的灰度值。

int main()
{
    Mat image = imread("buliding.png");
    Mat gray;
    cvtColor(image, gray, CV_BGR2GRAY);

    Mat cornerStrength;
    cornerHarris(gray, cornerStrength, 3, 3, 0.01);

    double maxStrength;
    double minStrength;
    // 找到圖像中的最大、最小值
    minMaxLoc(cornerStrength, &minStrength, &maxStrength);

    Mat dilated;
    Mat locaMax;
    // 膨脹圖像，最找出圖像中所有的局部最大值點
    dilate(cornerStrength, dilated, Mat());
    // compare是一個邏輯比較函數，返回兩幅圖像中對應點相同的二值圖像
    compare(cornerStrength, dilated, locaMax, CMP_EQ);

    Mat cornerMap;
    double qualityLevel = 0.01;
    double th = qualityLevel*maxStrength; // 閾值計算
    threshold(cornerStrength, cornerMap, th, 255, THRESH_BINARY);
    cornerMap.convertTo(cornerMap, CV_8U);
    // 逐點的位運算
    bitwise_and(cornerMap, locaMax, cornerMap);

    drawCornerOnImage(image, cornerMap);
    namedWindow("result");
    imshow("result", image);
    waitKey();

    return 0;
}
void drawCornerOnImage(Mat& image, const Mat&binary)
{
    Mat_<uchar>::const_iterator it = binary.begin<uchar>();
    Mat_<uchar>::const_iterator itd = binary.end<uchar>();
    for (int i = 0; it != itd; it++, i++)
    {
        if (*it)
            circle(image, Point(i%image.cols, i / image.cols), 3, Scalar(0, 255, 0), 1);
    }
}

如今咱們獲得的效果就比默認的函數獲得的結果有至關的改善，如上面最右邊效果圖。

3. 多尺度Harris角點

3.1 多尺度Harris角點的原理

雖然Harris角點檢測算子具備部分圖像灰度變化的不變性和旋轉不變性，但它不具備尺度不變性。可是尺度不變性對圖像特徵來講相當重要。人們在使用肉眼識別物體時，無論物體遠近，尺寸的變化都能認識物體，這是由於人的眼睛在辨識物體時具備較強的尺度不變性。在圖像特徵提取：尺度空間理論這篇文章裏就已經講到了高斯尺度空間的概念。下面將Harris角點檢測算子與高斯尺度空間表示相結合，使用Harris角點檢測算子具備尺度的不變性。

仿照Harris角點檢測中二階矩的表示方法，使用$M=\mu(x,\sigma_I,\sigma_D)$爲尺度自適應的二階矩：

$$\boldsymbol{M}=\mu(x,\sigma_I,\sigma_D)=\sigma_D^2g(\sigma_I)\otimes\begin{bmatrix}L_x^2(x,\sigma_D)&L_xL_y(x,\sigma_D)\\L_xL_y(x,\sigma_D)&L_y^2(x,\sigma_D)\end{bmatrix} $$

其中，$g(\sigma_I)$表示尺度爲$sigma_I$的高斯卷積核，$x$表示圖像的位置。與高斯測度空間相似，使用$L(x)$表示通過高斯平滑後的圖像，符號$\otimes$表示卷積，$L_x(x,\sigma_D)$和$L_y(x,\sigma_D)$表示對圖像使用高斯$g(\sigma_D)$函數進行平滑後，在$x$或$y$方向取其微分的結果，即$L_x=\partial_xL$和$L_y=\partial_yL$。一般將$\sigma_I$稱爲積分尺度，它是決定Harris角點當前尺度的變量，$\sigma_D$爲微分尺度或局部尺度，它是決定角點附近微分值變化的變量。顯然，積分尺度$\sigma_I$應該大於微分尺度$\sigma_D$。

3.2 多尺度Harris角點實現

首先，檢測算法從預先定義的一組尺度中進行積分尺度搜索，這一組尺度定義爲

$$\sigma_1\dots\sigma_n = \sigma_0\dots k^n\sigma_0$$

通常狀況下使用k=1.4。爲了減小搜索的複雜性，對於微分尺度$\sigma_D$的選擇，咱們採用在積分尺度的基礎上，乘以一個比例常數，即$\sigma_D=s\sigma_I$，通常取s=0.7。這樣，一般使用積分和微分的尺度，即可以生成$\mu(x,\sigma_I,\sigma_D)$，再利用Harris角點判斷準則，對角點進行搜索，具體能夠分兩步進行。

1. 與Harris角點搜索相似，對於給定的尺度空間值$\sigma_D$，進行以下角點響應值計算和判斷：

$$cornerness = det(\mu(x,\sigma_n)-\alpha trace^2(\mu(x,\sigma_n)))>threshold_H$$

2. 對於知足1中條件的點，在點的8鄰域內進行角點響應最大值搜索（即非最大值抑制）出在8鄰域內角點響應最大值的點。對於每一個尺度$\sigma_n(1,2,\dots,n)$都進行如上搜索。

因爲位置空間的候選點並不必定在尺度空間上也能成爲候選點，因此，咱們還要在尺度空間上進行搜索，找到該點的所謂特徵尺度值。搜索特徵尺度值也分兩步。

1. 對於位置空間搜索到的每一個候選點，進行拉普拉斯響應計算，並知足其絕對值大於給定的閾值條件：

$$F(x,\sigma_n) = \sigma_n^2|L_{xx}(x,\sigma_n)+L_{yy}(x,\sigma_n)| \ge threshold_L$$

2. 與鄰近的兩個尺度空間的拉普拉斯響應值進行比較，使其知足：

$$F(x,\sigma_n)>F(x,\sigma_l), \ \ \ l\in \{n-1.n+1\}$$

知足上述條件的尺度值就是該點的特徵尺度值。這樣，咱們就找到了在位置空間和尺度空間都知足條件的Harris角點。