本篇文章是對Harris角點檢測基於OpenCV 3.2.0的C++實現總結與記錄，程式程式碼力圖實現最基本的演算法步驟而無關優化，目的在於理解演算法，關於Harris角點的基本原理，本篇不予贅述，網上搜索即可，下面只列出出自Computer Vision:Algorithm and Application的演算法步驟：
1, Compute the horizontal and vertical derivatives of the image Ix and Iy by convolving the original image with derivatives of Gaussians.
2, Compute the three images corresponding to the outer products of these gradients.
3, Convolve each of these images with a larger Gaussian.
4, Compute a scalar interest measure using one of the formulas discussed above.
5, Find local maxima above a certain threshold and report them as detected feature point locations.

首先觀察這5個實現步驟，發現我們會需要x方向的高斯核、y方向的高斯核與原始高斯核，於是先實現這幾個高斯核再說。

/***
 * 引數ksize: 高斯核大小，若ksize==5，則高斯核為5*5大小
 * 引數sigma: 高斯標準差
 * 返回值：x方向上的高斯核
 * 注1：二維高斯一階導數公式為：G`(x, sigma) = -x/(sigma*sigma)*G(x, y, sigma)
 * 注2：x方向上的座標計算(5*5)為[-2 -1 0 1 2; -2 -1 0 1 2; -2 -1 0 1 2; -2 -1 0 1 2; -2 -1 0 1 2]
 */
Mat  myGetGaussianKernel_x(int
 
 ksize, float sigma)
{
    Mat kernel = myGetGaussianKernel(ksize, sigma); // 首先得到原始高斯核
    Mat kernel_x(ksize, ksize, CV_32F, Scalar(0.0)); // 定義一階高斯核
    for (int x = -ksize/2; x <= ksize/2; ++x) // 若為5*5大小，則x = (-2:1:2)
    {
        for (int i = 0; i < ksize; ++i)
        {
            kernel_x.at<float
 
>(i, x + ksize/2) = -x/(sigma * sigma) * kernel.at<float>(i, x + ksize/2);
        }
    }
    return kernel_x;
}
/***
 * 引數ksize: 高斯核大小，若ksize==5，則高斯核為5*5大小
 * 引數sigma: 高斯標準差
 * 返回值：y方向上的高斯核
 * 注1：二維高斯一階導數公式為：G`(y, sigma) = -y/(sigma*sigma)*G(x, y, sigma)
 * 注2：y方向上的座標計算(5*5)為[-2 -2 -2 -2 -2; -1 -1 -1 -1 -1; 0 0 0 0 0; 1 1 1 1 1; 2 2 2 2 2]
 */
Mat  myGetGaussianKernel_y(int ksize, float sigma)
{
    Mat kernel = myGetGaussianKernel(ksize, sigma); // 首先得到原始高斯核
    Mat kernel_y(ksize, ksize, CV_32F, Scalar(0.0)); // 定義一階高斯核
    for (int y = -ksize/2; y <= ksize/2; ++y) // 若為5*5大小，則y = (-2:1:2)
    {
        for (int i = 0; i < ksize; ++i)
        {
            kernel_y.at<float>(y + ksize/2, i) = -y/(sigma * sigma) * kernel.at<float>(y + ksize/2, i);
        }
    }
    return kernel_y;
}

/***
 * 引數ksize: 高斯核大小，若ksize==5，則高斯核為5*5大小
 * 引數sigma: 高斯標準差
 * 返回值：高斯核
 */
Mat  myGetGaussianKernel(int ksize, float sigma)
{
    if (ksize % 2 == 0)
    {
        cerr << "invalid kernel size." << endl;
        return Mat(1, 1, CV_32F, Scalar(-1));
    }
    Mat kernel(ksize, ksize, CV_32F, Scalar(0.0));
    Mat kernel_1d(ksize, 1, CV_32F, Scalar(0.0));
    for (int x = -ksize/2; x <= ksize/2; ++x)
    {
        kernel_1d.at<float>(x + ksize/2, 0) = exp(-(x * x)/(2 * sigma * sigma)) / (sigma * sqrt(2 * CV_PI));
    }
    kernel = kernel_1d * kernel_1d.t(); // 這裡用兩個一維相乘得到
    return kernel;
}

