原理

Harris算子是对Moravec算子的改进，包括：

（1）Harris算子用高斯函数代替Moravec算子的二值窗口函数，如下图所示，窗口函数应对离中心点越的像素赋予越大的权重，以减少噪声影响；

图1
对于图像I(x,y)，当在点(x,y)处平移(u,v)后的自相似性，可以通过自相关函数给出：
在这里插入图片描述
其中，w(x,y)是以点(x,y)为中心的窗口加权函数，它既可是常数（图1左），也可以是高斯加权函数（图1右）。高斯函数的表达式如下：

（2）Moravec算子只考虑每隔45度方向，Harris用Taylor展开去近似任意方向；
利用泰勒级数将I(x+u,y+v)展开得：
在这里插入图片描述
则E(u,v)表达式可以更新为：

其中，Ix为x方向的差分（一阶微分近似），Iy为y方向的差分，w(x,y)为高斯函数。通过推导，E(u,v)以矩阵的形式表示为：

矩阵M为实对称矩阵，且

则E(u,v)化简为二次项函数：
在这里插入图片描述
二次项函数本质上就是一个椭圆函数。椭圆的扁率和尺寸是由M的特征值\lambda _1,\lambda _2决定的，椭圆的方向是由M的特征向量决定的，如图2所示。

图2
椭圆函数特征值与图像中的角点、直线（边缘）和平面之间的关系如图3所示。共可分为三种情况：

1）图像上的直线。一个特征值比较大，另外一个特征值比较小；自相关函数在某一个方向上大，在另一个方向上小。

2）图像中的平滑区域。两个特征值都小，且近似相等；自相关函数值子各个方向上都小。

3）图像中的角点。两个特征值都比较大，且近似相等，自相关函数在所有方向上都增大。
在这里插入图片描述
图3
（3）角点响应函数
根据二次项函数特征值的计算公式，我们可以求M矩阵的特征值。但是Harris给出的角点差别方法并不需要计算具体的特征值，而是计算一个角点响应值R来判断角点。R的计算公式为：
在这里插入图片描述
其中，detM为矩阵M的行列式，traceM为矩阵M的迹，k为常数，取值范围为0.04~0.06。事实上，特征是隐含在detM和traceM中，因为

其实，角点量R的计算方式可以自由发挥，只要能反应角点的特征即可。例如，Nobel于1988年提出利用如下公式计算角点的响应值，无需设定参数k的值：
在这里插入图片描述
采用上述公式计算角点的CRF值，从而避免的参数k对角点选取的影响，在实际应用中，通常选用这个改进的Harris角点检测算法进行检测：当cim值大于预定的阈值，则该点为角点候选点，通过非极大值抑制挑选出最终的角点。

此外，或如式（1.12）计算R值，也无需考虑参数k。
在这里插入图片描述
4）最后设定R的阈值进行角点判断，以及角点的极大值抑制等。

算法步骤

输入：源单通道图像，参数k

输出：角点检测图

具体步骤：

（1）计算图像I(x,y)在X和Y方向的梯度Ix和Iy；
在这里插入图片描述
（2）计算梯度方向的乘积；

（3）使用高斯核进行加权，计算矩阵M的元素A，B，C；

（4）计算角点响应函数R，并设定阈值，当R小于阈值时，不是候选角点；

（5）进行局部极大值抑制。

算法结束。

优缺点

Harris角点检测算法有诸多优点，但也有不完善的地方。

（1）Harris角点检测算子具有旋转不变性

Harris角点检测算子使用的是角点附近的区域灰度二阶矩矩阵。而二阶矩矩阵可以表示成一个椭圆，椭圆的长短轴正是二阶矩矩阵特征值平方根的倒数。当特征椭圆转动时，特征值并不发生变化，所以判断角点响应值R也不发生变化，由此说明Harris角点检测算子具有旋转不变性。

（2）Harris角点检测算子对灰度平移和灰度尺度变化不敏感

