OpenCV 角点检测(二) Harrise

阳光穿透心脏的1/2处 2022-07-11 11:07 168阅读 0赞

Harrise算子是在Moravec算子的基础上改进得到的，Moravec角点检测算子见链接：http://blog.csdn.net/chaipp0607/article/details/54649235

**Harrise算子特点**

Harrise算子将比于Moravec具有更高的时间复杂度，对噪声同样比较敏感，且存在非均匀响应。前者应用更加广泛，且具有不错的检测率。

**Harrise算子计算步骤**  
(1).利用水平与竖直差分算子对图像进行卷积操作，计算的到相应的fx和fy，根据实对称矩阵，计算对应矩阵元素的值。  
(2).利用高斯函数对矩阵M进行平滑操作，得到引得矩阵M。  
(3).对每一个像素和给定的邻域窗口，计算局部特征结果矩阵M的特征值和响应函数H。  
(4).选取响应函数H的阈值，根据非极大值抑制原理，同时满足阈值及某邻域内的局部极大值为候选点。

**Harrise算子实现**  
opencv为Harrise算子提供了cornerHarris函数。API函数接口如下：

CV_EXPORTS_W void cornerHarris( 
    InputArray src,                                                 OutputArray dst, 
    int blockSize,
    int ksize, 
    double k,
    int borderType=BORDER_DEFAULT );

它的源码路径为：…opencv\\sources\\modules\\imgproc\\src\\thresh.cpp

源码如下：

void cv::cornerHarris( InputArray _src,OutputArray _dst, int blockSize, int ksize, double k, int borderType )  
    {  
       Mat src = _src.getMat();  
       _dst.create( src.size(), CV_32F );  
       Mat dst = _dst.getMat();  
       cornerEigenValsVecs( src, dst, blockSize, ksize, HARRIS, k, borderType);  
    } 
    //cornerEigenValsVecs函数源码
    static void  
    cornerEigenValsVecs( const Mat& src,Mat& eigenv, int block_size,  
                         int aperture_size, intop_type, double k=0.,  
                         intborderType=BORDER_DEFAULT )  
    {  
    #ifdef HAVE_TEGRA_OPTIMIZATION  
       if (tegra::cornerEigenValsVecs(src, eigenv, block_size, aperture_size,op_type, k, borderType))  
           return;  
    #endif  
       
       int depth = src.depth();  
       double scale = (double)(1 << ((aperture_size > 0 ?aperture_size : 3) - 1)) * block_size;  
       if( aperture_size < 0 )  
           scale *= 2.;  
       if( depth == CV_8U )  
           scale *= 255.;  
       scale = 1./scale;  
       
       CV_Assert( src.type() == CV_8UC1 || src.type() == CV_32FC1 );  
       
       Mat Dx, Dy;  
       if( aperture_size > 0 )  
        {  
           Sobel( src, Dx, CV_32F, 1, 0, aperture_size, scale, 0, borderType );  
           Sobel( src, Dy, CV_32F, 0, 1, aperture_size, scale, 0, borderType );  
        }  
       else  
        {  
           Scharr( src, Dx, CV_32F, 1, 0, scale, 0, borderType );  
           Scharr( src, Dy, CV_32F, 0, 1, scale, 0, borderType );  
        }  
       
       Size size = src.size();  
       Mat cov( size, CV_32FC3 );  
       int i, j;  
       
       for( i = 0; i < size.height; i++ )  
        {  
           float* cov_data = (float*)(cov.data + i*cov.step);  
           const float* dxdata = (const float*)(Dx.data + i*Dx.step);  
           const float* dydata = (const float*)(Dy.data + i*Dy.step);  
       
            for( j = 0; j < size.width; j++ )  
           {  
               float dx = dxdata[j];  
               float dy = dydata[j];  
       
               cov_data[j*3] = dx*dx;  
               cov_data[j*3+1] = dx*dy;  
               cov_data[j*3+2] = dy*dy;  
           }  
        }  
       
       boxFilter(cov, cov, cov.depth(), Size(block_size, block_size),  
           Point(-1,-1), false, borderType );  
       
       if( op_type == MINEIGENVAL )  
           calcMinEigenVal( cov, eigenv );  
       else if( op_type == HARRIS )  
           calcHarris( cov, eigenv, k );  
       else if( op_type == EIGENVALSVECS )  
           calcEigenValsVecs( cov, eigenv );  
    }  
       
    }

参考opencv中的源码，自己定义一个角点检测的函数：

