Harris角点检测原理及C++实现

（1）角点检测的核心思想：
　　使用一个固定窗口在图像上进行任意方向上的滑动，比较滑动前与滑动后两种情况，窗口中的像素灰度变化程度，如果存在任意方向上的滑动，都有着较大灰度变化，那么我们可以认为该窗口中存在角点。
（2）灰度变化描述
　　当窗口发生[u,v]移动时，那么滑动前与滑动后对应的窗口中的像素点灰度变化描述如下：
　　E(u,v)=∑(x,y)€Ww(x,y)[I(x+u,y+v)−I(x,y)]2
参数解释：

[u,v]是窗口W的偏移量；

(x,y)是窗口W所对应的像素坐标位置，窗口有多大，就有多少个位置；
I(x,y)是像素坐标位置(x,y)的图像灰度值；
I(x+u,y+v)是像素坐标位置(x+u,y+v)的图像灰度值；
w(x,y)是窗口函数，最简单情形就是窗口W内的所有像素所对应的w权重系数均为1.但有时候，我们会将w(x,y)函数设置为以窗口W中心为原点的二元正太分布。如果窗口W中心点是角点时，移动前与移动后，
该点在灰度变化贡献最大；而离窗口W中心(角点)较远的点，这些点的灰度变化几近平缓，这些点的权重系数，可以设定小值，以示该点对灰度变化贡献较小，那么我们自然而然想到使用二元高斯函数来表示窗口函数；

上述是Moravec角点检测算法，该算法是基于相邻像素之间的欧氏距离度量灰度变化量，缺点：不具有旋转不变性（同一张图片旋转一定角度后，检测到的角点不同），给实际工程应用带来不便。

针对上述问题：Harris运用泰勒公式做出改进，说不明白，上图如下：

至于算法的推导，可参考：https://www.cnblogs.com/zyly/p/9508131.html#_label5_0；

算法实现：

main.cpp

#include "iostream"
#include "opencv2/core.hpp"
#include "opencv2/highgui.hpp"
#include "my_harris.h"
using namespace cv;
using namespace std;

int main()
{
    my_Harris Harris;
    Mat img = imread("D:\qtproject\img\4.jpg");
    //imshow("img",img);
    //waitKey(0);
    /*0.图像预处理之灰度化RGB2GRAY*/
    int cols = img.cols;
    int rows =  img.rows;
    Mat img_gray = Mat(rows,cols,CV_8UC1,Scalar(0));
    Harris.RGB2GRAY(img,img_gray);
    /*1.计算图像在x,y方向的梯度*/
    Mat gx_img = Mat(rows,cols,CV_8UC1,Scalar(0));
    Mat gy_img = Mat(rows,cols,CV_8UC1,Scalar(0));
    Harris.IMG_GRAD(img_gray,gx_img,gy_img);

    /*2.计算图像在两个方向梯度的乘积*/
    Mat gx_p = Mat(rows,cols,CV_8UC1,Scalar(0));
    Mat gy_p = Mat(rows,cols,CV_8UC1,Scalar(0));
    Mat gxy = Mat(rows,cols,CV_8UC1,Scalar(0));
    Harris.GRAD_MULTI(gx_img,gy_img,gx_p,gy_p,gxy);

    /*3.使用高斯函数加权*/
    Mat A = Mat(rows,cols,CV_8UC1,Scalar(0));
    Mat B = Mat(rows,cols,CV_8UC1,Scalar(0));
    Mat C = Mat(rows,cols,CV_8UC1,Scalar(0));
    Harris.GAUSS_WEI(gx_p,gy_p,gxy,A,B,C);
    /*4.计算每个像素点的harris响应值*/
    Mat R = Mat(rows,cols,CV_8UC1,Scalar(0));
    Harris.GET_RESPONSE(A,B,C,R);
    /*5.过滤大于某一阈值t的R值*/
    Mat CORNER = Mat(rows,cols,CV_8UC1,Scalar(0));
    Harris.FILTER_THRESH(R,img,CORNER);
    return 0;
}

my_harris.cpp

