為便於理解，先簡要介紹角點的概念和角點檢測背景
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); }