特征點檢測廣泛應用到目標匹配、目標跟蹤、三維重建等應用中，在進行目標建模時會對圖像進行目標特征的提取，常用的有顏色、角點、特征點、輪廓、紋理等特征。現在開始講解常用的特征點檢測，其中Harris角點檢測是特征點檢測的基礎，提出了應用鄰近像素點灰度差值概念，從而進行判斷是否為角點、邊緣、平滑區域。Harris角點檢測原理是利用移動的窗口在圖像中計算灰度變化值，其中關鍵流程包括轉化為灰度圖像、計算差分圖像、高斯平滑、計算局部極值、確認角點。
一、基礎知識
圖像的變化類型：
在特征點檢測中經常提出尺度不變、旋轉不變、抗噪聲影響等，這些是判斷特征點是否穩定的指標。
性能較好的角點：
- 檢測出圖像中“真實”的角點
- 准確的定位性能
- 很高的重復檢測率
- 噪聲的魯棒性
- 較高的計算效率
角點的類型：
基於圖像灰度的方法通過計算點的曲率及梯度來檢測角點，避免了第一類方法存在的缺陷，此類方法主要有Moravec算子、Forstner算子、Harris算子、SUSAN算子等。
二、Harris角點原理
1、算法思想
角點原理來源於人對角點的感性判斷，即圖像在各個方向灰度有明顯變化。算法的核心是利用局部窗口在圖像上進行移動判斷灰度發生較大的變化，所以此窗口用於計算圖像的灰度變化為：[-1,0,1;-1,0,1;-1,0,1][-1,-1,-1;0,0,0;1,1,1]。人各個方向上移動這個特征的小窗口，如圖3中窗口內區域的灰度發生了較大的變化，那么就認為在窗口內遇到了角點。如圖1中，窗口內圖像的灰度沒有發生變化，那么窗口內就不存在角點；如果窗口在某一個方向移動時，窗口內圖像的灰度發生了較大的變化，而在另一些方向上沒有發生變化，那么，窗口內的圖像可能就是一條直線的線段。
2、數學模型
根據算法思想，構建數學模型，計算移動窗口的的灰度差值。
為了減小計算量，利用泰勒級數進行簡化公式：
上圖中W函數表示窗口函數，M矩陣為偏導數矩陣。對於矩陣可以進行對稱矩陣的變化，假設利用兩個特征值進行替代，其幾何含義類似下圖中的表達。在幾何
模型中通過判斷兩個特征值的大小，來判定像素的屬性。
M為梯度的協方差矩陣 ，在實際應用中為了能夠應用更好的編程，定義了角點響應函數R，通過判定R大小來判斷像素是否為角點。
R取決於M的特征值，對於角點|R|很大，平坦的區域|R|很小，邊緣的R為負值。
3、算法流程
1.利用水平，豎直差分算子對圖像的每個像素進行濾波以求得Ix,Iy,進而求得M中的四個元素的值。
算法源碼：
代碼中如果array為-1,0,1,-1,0,1,-1,0,1}則是求解X方向的，如果為{-1,-1,-1,0,0,0,1,1,1}為Y方向的，則Ix和Iy求解結束
//C代碼
//這里的size/2是為了不把圖像邊界算進去。
//array也就是3*3的窗口函數為：double dx[9]={-1,0,1,-1,0,1,-1,0,1};或者 //定義垂直方向差分算子並求Iy
double dy[9]={-1,-1,-1,0,0,0,1,1,1};
CvMat *mbys(CvMat *mat,int xwidth,int ywidth,double *array,int size1,int size2)//size
{
int i,j;
int i1,j1;
int px,py;
int m;
CvMat *mat1;
mat1=cvCloneMat(mat);
//為了去除邊界，從框體一半開始遍歷
for(i=size1/2;i<ywidth-size1/2;i++)
for(j=size2/2;j<xwidth-size2/2;j++)
{
m=0;
for(i1=0;i1<size1;i1++)
for(j1=0;j1<size2;j1++)
{
px=i-size1/2+i1;
py=j-size2/2+j1;
//CV_MAT_ELEM訪問矩陣元素
//每個元素和窗體函數遍歷相加
m+=CV_MAT_ELEM(*mat,double,px,py)*array[i1*size1+j1];
}
//賦值
CV_MAT_ELEM(*mat1,double,i,j)=m;
}
return mat1;
}
求解IX2相對比較簡單，像素相乘即可。下面為基於opencv的C++版本，很是簡單
//構建模板 Mat xKernel = (Mat_<double>(1,3) << -1, 0, 1); Mat yKernel = xKernel.t(); //計算IX和IY Mat Ix,Iy; //可參考filter2d的定義 filter2D(gray, Ix, CV_64F, xKernel); filter2D(gray, Iy, CV_64F, yKernel); //計算其平方 Mat Ix2,Iy2,Ixy; Ix2 = Ix.mul(Ix); Iy2 = Iy.mul(Iy); Ixy = Ix.mul(Iy);
2.對M的四個元素進行高斯平滑濾波，為的是消除一些不必要的孤立點和凸起，得到新的矩陣M。
//本例中使用5×5的高斯模板
//計算模板參數
//int gausswidth=5;
//double sigma=0.8;
double *g=new double[gausswidth*gausswidth];
for(i=0;i<gausswidth;i++)//定義模板
for(j=0;j<gausswidth;j++)
g[i*gausswidth+j]=exp(-((i-int(gausswidth/2))*(i-int(gausswidth/2))+(j-int(gausswidth/2))*(j-int(gausswidth/2)))/(2*sigma));
//歸一化：使模板參數之和為1（其實此步可以省略）
double total=0;
for(i=0;i<gausswidth*gausswidth;i++)
total+=g[i];
for(i=0;i<gausswidth;i++)
for(j=0;j<gausswidth;j++)
g[i*gausswidth+j]/=total;
//進行高斯平滑
mat_Ix2=mbys(mat_Ix2,cxDIB,cyDIB,g,gausswidth,gausswidth);
mat_Iy2=mbys(mat_Iy2,cxDIB,cyDIB,g,gausswidth,gausswidth);
mat_Ixy=mbys(mat_Ixy,cxDIB,cyDIB,g,gausswidth,gausswidth);
利用opencv接口則是很簡單
//構建高斯函數 Mat gaussKernel = getGaussianKernel(7, 1); //卷積計算 filter2D(Ix2, Ix2, CV_64F, gaussKernel); filter2D(Iy2, Iy2, CV_64F, gaussKernel); filter2D(Ixy, Ixy, CV_64F, gaussKernel);
3.接下來利用M計算對應每個像素的角點響應函數R，即：
也可以使用改進的R：
R=[Ix^2*Iy^2-(Ix*Iy)^2]/(Ix^2+Iy^2);里面沒有隨意給定的參數k，取值應當比第一個令人滿意。
//計算cim：即cornerness of image，我們把它稱做‘角點量’
CvMat *mat_cim;
mat_cim=mbcim(mat_Ix2,mat_Iy2,mat_Ixy,cxDIB,cyDIB);
//用來求得響應度
CvMat *mbcim(CvMat *mat1,CvMat *mat2,CvMat *mat3,int xwidth,int ywidth)
{
int i,j;
CvMat *mat;
mat=cvCloneMat(mat1);
for(i = 0; i <ywidth; i++)
{
for(j = 0; j < xwidth; j++)
{
//注意：要在分母中加入一個極小量以防止除數為零溢出
CV_MAT_ELEM(*mat,double,i,j)=(CV_MAT_ELEM(*mat1,double,i,j)*CV_MAT_ELEM(*mat2,double,i,j)-
CV_MAT_ELEM(*mat3,double,i,j)*CV_MAT_ELEM(*mat3,double,i,j))/
(CV_MAT_ELEM(*mat1,double,i,j)+CV_MAT_ELEM(*mat2,double,i,j)+0.000001);
}
}
return mat;
}
//opencv代碼
Mat cornerStrength(gray.size(), gray.type());
for (int i = 0; i < gray.rows; i++)
{
for (int j = 0; j < gray.cols; j++)
{
double det_m = Ix2.at<double>(i,j) * Iy2.at<double>(i,j) - Ixy.at<double>(i,j) * Ixy.at<double>(i,j);
double trace_m = Ix2.at<double>(i,j) + Iy2.at<double>(i,j);
cornerStrength.at<double>(i,j) = det_m - alpha * trace_m *trace_m;
}
}
4、局部極大值抑制，同時選取其極大值
//--------------------------------------------------------------------------
// 第四步：進行局部非極大值抑制
//--------------------------------------------------------------------------
CvMat *mat_locmax;
//const int size=7;
mat_locmax=mblocmax(mat_cim,cxDIB,cyDIB,size);
// cout<<CV_MAT_ELEM(*mat_locmax,double,cyDIB-1,cxDIB-1)<<endl;
//用來求得局部極大值
CvMat *mblocmax(CvMat *mat1,int xwidth,int ywidth,int size)
{
int i,j;
double max=-1000;
int i1,j1;
int px,py;
CvMat *mat;
mat=cvCloneMat(mat1);
for(i=size/2;i<ywidth-size/2;i++)
for(j=size/2;j<xwidth-size/2;j++)
{
max=-10000;
for(i1=0;i1<size;i1++)
for(j1=0;j1<size;j1++)
{
px=i-size/2+i1;
py=j-size/2+j1;
if(CV_MAT_ELEM(*mat1,double,px,py)>max)
max=CV_MAT_ELEM(*mat1,double,px,py);
}
if(max>0)
CV_MAT_ELEM(*mat,double,i,j)=max;
else
CV_MAT_ELEM(*mat,double,i,j)=0;
}
return mat;
}
在opencv中應用maxminloc函數
// threshold double maxStrength; minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL); Mat dilated; Mat localMax; //默認3*3核膨脹，膨脹之后，除了局部最大值點和原來相同，其它非局部最大值點被 //3*3鄰域內的最大值點取代 dilate(cornerStrength, dilated, Mat()); //與原圖相比，只剩下和原圖值相同的點，這些點都是局部最大值點，保存到localMax compare(cornerStrength, dilated, localMax, CMP_EQ);
5.在矩陣R中，同時滿足R(i,j)大於一定閾值threshold和R(i,j)是某領域內的局部極大值，則被認為是角點。
cv::Mat cornerMap; // 根據角點響應最大值計算閾值 threshold= qualityLevel*maxStrength; cv::threshold(cornerStrength,cornerTh, threshold,255,cv::THRESH_BINARY); // 轉為8-bit圖 cornerTh.convertTo(cornerMap,CV_8U);// 和局部最大值圖與，剩下角點局部最大值圖，即：完成非最大值抑制 cv::bitwise_and(cornerMap,localMax,cornerMap);
imgDst = cornerMap.clone();
for( int y = 0; y < cornerMap.rows; y++ )
{
const uchar* rowPtr = cornerMap.ptr<uchar>(y);
for( int x = 0; x < cornerMap.cols; x++ )
{
// 非零點就是角點
if (rowPtr[x])
{
points.push_back(cv::Point(x,y));
}
}
}
三、算法源碼
1、C版本
里面函數最基礎的代碼解釋和注釋
////////////////////////////////////////////////////////////////////// // Construction/Destruction ////////////////////////////////////////////////////////////////////// #define B(image,x,y) ((uchar *)(image->imageData+image->widthStep*(y)))[(x)*3] #define G(image,x,y) ((uchar *)(image->imageData+image->widthStep*(y)))[(x)*3+1] #define R(image,x,y) ((uchar *)(image->imageData+image->widthStep*(y)))[(x)*3+2] #define S(image,x,y) ((uchar *)(image->imageData+image->widthStep*(y)))[(x)] //卷積計算求Ix,Iy,以及濾波 //a指向的數組是size1*size2大小的...求導 CvMat *mbys(CvMat *mat,int xwidth,int ywidth,double *a,int size1,int size2) { int i,j; int i1,j1; int px,py; int m; CvMat *mat1; mat1=cvCloneMat(mat); for(i=size1/2;i<ywidth-size1/2;i++) for(j=size2/2;j<xwidth-size2/2;j++) { m=0; for(i1=0;i1<size1;i1++) for(j1=0;j1<size2;j1++) { px=i-size1/2+i1; py=j-size2/2+j1; //CV_MAT_ELEM訪問矩陣元素 m+=CV_MAT_ELEM(*mat,double,px,py)*a[i1*size1+j1]; } CV_MAT_ELEM(*mat1,double,i,j)=m; } return mat1; } //計算Ix2,Iy2,Ixy CvMat *mbxy(CvMat *mat1,CvMat *mat2,int xwidth,int ywidth) { int i,j; CvMat *mat3; mat3=cvCloneMat(mat1); for(i=0;i<ywidth;i++) for (j=0;j<xwidth;j++) { CV_MAT_ELEM(*mat3,double,i,j)=CV_MAT_ELEM(*mat1,double,i,j)*CV_MAT_ELEM(*mat2,double,i,j); } return mat3; } //用來求得響應度 CvMat *mbcim(CvMat *mat1,CvMat *mat2,CvMat *mat3,int xwidth,int ywidth) { int i,j; CvMat *mat; mat=cvCloneMat(mat1); for(i = 0; i <ywidth; i++) { for(j = 0; j < xwidth; j++) { //注意：要在分母中加入一個極小量以防止除數為零溢出 CV_MAT_ELEM(*mat,double,i,j)=(CV_MAT_ELEM(*mat1,double,i,j)*CV_MAT_ELEM(*mat2,double,i,j)- CV_MAT_ELEM(*mat3,double,i,j)*CV_MAT_ELEM(*mat3,double,i,j))/ (CV_MAT_ELEM(*mat1,double,i,j)+CV_MAT_ELEM(*mat2,double,i,j)+0.000001); } } return mat; } //用來求得局部極大值 CvMat *mblocmax(CvMat *mat1,int xwidth,int ywidth,int size) { int i,j; double max=-1000; int i1,j1; int px,py; CvMat *mat; mat=cvCloneMat(mat1); for(i=size/2;i<ywidth-size/2;i++) for(j=size/2;j<xwidth-size/2;j++) { max=-10000; for(i1=0;i1<size;i1++) for(j1=0;j1<size;j1++) { px=i-size/2+i1; py=j-size/2+j1; if(CV_MAT_ELEM(*mat1,double,px,py)>max) max=CV_MAT_ELEM(*mat1,double,px,py); } if(max>0) CV_MAT_ELEM(*mat,double,i,j)=max; else CV_MAT_ELEM(*mat,double,i,j)=0; } return mat; } //用來確認角點 CvMat *mbcorner(CvMat *mat1,CvMat *mat2,int xwidth,int ywidth,int size,double thresh) { CvMat *mat; int i,j; mat=cvCreateMat(ywidth,xwidth,CV_32FC1); for(i=size/2;i<ywidth-size/2;i++) for(j=size/2;j<xwidth-size/2;j++) { if(CV_MAT_ELEM(*mat1,double,i,j)==CV_MAT_ELEM(*mat2,double,i,j))//首先取得局部極大值 if(CV_MAT_ELEM(*mat1,double,i,j)>thresh)//然后大於這個閾值 CV_MAT_ELEM(*mat,int,i,j)=255;//滿足上兩個條件，才是角點！ else CV_MAT_ELEM(*mat,int,i,j)=0; } return mat; } CvPoint* CHarris::harris_features(IplImage *src,int gausswidth,double sigma,int size,int threshold) { CvMat *mat_I,*mat_Ix,*mat_Iy,*mat_Ixy,*mat_Ix2,*mat_Iy2;//相應的矩陣 IplImage *pImgGray=NULL; //灰度圖像 IplImage *dst=NULL; //目標圖像 IplImage *pImgDx=NULL; //水平梯度卷積后的圖像 IplImage *pImgDy=NULL; //豎起梯度卷積后的圖像 IplImage *pImgDx2=NULL;//Ix2圖像 IplImage *pImgDy2=NULL;//Iy2圖像 IplImage *pImgDxy=NULL;//Ixy圖像 pImgGray=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1); dst=cvCreateImage(cvGetSize(src),src->depth,3); pImgDx=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);//創建圖像 pImgDy=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1); pImgDx2=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1); pImgDy2=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1); pImgDxy=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1); const int cxDIB=src->width ; // 圖像寬度 const int cyDIB=src->height; // 圖像高度 double *I=new double[cxDIB*cyDIB]; cvCvtColor(src,pImgGray,CV_RGB2GRAY);//灰度化 dst=cvCloneImage(src); int i,j; for(j=0;j<cyDIB;j++) for(i=0;i<cxDIB;i++) { I[j*cxDIB+i]=S(pImgGray,i,j);//將灰度圖像數值存入I中 } mat_I=cvCreateMat(cyDIB,cxDIB,CV_64FC1); cvInitMatHeader(mat_I,cyDIB,cxDIB,CV_64FC1,I);//用I來初始化相應的矩陣 // cout<<CV_MAT_ELEM(*mat_I,double,200,200)<<endl; //-------------------------------------------------------------------------- // 第一步：利用差分算子對圖像進行濾波 //-------------------------------------------------------------------------- //定義水平方向差分算子並求Ix double dx[9]={-1,0,1,-1,0,1,-1,0,1}; mat_Ix=mbys(mat_I,cxDIB,cyDIB,dx,3,3); //求Ix矩陣 // cout<<CV_MAT_ELEM(*mat_Ix,double,200,200)<<endl; //定義垂直方向差分算子並求Iy double dy[9]={-1,-1,-1,0,0,0,1,1,1}; mat_Iy=mbys(mat_I,cxDIB,cyDIB,dy,3,3);//求Iy矩陣 // cout<<CV_MAT_ELEM(*mat_Iy,double,200,200)<<endl; for(j=0;j<cyDIB;j++) for(i=0;i<cxDIB;i++) { S(pImgDx,i,j)=CV_MAT_ELEM(*mat_Ix,double,j,i);//為相應圖像賦值 S(pImgDy,i,j)=CV_MAT_ELEM(*mat_Iy,double,j,i); } mat_Ix2=mbxy(mat_Ix,mat_Ix,cxDIB,cyDIB);//計算Ix2,Iy2,Ixy矩陣 mat_Iy2=mbxy(mat_Iy,mat_Iy,cxDIB,cyDIB); mat_Ixy=mbxy(mat_Ix,mat_Iy,cxDIB,cyDIB); for(j=0;j<cyDIB;j++) for(i=0;i<cxDIB;i++) { S(pImgDxy,i,j)=CV_MAT_ELEM(*mat_Ixy,double,j,i);//為相應圖像賦值 S(pImgDx2,i,j)=CV_MAT_ELEM(*mat_Ix2,double,j,i); S(pImgDy2,i,j)=CV_MAT_ELEM(*mat_Iy2,double,j,i); } //-------------------------------------------------------------------------- // 第二步：對Ix2/Iy2/Ixy進行高斯平滑，以去除噪聲 //-------------------------------------------------------------------------- //本例中使用5×5的高斯模板 //計算模板參數 //int gausswidth=5; //double sigma=0.8; double *g=new double[gausswidth*gausswidth]; for(i=0;i<gausswidth;i++)//定義模板 for(j=0;j<gausswidth;j++) g[i*gausswidth+j]=exp(-((i-int(gausswidth/2))*(i-int(gausswidth/2))+(j-int(gausswidth/2))*(j-int(gausswidth/2)))/(2*sigma)); //歸一化：使模板參數之和為1（其實此步可以省略） double total=0; for(i=0;i<gausswidth*gausswidth;i++) total+=g[i]; for(i=0;i<gausswidth;i++) for(j=0;j<gausswidth;j++) g[i*gausswidth+j]/=total; //進行高斯平滑 mat_Ix2=mbys(mat_Ix2,cxDIB,cyDIB,g,gausswidth,gausswidth); mat_Iy2=mbys(mat_Iy2,cxDIB,cyDIB,g,gausswidth,gausswidth); mat_Ixy=mbys(mat_Ixy,cxDIB,cyDIB,g,gausswidth,gausswidth); //-------------------------------------------------------------------------- // 第三步：計算角點量 //-------------------------------------------------------------------------- //計算cim：即cornerness of image，我們把它稱做‘角點量’ CvMat *mat_cim; mat_cim=mbcim(mat_Ix2,mat_Iy2,mat_Ixy,cxDIB,cyDIB); // cout<<CV_MAT_ELEM(*mat_cim,double,cyDIB-1,cxDIB-1)<<endl; //-------------------------------------------------------------------------- // 第四步：進行局部非極大值抑制 //-------------------------------------------------------------------------- CvMat *mat_locmax; //const int size=7; mat_locmax=mblocmax(mat_cim,cxDIB,cyDIB,size); // cout<<CV_MAT_ELEM(*mat_locmax,double,cyDIB-1,cxDIB-1)<<endl; //-------------------------------------------------------------------------- // 第五步：獲得最終角點 //-------------------------------------------------------------------------- CvMat *mat_corner; //const double threshold=4500; //int cornernum=0; mat_corner=mbcorner(mat_cim,mat_locmax,cxDIB,cyDIB,gausswidth,threshold); //CCommon CommonClass; CvPoint pt[5000]; for(j=size/2;j<cyDIB-size/2;j++) for(i=size/2;i<cxDIB-size/2;i++) { if(CV_MAT_ELEM(*mat_corner,int,j,i)==255) { pt[cornerno].x=i; pt[cornerno].y=j; cornerno++; // CommonClass.DrawCross(showImg2,pt,CV_RGB(0,0,255),1,4); // cvCircle(dst,pt,2,CV_RGB(255,0,0),1,8,0); // cout<<i<<" "<<j<<endl; } } return pt; }
2、基於opencv源碼
1 #include "imgFeat.h" 2 3 void feat::detectHarrisCorners(const Mat& imgSrc, Mat& imgDst, double alpha) 4 { 5 Mat gray; 6 if (imgSrc.channels() == 3) 7 { 8 cvtColor(imgSrc, gray, CV_BGR2GRAY); 9 } 10 else 11 { 12 gray = imgSrc.clone(); 13 } 14 gray.convertTo(gray, CV_64F); 15 16 Mat xKernel = (Mat_<double>(1,3) << -1, 0, 1); 17 Mat yKernel = xKernel.t(); 18 19 Mat Ix,Iy; 20 filter2D(gray, Ix, CV_64F, xKernel); 21 filter2D(gray, Iy, CV_64F, yKernel); 22 23 Mat Ix2,Iy2,Ixy; 24 Ix2 = Ix.mul(Ix); 25 Iy2 = Iy.mul(Iy); 26 Ixy = Ix.mul(Iy); 27 28 Mat gaussKernel = getGaussianKernel(7, 1); 29 filter2D(Ix2, Ix2, CV_64F, gaussKernel); 30 filter2D(Iy2, Iy2, CV_64F, gaussKernel); 31 filter2D(Ixy, Ixy, CV_64F, gaussKernel); 32 33 34 Mat cornerStrength(gray.size(), gray.type()); 35 for (int i = 0; i < gray.rows; i++) 36 { 37 for (int j = 0; j < gray.cols; j++) 38 { 39 double det_m = Ix2.at<double>(i,j) * Iy2.at<double>(i,j) - Ixy.at<double>(i,j) * Ixy.at<double>(i,j); 40 double trace_m = Ix2.at<double>(i,j) + Iy2.at<double>(i,j); 41 cornerStrength.at<double>(i,j) = det_m - alpha * trace_m *trace_m; 42 } 43 } 44 // threshold 45 double maxStrength; 46 minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL); 47 Mat dilated; 48 Mat localMax; 49 dilate(cornerStrength, dilated, Mat()); 50 compare(cornerStrength, dilated, localMax, CMP_EQ); 51 52 53 Mat cornerMap; 54 double qualityLevel = 0.01; 55 double thresh = qualityLevel * maxStrength; 56 cornerMap = cornerStrength > thresh; 57 bitwise_and(cornerMap, localMax, cornerMap); 58 59 imgDst = cornerMap.clone(); 60 61 } 62 63 void feat::drawCornerOnImage(Mat& image, const Mat&binary) 64 { 65 Mat_<uchar>::const_iterator it = binary.begin<uchar>(); 66 Mat_<uchar>::const_iterator itd = binary.end<uchar>(); 67 for (int i = 0; it != itd; it++, i++) 68 { 69 if (*it) 70 circle(image, Point(i%image.cols, i / image.cols), 3, Scalar(0, 255, 0), 1); 71 } 72 }
3、opencv中Harris接口
OpenCV的Hairrs角點檢測的函數為cornerHairrs()，但是它的輸出是一幅浮點值圖像，浮點值越高，表明越可能是特征角點，我們需要對圖像進行閾值化。
C++: void cornerHarris(InputArray src, OutputArray dst, int blockSize, int apertureSize, double k, int borderType = BORDER_DEFAULT);
- src – 輸入的單通道8-bit或浮點圖像。
- dst – 存儲着Harris角點響應的圖像矩陣，大小與輸入圖像大小相同，是一個浮點型矩陣。
- blockSize – 鄰域大小。
- apertureSize – 擴展的微分算子大。
- k – 響應公式中的，參數$\alpha$。
- boderType – 邊界處理的類型。
int main()
{
Mat image = imread("../buliding.png");
Mat gray;
cvtColor(image, gray, CV_BGR2GRAY);
Mat cornerStrength;
cornerHarris(gray, cornerStrength, 3, 3, 0.01);
threshold(cornerStrength, cornerStrength, 0.001, 255, THRESH_BINARY);
return 0;
}
從上面上間一幅圖像我們可以看到，有很多角點都是粘連在一起的，我們下面通過加入非極大值抑制來進一步去除一些粘在一起的角點。
非極大值抑制原理是，在一個窗口內，如果有多個角點則用值最大的那個角點，其他的角點都刪除，窗口大小這里我們用3*3，程序中通過圖像的膨脹運算來達到檢測極大值的目的，因為默認參數的膨脹運算就是用窗口內的最大值替代當前的灰度值。
int main() { Mat image = imread("buliding.png"); Mat gray; cvtColor(image, gray, CV_BGR2GRAY); Mat cornerStrength; cornerHarris(gray, cornerStrength, 3, 3, 0.01); double maxStrength; double minStrength; // 找到圖像中的最大、最小值 minMaxLoc(cornerStrength, &minStrength, &maxStrength); Mat dilated; Mat locaMax; // 膨脹圖像，最找出圖像中全部的局部最大值點 dilate(cornerStrength, dilated, Mat()); // compare是一個邏輯比較函數，返回兩幅圖像中對應點相同的二值圖像 compare(cornerStrength, dilated, locaMax, CMP_EQ); Mat cornerMap; double qualityLevel = 0.01; double th = qualityLevel*maxStrength; // 閾值計算 threshold(cornerStrength, cornerMap, th, 255, THRESH_BINARY); cornerMap.convertTo(cornerMap, CV_8U); // 逐點的位運算黑色圖片顯示的結果 bitwise_and(cornerMap, locaMax, cornerMap); drawCornerOnImage(image, cornerMap); namedWindow("result"); imshow("result", image); waitKey(); return 0; } void drawCornerOnImage(Mat& image, const Mat&binary) { Mat_<uchar>::const_iterator it = binary.begin<uchar>(); Mat_<uchar>::const_iterator itd = binary.end<uchar>(); for (int i = 0; it != itd; it++, i++) { if (*it) circle(image, Point(i%image.cols, i / image.cols), 3, Scalar(0, 255, 0), 1); } }
四、Harris角點性質
1、閾值決定檢測點數量
增大$\alpha$的值，將減小角點響應值$R$，降低角點檢測的靈性，減少被檢測角點的數量；減小$\alpha$值，將增大角點響應值$R$，增加角點檢測的靈敏性，增加被檢測角點的數量
2、Harris角點檢測算子對亮度和對比度的變化不敏感
這是因為在進行Harris角點檢測時，使用了微分算子對圖像進行微分運算，而微分運算對圖像密度的拉升或收縮和對亮度的抬高或下降不敏感。換言之，對亮度和對比度的仿射變換並不改變Harris響應的極值點出現的位置，但是，由於閾值的選擇，可能會影響角點檢測的數量。
3. Harris角點檢測算子具有旋轉不變性
Harris角點檢測算子使用的是角點附近的區域灰度二階矩矩陣。而二階矩矩陣可以表示成一個橢圓，橢圓的長短軸正是二階矩矩陣特征值平方根的倒數。當特征橢圓轉動時，特征值並不發生變化，所以判斷角點響應值也不發生變化，由此說明Harris角點檢測算子具有旋轉不變性。
4. Harris角點檢測算子不具有尺度不變性
如下圖所示，當右圖被縮小時，在檢測窗口尺寸不變的前提下，在窗口內所包含圖像的內容是完全不同的。左側的圖像可能被檢測為邊緣或曲線，而右側的圖像則可能被檢測為一個角點。
五、函數備注
1、compare
功能：兩個數組之間或者一個數組和一個常數之間的比較
結構：
void compare(InputArray src1, InputArray src2, OutputArray dst, int cmpop)
src1 ：第一個數組或者標量，如果是數組，必須是單通道數組。
src2 ：第二個數組或者標量，如果是數組，必須是單通道數組。
dst ：輸出數組，和輸入數組有同樣的size和type=CV_8UC1
cmpop ：
標志指明了元素之間的對比關系
CMP_EQ src1 相等 src2.
CMP_GT src1 大於 src2.
CMP_GE src1 大於或等於 src2.
CMP_LT src1 小於 src2.
CMP_LE src1 小於或等於 src2.
CMP_NE src1 不等於 src2.
如果對比結果為true，那么輸出數組對應元素的值為255，否則為0
//獲取角點圖
Mat getCornerMap(double qualityLevel) {
Mat cornerMap;
// 根據角點響應最大值計算閾值
thresholdvalue= qualityLevel*maxStrength;
threshold(cornerStrength,cornerTh,
thresholdvalue,255,cv::THRESH_BINARY);
// 轉為8-bit圖
cornerTh.convertTo(cornerMap,CV_8U);
// 和局部最大值圖與，剩下角點局部最大值圖，即：完成非最大值抑制
bitwise_and(cornerMap,localMax,cornerMap);
return cornerMap;
}
2、bitwise_and（位運算函數）
功能：計算兩個數組或數組和常量之間與的關系
結構：
void bitwise_and(InputArray src1, InputArray src2, OutputArray dst, InputArray mask=noArray())
src1 ：第一個輸入的數組或常量
src2 ：第二個輸入的數組或常量
dst ：輸出數組，和輸入數組有同樣的size和type
mask ：可選擇的操作掩碼，為8位單通道數組，指定了輸出數組哪些元素可以被改變，哪些不可以
操作過程為：
如果為多通道數組，每個通道單獨處理
3、filter2D
OpenCV中的內部函數filter2D可以直接用來做圖像卷積操作
void filter2D(InputArray src, OutputArray dst, int ddepth, InputArray kernel, Point anchor=Point(-1,-1), double delta=0, int borderType=BORDER_DEFAULT )
六、參考文章
2、尺度空間理論