这是因为在进行Harris角点检测时，使用了微分算子对图像进行微分运算，而微分运算对图像亮度的抬高或下降（I=I+a）、密度的拉升或收缩（I=bI）不敏感。换言之，对亮度和对比度的仿射变换并不改变Harris响应的极值点出现的位置，但是，由于阈值的选择，可能会影响角点检测的数量。

（3）Harris角点检测算子不具有尺度不变性

如图4所示，当图像被缩小时，在检测窗口尺寸不变的前提下，窗口内所包含图像的内容可能是完全不同的。左侧的图像可能被检测为边缘或曲线，而右侧的图像则可能被检测为一个角点。或者说如果图像尺度发生变化，原来是角点的点在新的尺度可能就不是角点了。
在这里插入图片描述
注：尺度不变性问题可通过图像金字塔解决，例如，在运算的开始先将图像转化到尺度空间表示，即将原图像进行尺度变换，而尺度变换的方式就是源图像与尺度核函数做卷积运算：

其中，sigma表示尺度。然后使用L代替原图像去进行运算，尺度为运算的参数。
Harris角点本身就不受光照，旋转的影响，现在又使其满足尺度不变性，至此，Harris角点可以成为一个优秀的特征了。

源码

/**********************************************************************************
*函数 Mat detectHarrisCorners(const Mat& imgSrc, double alpha)
*输入:
*imgSrc : 源单通道图像
*alpha  : Harris响应函数参数
*输出
*imgDst : 提取到角点的图像
***************************************************************************************/#include <iostream> 
#include "opencv2/opencv.hpp"using namespace cv;
using namespace std;Mat detectHarrisCorners(const Mat& imgSrc, double alpha)
{Mat gray;gray = imgSrc.clone();gray.convertTo(gray, CV_64F);Mat xKernel = (Mat_<double>(1, 3) << -1, 0, 1);//水平方向模板计算Ix Mat yKernel = xKernel.t();Mat Ix, Iy;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);Mat gaussKernel = getGaussianKernel(5, 1);//获得高斯核size=5，sigma=1filter2D(Ix2, Ix2, CV_64F, gaussKernel);filter2D(Iy2, Iy2, CV_64F, gaussKernel);filter2D(Ixy, Ixy, CV_64F, gaussKernel);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;//响应函数值R}}double maxStrength;minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);double qualityLevel = 0.1;double thresh = qualityLevel * maxStrength;// 设置thresholdMat dilated, localMax;//默认3 * 3核膨胀，膨胀之后，除了局部最大值点和原来相同，其它非局部最大值点被3*3邻域内的最大值点取代dilate(cornerStrength, dilated, Mat());//与原图相比，只剩下和原图值相同的点，这些点都是局部最大值点，保存到localMax compare(cornerStrength, dilated, localMax, CMP_EQ);//和局部最大值图与，剩下角点局部最大值图，即：完成非最大值抑制Mat cornerMap;cornerMap = cornerStrength > thresh;bitwise_and(cornerMap, localMax, cornerMap);vector<Point> points;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(Point(x, y));}}//画角点Mat imgRGB = imgSrc.clone();cvtColor(imgRGB, imgRGB, COLOR_GRAY2RGB);vector<Point>::const_iterator it = points.begin();while (it != points.end()){circle(imgRGB, *it, 3, Scalar(0, 0, 255), 1);++it;}return imgRGB;
}int main()
{Mat gray = imread("D:\\Work\\OpenCV\\Workplace\\Test_1\\3.jpg", 0);if (!gray.data)return -1;imshow("grayImage", gray);double time0 = static_cast<double>(getTickCount());Mat HarrisImage = detectHarrisCorners(gray, 0.05);time0 = ((double)getTickCount() - time0) / getTickFrequency();cout << "runtime ：" << time0 << "s" << endl;imshow("HarrisImage", HarrisImage);waitKey(0);return 0;
}