#include "my_harris.h"
#define pi 3.14
my_Harris::my_Harris()
{
}
void my_Harris::RGB2GRAY(Mat rgb_img, Mat &gray_img)
{
    Mat img_src = rgb_img.clone();
    int rows =  img_src.rows;
    int cols = img_src.cols;
    Mat img_gray = Mat(rows,cols,CV_8UC1,Scalar(0));
    for(int i = 0; i < rows; i++)
    {
        for(int j = 0; j<cols; j++)
        {
            img_gray.at<uchar>(i,j) = (0.1*img_src.at<Vec3b>(i,j)[0])+(0.6*img_src.at<Vec3b>(i,j)[1])+(0.3*img_src.at<Vec3b>(i,j)[2]);//粗略参数
        }
    }
    gray_img = img_gray;
}


void my_Harris::IMG_GRAD(Mat img, Mat &x_grad, Mat &y_grad)
{
    Mat img_src = img.clone();
    int rows = img_src.rows;
    int cols = img_src.cols;
    Mat x_g = Mat(rows,cols,CV_8UC1,Scalar(0));
    Mat y_g = Mat(rows,cols,CV_8UC1,Scalar(0));
    for(int i = 1; i <  rows-1; i++)
    {
        for(int j = 1; j < cols -1; j++)
        {
　　　　　　　//使用sobel算子求梯度
            x_g.at<uchar>(i,j) = abs((img_src.at<uchar>(i+1,j-1) - img_src.at<uchar>(i-1,j-1))+ 2*(img_src.at<uchar>(i+1,j) - img_src.at<uchar>(i-1,j)) + (img_src.at<uchar>(i+1,j+1) - img_src.at<uchar>(i-1,j+1)))/3;
            y_g.at<uchar>(i,j) = abs((img_src.at<uchar>(i-1,j+1) - img_src.at<uchar>(i-1,j-1))+2*(img_src.at<uchar>(i,j+1) - img_src.at<uchar>(i,j-1)) + (img_src.at<uchar>(i+1,j+1) - img_src.at<uchar>(i+1,j-1)))/3;
        }
    }
    x_grad = x_g;
    y_grad = y_g;
}

void my_Harris::GRAD_MULTI(Mat x_grad, Mat y_grad, Mat &p_x, Mat &p_y, Mat &p_xy)
{
    Mat x_g = x_grad.clone();
    Mat y_g = y_grad.clone();
    int rows = x_g.rows;
    int cols = x_g.cols;
    Mat px = Mat(rows,cols,CV_8UC1,Scalar(0));
    Mat py = Mat(rows,cols,CV_8UC1,Scalar(0));
    Mat pxy = Mat(rows,cols,CV_8UC1,Scalar(0));
    for(int i = 0; i < rows; i++)
    {
        for(int j = 0; j < cols; j++)
        {
            px.at<uchar>(i,j) = x_g.at<uchar>(i,j)*x_g.at<uchar>(i,j);
            py.at<uchar>(i,j) = y_g.at<uchar>(i,j)*y_g.at<uchar>(i,j);
            pxy.at<uchar>(i,j) = x_g.at<uchar>(i,j)*y_g.at<uchar>(i,j);
        }
    }
    p_x = px;
    p_y = py;
    p_xy = pxy;
}

void my_Harris::GAUSS_WEI(Mat p_x, Mat p_y, Mat p_xy, Mat &A, Mat &B, Mat &C)
{
    Mat Ix = p_x.clone();
    Mat Iy = p_y.clone();
    Mat Ixy = p_xy.clone();
    int rows = Ix.rows;
    int cols = Ix.cols;
    Mat g_Ix = Mat(rows,cols,CV_8UC1,Scalar(0));
    Mat g_Iy = Mat(rows,cols,CV_8UC1,Scalar(0));
    Mat g_Ixy = Mat(rows,cols,CV_8UC1,Scalar(0));
    int k = 3;
    Mat gauss_plate = Mat(k,k,CV_8UC1,Scalar(0));
    //GET_GAUSS(3,0.5, gauss_plate);
    g_Ix = GAUSS_PRO(Ix);
    g_Iy = GAUSS_PRO(Iy);
    g_Ixy = GAUSS_PRO(Ixy);
    A = g_Ix;
    B = g_Iy;
    C = g_Ixy;
}


