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；