Harris角点检测算法-CSDN博客

本文的内容基本来源于以下博客，主要是为了方便以后再次学习。

Harris原理分析：

人眼对角点的识别通常是在一个局部的小区域或小窗口完成的。如果这个特定的小窗口在各个方向上移动时，窗口内图像的灰度发生了较大的变化，那么就认为在窗口内遇到了角点。如果这个特定的窗口在图像各个方向上移动时，窗口内图像的灰度没有发生变化，那么窗口内就不存在角点；如果窗口在某一个方向移动时，窗口内图像的灰度发生了较大的变化，而在另一些方向上没有发生变化，那么，窗口内的图像可能就是一条线段。如下图：

首先，对于图像，在点(x, y)处平移(u, v)后产生的灰度变化的自相关函数如下：

其中，窗口函数可以是平坦的，也可以是高斯的，如下图：

将 I(x+u, y+v) 进行泰勒展开如下：

则有

由于是局部微小的移动量(u, v)，所以可以近似得到下面忽略余项之后的表达式，是一个二项式函数：

所以，

其中，M的表达式如下，可由图像的导数求得：

上面说到，忽略余项之后的表达式为一个二项式函数。本质上，二项式函数就是一个椭圆函数，椭圆的扁率和尺寸是由M(x,y)的特征值λ1、λ2决定的，椭圆的方向是由M(x,y)的特征矢量决定的，如下图所示，椭圆方程为：

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

a. 图像中的直线：一个特征值大，另一个特征值小，λ1>λ2或λ2>λ1；自相关函数值在某一方向上大，在其他方向上小。

b. 图像中的平面：两个特征值都小，且近似相等；自相关函数数值在各个方向上都小。

c. 图像中的角点：两个特征值都大，且近似相等；自相关函数在所有方向都增大。

我们是通过M的两个特征值的大小对图像进行分类，所以，定义角点相应函数R：

其中k为经验常数，一般是一个远小于1的数，在opencv中的默认值是0.04。

对于边缘部分，λ1>>λ2且 λ2< k*λ1，角点相应函数R: λ1λ2 - k(λ1+λ2)^2 ≈ λ1λ2 - k(λ1)^2 < λ1 k λ1 - k(λ1)^2 = 0，因此，在边缘处的角点响应值小于0。

对于非边缘部分，λ1、λ2相差不大，可以认为(λ1+λ2)^2 ≈ 4λ1λ2，角点相应函数R: λ1λ2 - k(λ1+λ2)^2 ≈ λ1λ2- 4kλ1λ2 = (1- 4k)λ1λ2，由于k远小于1，因此， (1- 4k) ≈ 1，所以λ1λ2 - k(λ1+λ2)^2 ≈ λ1λ2

在角点处，λ1、λ2都比较大，故角点响应值也比较大；在平坦区域λ1、λ2都很小，故角点响应值也很小。

所以，上图可以转化为：

其中：

• R 只与M的特征值有关

• 角点：R 为大数值正数

• 边缘：R 为大数值负数

• 平坦区：R 为小数值

在判断角点的时候，对角点响应函数R进行阈值处理：R > threshold，并且提取R的局部极大值。

Harris角点的缺点：

由Harris角点的定义可知，该角点具有旋转不变性。

Harris角点不具有尺度不变性。

Harris代码实现：

void detectHarrisCorners(const Mat& imgSrc, Mat& imgDst, double alpha)
{
	Mat gray;
	if (imgSrc.channels() == 3)
	{
		cvtColor(imgSrc, gray, CV_BGR2GRAY);
	}
	else
	{
		gray = imgSrc.clone();
	}
	gray.convertTo(gray, CV_64F);

	Mat xKernel = (Mat_<double>(1, 3) << -1, 0, 1);
	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(7, 1);
	filter2D(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;
		}
	}
	
	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);


	Mat cornerMap;
	double qualityLevel = 0.01;
	// 根据角点响应最大值计算阈值
	double thresh = qualityLevel * maxStrength;
	// 根据阈值得到角点图
	cornerMap = cornerStrength > thresh;
	// 和局部最大值图进行与操作，完成非最大值抑制
	bitwise_and(cornerMap, localMax, cornerMap);

	imgDst = cornerMap.clone();

}

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);
	}
}