Mat my_Harris::GAUSS_PRO(Mat g)
{
    Mat Ix = g.clone();
    float gauss_plate[3][3]={0.0113437, 0.0838195, 0.0113437, 0.0838195, 0.619347, 0.0838195,0.0113437,0.0838195,0.0113437};
    int rows = Ix.rows;
    int cols = Ix.cols;
    Mat g_f = Mat(rows,cols,CV_8UC1,Scalar(0));
    for(int i = 1; i < rows-1; i++)
    {
        for(int j = 1; j < cols-1; j++)
        {
            g_f.at<uchar>(i,j) = Ix.at<uchar>(i-1,j-1)*gauss_plate[0][0] + Ix.at<uchar>(i,j-1)*gauss_plate[1][0] + Ix.at<uchar>(i+1,j-1)*gauss_plate[2][0]
                    + Ix.at<uchar>(i-1,j)*gauss_plate[0][1] + Ix.at<uchar>(i,j-1)*gauss_plate[1][1] + Ix.at<uchar>(i+1,j-1)*gauss_plate[2][1]
                    + Ix.at<uchar>(i-1,j+1)*gauss_plate[0][2] + Ix.at<uchar>(i,j+1)*gauss_plate[1][2] + Ix.at<uchar>(i+1,j+1)*gauss_plate[2][2];

         }
    }
    return g_f;
}
void my_Harris::GET_RESPONSE(Mat A_, Mat B_, Mat C_, Mat &R)
{
    Mat A = A_.clone();
    Mat B = B_.clone();
    Mat C = C_.clone();
    int rows = A.rows;
    int cols = B.cols;
    Mat M_R = Mat(rows, cols, CV_8UC1, Scalar(0));
    for(int i = 0; i < rows; i++)
    {
        for(int j = 0; j < cols; j++)
        {
            float k = 0.15;
            float a = A.at<uchar>(i,j);
            float b = B.at<uchar>(i,j);
            float c = C.at<uchar>(i,j);
            float r = (a*b-c*c-k*((a+b)*(a+b)));
            M_R.at<uchar>(i,j) = int(r);
            if(abs(int(r)) > 48000)
            {
                M_R.at<uchar>(i,j) = 255;
            }
            else
            {
               M_R.at<uchar>(i,j) = 0;
            }
        }
    }
    R = M_R;
}

void my_Harris::FILTER_THRESH(Mat R, Mat img, Mat &F_R)
{
    Mat img_src = R.clone();
    int rows = img_src.rows;
    int cols = img_src.cols;
    for(int i = 0; i < rows; i++)
    {
        for(int j = 0; j < cols; j++)
        {
            if(abs(img_src.at<uchar>(i,j))>254)
            {
                img.at<Vec3b>(i,j)[0]=0;
                img.at<Vec3b>(i,j)[1]=0;
                img.at<Vec3b>(i,j)[2]=255;
            }


        }
    }
    imshow("img_src",img);
    waitKey(0);
}

my_harris.h

#ifndef MY_HARRIS_H
#define MY_HARRIS_H
#include "iostream"
#include "opencv2/core.hpp"
#include "opencv2/highgui.hpp"

using namespace cv;
using namespace std;
class my_Harris
{
public:
    my_Harris();
    void RGB2GRAY(Mat rgb_img, Mat &gray_img);
    void IMG_GRAD(Mat img, Mat &x_grad, Mat &y_grad);
    void GRAD_MULTI(Mat x_grad, Mat y_grad, Mat &p_x, Mat &p_y, Mat &p_xy);
    void GAUSS_WEI(Mat p_x, Mat p_y, Mat p_xy, Mat &A, Mat &B, Mat &C);
    Mat GAUSS_PRO( Mat g);
    void GET_RESPONSE(Mat A_, Mat B_, Mat C_, Mat &R);
    void FILTER_THRESH(Mat R, Mat img, Mat &F_R);
};

#endif // MY_HARRIS_H

效果图：

      

实现方法简单，可继续改进！

参考：https://www.cnblogs.com/zyly/p/9508131.html#_label5_0；