2、Harris角点检测算法

最新推荐文章于 2020-07-17 21:02:51 发布

诺亚方舟369

最新推荐文章于 2020-07-17 21:02:51 发布

阅读量3.4k

点赞数 10

分类专栏：图像处理与计算机视觉文章标签：图像处理计算机视觉角点检测

本文链接：https://blog.csdn.net/a1056258349/article/details/81220400

版权

图像处理与计算机视觉专栏收录该内容

2 篇文章 0 订阅

订阅专栏

1、基本原理

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

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

图1

对于图像I(x,y)，当在点(x,y)处平移(u,v)后的自相似性，可以通过自相关函数给出：

$E(u,v)=\sum_{x,y}\omega (x,y)[I(x+u,y+v)-I(x,y)]^2$ (1.1)

其中，w(x,y)是以点(x,y)为中心的窗口加权函数，它既可是常数（图1左），也可以是高斯加权函数（图1右）。高斯函数的表达式如下：

$\omega (x,y)=\frac{1}{2\pi \sigma ^2}e^{-(x^2+y^2)/2\sigma ^2}$ (1.2)

（2）Moravec算子只考虑每隔45度方向，Harris用Taylor展开去近似任意方向；

利用泰勒级数将I(x+u,y+v)展开得：

$I(x+u,y+v)=I(x,y)+I_xu+I_yv+o(u^2,v^2)$ (1.3)

则E(u,v)表达式可以更新为：

$E(u,v)\cong \sum_{x,y}\omega (x,y)[I_xu+I_yv]^2$ (1.4)

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

$E(u,v)\cong [u,v]M[u,v]^T$ (1.5)

矩阵M为实对称矩阵，且

$M=\omega (x,y)\bigotimes \begin{bmatrix} I_x^2 & I_xI_y\\ I_xI_y& I_y^2 \end{bmatrix}=\begin{bmatrix} A & C\\ C& B \end{bmatrix}$ (1.6)

则E(u,v)化简为二次项函数：

$E(u,v)=Au^2+Bv^2+2Cuv$ (1.7)

二次项函数本质上就是一个椭圆函数。椭圆的扁率和尺寸是由M的特征值 $\lambda _1,\lambda _2$ 决定的，椭圆的方向是由M的特征向量决定的，如图2所示。

图2

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

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

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

3）图像中的角点。两个特征值都比较大，且近似相等，自相关函数在所有方向上都增大。

图3

（3）角点响应函数

根据二次项函数特征值的计算公式，我们可以求M矩阵的特征值。但是Harris给出的角点差别方法并不需要计算具体的特征值，而是计算一个角点响应值R来判断角点。R的计算公式为：

$R=detM-k(traceM)^2$ (1.8)

其中，detM为矩阵M的行列式，traceM为矩阵M的迹，k为常数，取值范围为0.04~0.06。事实上，特征是隐含在detM和traceM中，因为

$detM=\lambda _1\lambda _2=AB-C^2$ (1.9)

$traceM=\lambda _1+\lambda _2=A+B$ (1.10)

其实，角点量R的计算方式可以自由发挥，只要能反应角点的特征即可。例如，Nobel于1988年提出利用如下公式计算角点的响应值，无需设定参数k的值：

$cim=\frac{I_x^2I_y^2-(I_xI_y)^2}{I_x^2+I_y^2}$ (1.11)

采用上述公式计算角点的CRF值，从而避免的参数k对角点选取的影响，在实际应用中，通常选用这个改进的Harris角点检测算法进行检测：当cim值大于预定的阈值，则该点为角点候选点，通过非极大值抑制挑选出最终的角点。

此外，或如式（1.12）计算R值，也无需考虑参数k。

$R=\frac{detM}{traceM}$ (1.12)

（4）最后设定R的阈值进行角点判断，以及角点的极大值抑制等。

二、算法步骤

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

输出：角点检测图

具体步骤：

（1）计算图像I(x,y)在X和Y方向的梯度Ix和Iy；

$I_x=\frac{\partial I}{\partial x}=I\bigotimes (-1,0,1)$

$I_y=\frac{\partial I}{\partial x}=I\bigotimes (-1,0,1)^T$

（2）计算梯度方向的乘积；

$I_x^2=I_x\cdot I_x$

$I_y^2=I_y\cdot I_y$

$I_x_y=I_x\cdot I_y$

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

$A=g(I_x^2)=I_x^2\bigotimes \omega$

$B=g(I_y^2)=I_y^2\bigotimes \omega$

$C=g(I_x_y)=I_x_y\bigotimes \omega$

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

$R=(AB-C^2)-k(A+B)^2$

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

算法结束。

三、优缺点

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

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

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

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

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

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

如图4所示，当图像被缩小时，在检测窗口尺寸不变的前提下，窗口内所包含图像的内容可能是完全不同的。左侧的图像可能被检测为边缘或曲线，而右侧的图像则可能被检测为一个角点。或者说如果图像尺度发生变化，原来是角点的点在新的尺度可能就不是角点了。

图 4

注：尺度不变性问题可通过图像金字塔解决，例如，在运算的开始先将图像转化到尺度空间表示，即将原图像进行尺度变换，而尺度变换的方式就是源图像与尺度核函数做卷积运算：

$L(x,y,\sigma )=g(X,Y,\sigma )\bigotimes I(x,y)$

其中，sigma表示尺度。然后使用L代替原图像去进行运算，尺度为运算的参数。关于图像金字塔的具体内容回继续学习！

Harris角点本身就不受光照，旋转的影响，现在又使其满足尺度不变性，至此，Harris角点可以成为一个优秀的特征了。

四、代码分析

/**********************************************************************************
*函数 Mat detectHarrisCorners(const Mat& imgSrc, double alpha)
*输入:
      *imgSrc : 源单通道图像
      *alpha  : Harris响应函数参数
*输出
     *imgDst : 提取到角点的图像
***************************************************************************************/

#include <iostream> 
#include <vector>
#include "opencv2/opencv.hpp"
#include "opencv2/imgproc/imgproc.hpp"  
#include "opencv2/highgui/highgui.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=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;//响应函数值R
		}
	}
	
	double maxStrength;
	minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);
	double qualityLevel = 0.1;
	double thresh = qualityLevel * maxStrength;// 设置threshold
 
	Mat 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 gray1 = imgSrc.clone();
	vector<Point>::const_iterator it = points.begin();
	while (it != points.end())
	{
		circle(gray1, *it, 3, Scalar(255, 255, 255), 1);
		++it;
	}

	return gray1;
}

int main()
{
	Mat SrcImage = imread("lena.png");
	if (!SrcImage.data)
		return -1;
    Mat gray;
	if (SrcImage.channels() == 3)
	{
		cvtColor(SrcImage, gray, CV_BGR2GRAY);  //RGB转换成灰度图像
	}
	else
	{
		gray = SrcImage.clone();
	}

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

实验结果图如下，代码是拼凑的并未优化，运行时间还是比较长的。阈值设置的比较大，检测到的角点数量也比较少。

图5

诺亚方舟369

关注

10
点赞
踩
26

收藏

觉得还不错? 一键收藏
2
评论
2、Harris角点检测算法

1、基本原理Harris算子是对Moravec算子的改进，包括：（1）Harris算子用高斯函数代替Moravec算子的二值窗口函数，如图1所示，窗口函数应对离中心点越的像素赋予越大的权重，以减少噪声影响；图1对于图像I(x,y)，当在点(x,y)处平移(u,v)后的自相似性，可以通过自相关函数给出： (1.1)其中，w(x,y)是以点(x,y)为中心的窗口加权函数，它...
复制链接

扫一扫