OpenCV 角点检测(二) Harrise

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;
}

原图：



结果图：