harris角點檢測的簡要總結

時間 2019-12-09

標籤 harris 檢測簡要總結简体版

原文原文鏈接

目錄html

1. 概述相關

harris角點檢測是一種特徵提取的方法，而特徵提取正是計算機視覺的一種重要手段。儘管它看起來很複雜，其實也是基於數學原理和簡單的圖像處理來實現的。
本文以前能夠參看筆者寫的幾篇圖像處理的文章，將會有助於更深刻了解harris角點檢測的實現。ios

2. 原理詳解

1) 算法思想

爲了判斷圖像的角點，能夠利用卷積窗口滑動的思想，讓以該點爲中心的窗口在附近滑動。以下圖是全部描述角點文章的初始圖例，它表徵的正是這一特性：當滑動窗口在全部方向移動時，窗口內的像素灰度出現了較大的變化，就多是角點。
算法

2) 數學模型

根據上述的算法思想，能夠構建數學模型，圖像窗口平移[u,v]產生灰度變化E(u,v)爲：

其中w(x,y)是一種加權函數，幾乎全部的應用都把它設爲高斯函數。由上述公式，進行推導以下：

最後獲得的公式(6)，在幾何意義上表徵的是一個橢圓。橢圓的長短軸分別沿着矩陣M的兩個特徵向量的方向，而兩個與之對應的特徵值分別是半長軸和半短軸的長度的平方的倒數。

那麼根據矩陣M的兩個特徵值λ1和λ2，能夠將圖像上的像素點分類成直線、平面與角點：當λ1和λ2 都比較大，且近似相等時，能夠認爲是角點。以下圖所示：
函數

3) 優化推導

而上述表達不太方便使用，又定義了一個角點響應函數R，經過R的大小來判斷像素是否爲角點：

式中，detM爲矩陣M的行列式，traceM爲矩陣M的直跡。α爲常常常數，取值範圍爲0.04~0.06。對於R公式，有推導以下：

能夠知道，角點響應值R仍然表徵了矩陣M兩個特徵值λ1和λ2，一樣能夠進行上述分類：當R爲大數值正數的時候，表示爲角點。以下圖所示：
優化

3. 具體實現

在OpenCV中，已經提供了Harris角點檢測函數cornerHarris()。爲了更好地理解Harris角點提取的原理，這裏參考了網上代碼，本身實現了其算法，不過也調用了OpenCV中一些基本函數。
根據上述原理，Harris圖像角點檢測算法的關鍵是計算M矩陣，M矩陣是圖像I(x,y)的偏導數矩陣，也就是要先求出圖像的梯度。spa

1) 詳細步驟

1.計算圖像I(x,y)在X，Y方向的梯度。在這裏是經過卷積函數filter2D實現的，具體原理能夠看(1)中提到的相關文章。3d

Mat gray;
imgSrc.convertTo(gray, CV_64F);

Mat xKernel = (Mat_<double>(1, 3) << -1, 0, 1);
Mat yKernel = xKernel.t();

Mat Ix, Iy;
filter2D(gray, Ix, CV_64F, xKernel);
filter2D(gray, Iy, CV_64F, yKernel);

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

Mat Ix2, Iy2, Ixy;
Ix2 = Ix.mul(Ix);
Iy2 = Iy.mul(Iy);
Ixy = Ix.mul(Iy);

3.對Ix二、Iy2和Ixy進行高斯濾波，生成矩陣M的元素A、B和C。htm

Mat gaussKernel = getGaussianKernel(7, 1);
filter2D(Ix2, Ix2, CV_64F, gaussKernel);
filter2D(Iy2, Iy2, CV_64F, gaussKernel);
filter2D(Ixy, Ixy, CV_64F, gaussKernel);

4.根據公式計算每一個像素的Harris響應值R，獲得圖像對應的響應值矩陣。blog

Mat cornerStrength(gray.size(), gray.type());
for (int i = 0; i < gray.rows; i++)
{
    for (int j = 0; j < gray.cols; j++)
    {
        double det_m = Ix2.at<double>(i, j) * Iy2.at<double>(i, j) - Ixy.at<double>(i, j) * Ixy.at<double>(i, j);
        double trace_m = Ix2.at<double>(i, j) + Iy2.at<double>(i, j);
        cornerStrength.at<double>(i, j) = det_m - alpha * trace_m *trace_m;
    }
}

