Harris角点算法

特征点检测广泛应用到目标匹配、目标跟踪、三维重建等应用中，在进行目标建模时会对图像进行目标特征的提取，常用的有颜色、角点、特征点、轮廓、纹理等特征。现在开始讲解常用的特征点检测，其中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};或者 //定义垂直方向差分算子并求Iydouble 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;elseCV_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版本

里面函数最基础的代码解释和注释

   C源码 

2、基于opencv源码

   opencv+C++推荐 

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，程序中通过图像的膨胀运算来达到检测极大值的目的，因为默认参数的膨胀运算就是用窗口内的最大值替代当前的灰度值。

   View Code 

四、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 )