#include "opencv2/highgui/highgui.hpp"
    #include "opencv2/imgproc/imgproc.hpp"
    #include <iostream>
    #include <stdio.h>
    #include <stdlib.h>
    using namespace cv;
    using namespace std;
    void CornerHarris(const Mat& srcImage, Mat& result, 
    	int blockSize, int kSize, double k)
    {
    	Mat src;
    	srcImage.copyTo(src);
    	result.create(src.size(), CV_32F);
    	int depth = src.depth();
    	// 检测掩膜尺寸
    	double scale = (double)(1 << ((kSize > 0 ?
    kSize : 3) - 1)) * blockSize;
    	if (depth == CV_8U)
    		scale *= 255.;
    	scale = 1. / scale;
    	// sobel滤波
    	Mat dx, dy;
    	Sobel(src, dx, CV_32F, 1, 0, kSize, scale, 0);
    	Sobel(src, dy, CV_32F, 0, 1, kSize, scale, 0);
    	Size size = src.size();
    	cv::Mat cov(size, CV_32FC3);
    	int i, j;
    	// 求解水平与竖直梯度
    	for (i = 0; i < size.height; i++){
    		float *covData = (float*)(cov.data + i*cov.step);
    		const float *dxData = (const float*)(dx.data + i*dx.step);
    		const float *dyData = (const float*)(dy.data + i*dy.step);
    		for (j = 0; j < size.width; j++)
    		{
    			float dx_ = dxData[j];
    			float dy_ = dyData[j];
    			covData[3 * j] = dx_*dx_;
    			covData[3 * j + 1] = dx_*dy_;
    			covData[3 * j + 2] = dy_*dy_;
    		}
    	}
    	// 计算窗口内求和
    	boxFilter(cov, cov, cov.depth(), 
    		Size(blockSize, blockSize), Point(-1, -1), false);
    	// 判断图像连续性
    	if (cov.isContinuous() && result.isContinuous())
    	{
    		size.width *= size.height;
    		size.height = 1;
    	}
    	else
    		size = result.size();
    	// 计算响应函数 
    	for (i = 0; i < size.height; i++)
    	{
    		// 获取图像矩阵指针
    		float *resultData = (float*)(result.data + i*result.step);
    		const float *covData = (const float*)(cov.data + i*cov.step);
    		for (j = 0; j < size.width; j++)
    		{
    			// 焦点响应生成
    			float a = covData[3 * j];
    			float b = covData[3 * j + 1];
    			float c = covData[3 * j + 2];
    			resultData[j] = a*c - b*b - k*(a + c)*(a + c);
    		}
    	}
    }
    int main()
    {
    	cv::Mat srcImage = cv::imread("1.jpg");
    	if (!srcImage.data)
    		return -1;
    	cv::imshow("srcImage", srcImage);
    	cv::Mat srcGray, result;
    	cvtColor(srcImage, srcGray, CV_BGR2GRAY);
    	result = Mat::zeros(srcImage.size(), CV_32FC1);
    	// 角点检测参数
    	int blockSize = 2;
    	int apertureSize = 3;
    	double k = 0.04;
    	// 角点检测
    	// cornerHarris( srcGray, result, blockSize, apertureSize, k, BORDER_DEFAULT );
    	CornerHarris(srcGray, result, blockSize, apertureSize, k);
    	// 矩阵归一化
    	normalize(result, result, 0, 255, NORM_MINMAX, CV_32FC1, Mat());
    	convertScaleAbs(result, result);
    	// 绘图角点检测结果
    	for (int j = 0; j < result.rows; j++)
    	{
    		for (int i = 0; i < result.cols; i++)
    		{
    			if ((int)(result.at<uchar>(j, i)) > 150)
    			{
    				circle(srcImage, Point(i, j), 5, Scalar(0), 2, 8, 0);
    			}
    		}
    	}
    	cv::imshow("result", srcImage);
    	cv::waitKey(0);
    	return 0;
    }

原图：  
![这里写图片描述][e044121f02c0c00f8e3a2f3bd4de24dd.png]

结果图：  
![这里写图片描述][704fefaa7bddbdd442216b5228df7781.png]

[e044121f02c0c00f8e3a2f3bd4de24dd.png]: /images/20220711/fd8b3c87fa2b490eb20308ed95cea645.png
[704fefaa7bddbdd442216b5228df7781.png]: /images/20220711/639f5d6888594577a93dba65abf4e656.png