Harris角点检测原理及实现
为便于理解，先简要介绍角点的概念和角点检测背景
1 背景
角点检测大致可分为三类：基于灰度图的角点检测、基于二值化图像的角点检测和基于轮廓曲线的角点检测。Harris属于基于灰度图的角点检测。
2 Harris特征原理
2.1 概述
Harris角点检测根据窗口向多个方向，通过判断窗口内像素值有无明显变化判断有无角点。如下图：
第一幅图像中，窗口内像素值无明显变化，无角点。
第二幅图像中，窗口水平移动时有明显变化，无角点。
第三幅图中，窗口多个方向移动时有明显变化，有角点。
Harris角点检测可分为三步：梯度计算、响应值计算、角点提取。下面按步骤介绍。
2.2梯度计算：
对图像中的任意一像素点I（x,y），进行自相关平移w(x+∆x、y+∆y)得到自相关函数：
c(x,y,∆x,∆y) = ∑wh(x,y)(I(x,y)-I(x+∆x,y+∆y))2
其中 ∑w表示窗口内的点，h(x,y)表示加权函数，加权函数可根据自己需要进行修改（通过修改源代码）。
由泰勒可得：
I(x+∆x,y+∆y) = I(x,y)+∆xIx(x,y)+∆yIy(x,y)+p ≈I(x,y)+∆xIx(x,y)+∆yIy(x,y)
代入自相关函数可得（加权函数暂时忽略）：
c(x,y,∆x,∆y) = ∑w(I(x,y)-I(x+∆x,y+∆y))2 ≈ ∑w((∆xIx(x,y))2+2∆x∆yIx(x,y)Iy(x,y)+(∆yIy(x,y))2)
将上公式用图表示如下：
其中，u和v分别表示∆x和∆y，w(x,y)表示加权函数。Harris算法是通过判断像素值是否在多个方向上有明显变化可转换为为是否在x和y方向上像素值均有明显变化，再转换为Ix或Iy的变化，再转换为M矩阵的特征值λ1,λ2的变化，如下图：
2.3响应值计算 ：
上面计算不易于通过编程实现，Harris通过定义角点响应函数R的方式，用于表示一个角点的Harris响应值：
trace表示为矩阵的迹，det为矩阵的行列式（矩阵的迹：主对角线上的值相加即所有特征值的和），k为经验常数，一般取0.04~0.06。
R = λ1λ2-k(λ1+λ2)2
2.4角点提取 ：
通过2.3计算出一角点的Harris响应值大小，后通过设定的阈值t，进行判断该点是否为角点（大于t为角点，小于t则不为角点）。
2.5Harris角点检测优点：
1.亮度和对比度的变化对角点无影响。
2.Harris角点检测算子具有旋转不变性。
3.Harris角点检测具有尺度不变性。
3 代码实现
3.1 API和demo：
API为网上一个版本，现可能发生改变。
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
四、代码应用：
OpenCV内的API：
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); }
