Harris角点检测原理及实现

void goodFeaturesToTrack( InputArray image, OutputArray corners,
                          int maxCorners, double qualityLevel, 
                          double minDistance,InputArray mask=noArray(), 
                          int blockSize=3,  bool useHarrisDetector=false, 
                          double k=0.04 );　　//参数分别为：输入图像，输出角点,检测到的角点的数量的最大值，角点的质量等级（可以看做R的一个转化）,　　　　
                                　　　　　　　　　　两个角点间的最小距离，检测区域（即有无roi区域），窗口的大小，是否使用Harris角点检测(False为Shi-Tomasi算子)
　　　　　　　　　　　　　　　　　　　　　　　　　　　　　参数k

#include <opencv2/core/core.hpp>

#include <opencv2/highgui/highgui.hpp>

#include <opencv2/imgproc/imgproc.hpp>

#include <iostream>

 

using namespace cv;

using namespace std;

 

 

Mat image;

Mat imageGray;

int thresh=200;

int MaxThresh=255;

 

void Trackbar(int,void*);  //阈值控制

 

int main()  

{  

    image=imread("Test.bmp");

    cvtColor(image,imageGray,CV_RGB2GRAY);

    GaussianBlur(imageGray,imageGray,Size(5,5),1); // 滤波

    namedWindow("Corner Detected");

    createTrackbar("threshold：","Corner Detected",&thresh,MaxThresh,Trackbar);

    imshow("Corner Detected",image);

    Trackbar(0,0);

    waitKey();

    return 0;

}  

 

void Trackbar(int,void*)

{

    Mat dst,dst8u,dstshow,imageSource;

    dst=Mat::zeros(image.size(),CV_32FC1);  

    imageSource=image.clone();

    cornerHarris(imageGray,dst,3,3,0.04,BORDER_DEFAULT);

    normalize(dst,dst8u,0,255,CV_MINMAX);  //归一化

    convertScaleAbs(dst8u,dstshow);

    imshow("dst",dstshow);  //dst显示

    for(int i=0;i<image.rows;i++)

    {

        for(int j=0;j<image.cols;j++)

        {

            if(dstshow.at<uchar>(i,j)>thresh)  //阈值判断

            {

                circle(imageSource,Point(j,i),2,Scalar(0,0,255),2); //标注角点

            }

        }

    }

    imshow("Corner Detected",imageSource);

}

//--------------------【cornerHarris源码分析】------------------------------------

//Harries角点实现函数,截取cornerHarries中的关键代码做了简化

void myHarries( const Mat& src, Mat& eigenv, int block_size, int aperture_size, double k=0.)

{

    eigenv.create(src.size(), CV_32F);

    Mat Dx, Dy;

    //soble operation get Ix, Iy

    Sobel(src, Dx, CV_32F, 1, 0, aperture_size);

    Sobel(src, Dy, CV_32F, 0, 1, aperture_size);

 

    //get covariance matrix

    Size size = src.size();

    Mat cov(size, CV_32FC3);     //创建一个三通道cov矩阵分别存储[Ix*Ix, Ix*Iy, Iy*Iy];

 

    for( int i=0; i < size.height;i++)

    {

        float* cov_data = cov.ptr<float>(i);

        const float* dxdata = Dx.ptr<float>(i);

        const float* dydata = Dy.ptr<float>(i);

 

        for( int j=0; j < size.width; j++ )

        {

            float dx = dxdata[j];

            float dy = dydata[j];

 

            cov_data[j*3] = dx * dx;    //即  Ix*Ix

            cov_data[j*3 + 1] = dx*dy;  //即  Ix*Iy

            cov_data[j*3 + 2] = dy*dy;  //即  Iy*Iy

        }

    }

 

 

    //方框滤波W(x,y)卷积, 也可用高斯核加权...

    //W(Y,Y)与矩阵cov卷积运算得到 H 矩阵,后面通过H矩阵的特征值决定是否是角点

    boxFilter(cov, cov, cov.depth(), Size(block_size,block_size),Point(-1,-1),false);

 

    //cale Harris

    size = cov.size();

    if( cov.isContinuous() && eigenv.isContinuous())

    {

        size.width *= size.width;

        size.height = 1;

        //cout << "yes"<<endl;

    }

 

    //此处计算响应R= det(H) - k*trace(H)*trace(H);

    for (int i = 0; i < size.height; i++)

    {

        const float* covPtr = cov.ptr<float>(i);

        float* dstPtr = eigenv.ptr<float>(i);

 

        for( int j = 0; j < size.width; j++)

        {

            float a = covPtr[j*3];

            float b = covPtr[j*3 + 1];

            float c = covPtr[j*3 + 2];

 

            //根据公式 R = det(H) - k*trace(H)*trace(H);

            dstPtr[j] = (float)(a*c - b*b - k*(a + c)(a + c));

        }

    }

 

    double max, min;

    minMaxLoc(eigenv, &min, &max);

    //cout<< max << endl;

    double threshold = 0.1*max;

    cv::threshold(eigenv, eigenv,threshold, 1, CV_THRESH_BINARY);   //eigenv的类型是CV_32F,

    imshow("eigenv", eigenv);

}