5.在3×3的鄰域內進行非最大值抑制，找到局部最大值點，即爲圖像中的角點。在這裏非最大值抑制是經過圖像膨脹的實現的。比較膨脹先後的響應值矩陣，獲得局部最大值。

//在3×3的鄰域內進行非最大值抑制，找到局部最大值點，即爲圖像中的角點
double maxStrength;
minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);
Mat dilated;
Mat localMax;
dilate(cornerStrength, dilated, Mat());             //膨脹
compare(cornerStrength, dilated, localMax, CMP_EQ);      //比較保留最大值的點

//獲得角點的位置
Mat cornerMap;
double qualityLevel = 0.01;
double thresh = qualityLevel * maxStrength;   
cornerMap = cornerStrength > thresh;                //小於閾值t的R置爲零。
bitwise_and(cornerMap, localMax, cornerMap);            //位與運算，有0則爲0, 全爲1則爲1

imgDst = cornerMap.clone();

2) 最終實現

合併以上步驟，傳入參數，最終的實現代碼：

#include <iostream>
#include <algorithm>
#include <opencv2\opencv.hpp>

using namespace cv;
using namespace std;

void detectHarrisCorners(const Mat& imgSrc, Mat& imgDst, double alpha)
{
    //
    Mat gray;
    imgSrc.convertTo(gray, CV_64F);

    //計算圖像I(x,y)在X，Y方向的梯度
    Mat xKernel = (Mat_<double>(1, 3) << -1, 0, 1);
    Mat yKernel = xKernel.t();

    Mat Ix, Iy;
    filter2D(gray, Ix, CV_64F, xKernel);
    filter2D(gray, Iy, CV_64F, yKernel);

    //計算圖像兩個方向梯度的乘積。
    Mat Ix2, Iy2, Ixy;
    Ix2 = Ix.mul(Ix);
    Iy2 = Iy.mul(Iy);
    Ixy = Ix.mul(Iy);

    //對Ix二、Iy2和Ixy進行高斯濾波，生成矩陣M的元素A、B和C。
    Mat gaussKernel = getGaussianKernel(7, 1);
    filter2D(Ix2, Ix2, CV_64F, gaussKernel);
    filter2D(Iy2, Iy2, CV_64F, gaussKernel);
    filter2D(Ixy, Ixy, CV_64F, gaussKernel);

    //根據公式計算每一個像素的Harris響應值R，獲得圖像對應的響應值矩陣。
    Mat cornerStrength(gray.size(), gray.type());
    for (int i = 0; i < gray.rows; i++)
    {
        for (int j = 0; j < gray.cols; j++)
        {
            double det_m = Ix2.at<double>(i, j) * Iy2.at<double>(i, j) - Ixy.at<double>(i, j) * Ixy.at<double>(i, j);
            double trace_m = Ix2.at<double>(i, j) + Iy2.at<double>(i, j);
            cornerStrength.at<double>(i, j) = det_m - alpha * trace_m *trace_m;
        }
    }

    //在3×3的鄰域內進行非最大值抑制，找到局部最大值點，即爲圖像中的角點
    double maxStrength;
    minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);
    Mat dilated;
    Mat localMax;
    dilate(cornerStrength, dilated, Mat());             //膨脹
    compare(cornerStrength, dilated, localMax, CMP_EQ);      //比較保留最大值的點
    
    //獲得角點的位置
    Mat cornerMap;
    double qualityLevel = 0.01;
    double thresh = qualityLevel * maxStrength;   
    cornerMap = cornerStrength > thresh;                //小於閾值t的R置爲零。
    bitwise_and(cornerMap, localMax, cornerMap);            //位與運算，有0則爲0, 全爲1則爲1

    imgDst = cornerMap.clone();
}

//在角點位置繪製標記
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);
    }
}

int main()
{
    //從文件中讀取成灰度圖像
    const char* imagename = "D:\\Data\\imgDemo\\whdx.jpg";
    Mat img = imread(imagename, IMREAD_GRAYSCALE);
    if (img.empty())
    {
        fprintf(stderr, "Can not load image %s\n", imagename);
        return -1;
    }

    //
    Mat imgDst;
    double alpha = 0.05;    
    detectHarrisCorners(img, imgDst, alpha);
    
    //在角點位置繪製標記
    drawCornerOnImage(img, imgDst);

    //
    imshow("Harris角點檢測", img);
    waitKey();
    
    return 0;
}

其運行結果爲：