Harris角点检测实现详解

本篇文章是对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;}

以下为测试结果：