各種高斯核得到之後，就開始進行Harris角點檢測的步驟了，下面上程式碼以及註釋：

/***
 *引數img：要進行角點檢測的輸入
 *引數featurePts：儲存角點位置
 */
void HarrisCornerDetector(Mat img, vector<Point> &featurePts)
{
    img.convertTo(img, CV_64FC1);
    Mat img_x, img_y;
    // 1, Compute the horizontal and vertical derivatives of the image Ix and Iy by convolving the original image with derivatives of Gaussians.
    Mat gau_x = myGetGaussianKernel_x(5, 1.0);
    Mat gau_y = myGetGaussianKernel_y(5, 1.0);
    filter2D(img, img_x, img.depth(), gau_x);
    filter2D(img, img_y, img.depth(), gau_y);

    // 2, Compute the three images corresponding to the outer products of these gradients.
    // 3, Convolve each of these images with a larger Gaussian.
    Mat img_x2, img_y2, img_xy;
    Mat gau = myGetGaussianKernel(9, 2.0); // 此處的高斯核要larger，所以我定義了9*9, sigma = 2.0的高斯核，larger的影響是這兩個引數的增大使角點更為稀疏
    filter2D(img_x.mul(img_x), img_x2, img_x.depth(), gau);
    filter2D(img_y.mul(img_y), img_y2, img_y.depth(), gau);
    filter2D(img_x.mul(img_y), img_xy, img_x.depth(), gau);

    //4, Compute a scalar interest measure using one of the formulas discussed above.
    Mat cim(img.rows, img.cols, CV_64FC1, Scalar(0.0));
    for (int i = 0; i < cim.rows; ++i)
    {
        for (int j = 0; j < cim.cols; ++j)
        {
            // 這個公式採用了det(A)/tr(A)
            cim.at<double>(i, j) = (img_x2.at<double>(i, j) * img_y2.at<double>(i, j)
                                    - img_xy.at<double>(i, j) * img_xy.at<double>(i, j))
                                    /(img_x2.at<double>(i, j) + img_y2.at<double>(i, j));
        }
    }
    // 5, Find local maxima above a certain threshold and report them as detected feature point locations.
    for (int i = 1; i < cim.rows - 1; ++i)
    {
        for (int j = 1; j < cim.cols - 1; ++j)
        {
            // 下面的for迴圈就是進行3*3視窗的非最大值抑制以及閾值處理
            double lmax = cim.at<double>(i, j);
            for (int x = -1; x <= 1; ++x)
            {
                for (int y = -1; y <= 1; ++y)
                {
                    if (cim.at<double>(i + x, j + y) > lmax)
                    {
                        lmax = cim.at<double>(i + x, j + y);
                    }
                }
            }
            double thresh = 30.0;
            // 選擇區域性最大值且最大值大於閾值的點作為角點
            if (cim.at<double>(i, j) == lmax && lmax > thresh)
            {
                // 這裡要注意一下i和j互換一下，這裡的i表示行數，在畫角點時代表的是OpenCV座標系的y座標，相應的j代表的是x座標
                featurePts.push_back(Point(j, i));
            }
        }
    }
}

下面進行測試：

#include <iostream>
#include <vector>
#include <cmath>
#include <opencv2/opencv.hpp>

using namespace std;
using namespace cv;

Mat  myGetGaussianKernel(int ksize = 5, float sigma = 1.0);
Mat  myGetGaussianKernel_x(int ksize = 5, float sigma = 1.0);
Mat  myGetGaussianKernel_y(int ksize = 5, float sigma = 1.0);
void HarrisCornerDetector(Mat img, vector<Point> &featurePts);

int main()
{
    Mat img = imread("lena.tiff", IMREAD_GRAYSCALE);
    vector<Point> featurePts;
    HarrisCornerDetector(img, featurePts);
    for (int i = 0; i < featurePts.size(); ++i)
    {
        circle(img, featurePts[i], 5, Scalar(255));
    }
    imshow("Harris Corner", img);
    waitKey(0);
    return 0;
}

以下為測試結果：