opencv harris 角点检测

一、图像特征的分类

【OpenCV入门教程之十六】OpenCV角点检测之Harris角点检测_【浅墨的游戏编程Blog】毛星云（浅墨）的专栏-CSDN博客_基于轮廓曲线的角点检测

下面部分图像来自与B站上将opencv的一位叫贾志刚讲师课件的截图

1、边缘

2、角点（兴趣点）：如果某一点在任意方向的一个微小变动都会引起灰度很大的变化，那么我们就把它称之为角点

角点位于两条边缘的交点处，代表了两个边缘变化的方向上的点,，所以他们是可以精确定位的二维特征，甚至可以达到亚像素的精度，且其图像梯度有很高的变化，这种变化是可以用来帮助检测角点的

3、斑点（blobs）

二、角点检测算法

思想：

使用一个固定窗口在图像上进行任意方向上的滑动，比较滑动前与滑动后两种情况，窗口中的像素灰度变化程度，如果存在任意方向上的滑动，都有着较大灰度变化，那么我们可以认为该窗口中存在角点。

在当前的图像处理领域，角点检测算法可归纳为三类：

• 基于灰度图像的角点检测：基于梯度、基于模板和基于模板梯度组合
• 基于二值图像的角点检测 ：Harris角点检测算法、KLT角点检测算法及SUSAN角点检测算法
• 基于轮廓曲线的角点检测

• 关于角点的具体描述可以有几种：

  • 一阶导数(即灰度的梯度)的局部最大所对应的像素点；
  • 两条及两条以上边缘的交点；
  • 图像中梯度值和梯度方向的变化速率都很高的点；
  • 角点处的一阶导数最大，二阶导数为零，指示物体边缘变化不连续的方向。

三、Harris角点检测原理

从上面可以知道什么是角点了，那么我们先移动窗口，对其描述如下：

w(x,y)是窗口函数，最简单情形就是窗口W内的所有像素所对应的w权重系数均为1.但有时候，我们会将w(x,y)函数设置为以窗口W中心为原点的二元正太分布。如果窗口W中心点是角点时，移动前与移动后，该点在灰度变化贡献最大；而离窗口W中心(角点)较远的点，这些点的灰度变化几近平缓，这些点的权重系数，可以设定小值，以示该点对灰度变化贡献较小，那么我们自然而然想到使用二元高斯函数来表示窗口函数。

接着我们要复习一下数学分析关于泰勒级数的内容，就像任何连续周期的信号都可以由一组适当的正弦曲线组合而成一样，任何一个表达式都可以用泰勒公式来无限逼近：

 那么我们将泰勒公式带入到上面计算E（u，v）中：

 接下啦就是如何求解M这个矩阵

难道我们是直接求上述的E(u,v)值来判断角点吗？Harris角点检测并没有这样做，而是通过对窗口内的每个像素的x方向上的梯度与y方向上的梯度进行统计分析。这里以Ix和Iy为坐标轴，因此每个像素的梯度坐标可以表示成(Ix,Iy)。针对平坦区域，边缘区域以及角点区域三种情形进行分析：

 请注意看下面的这幅图

注意到这三种区域的特点:

• 平坦区域上的每个像素点所对应的(Ix,Iy)坐标分布在原点附近，其实也很好理解，针对平坦区域的像素点，他们的梯度方向虽然各异，但是其幅值都不是很大，所以均聚集在原点附近；
• 边缘区域有一坐标轴分布较散，至于是哪一个坐标上的数据分布较散不能一概而论，这要视边缘在图像上的具体位置而定，如果边缘是水平或者垂直方向，那么Iy轴方向或者Ix方向上的数据分布就比较散；
• 角点区域的x、y方向上的梯度分布都比较散。

我们是不是可以根据这些特征来判断哪些区域存在角点呢

 

 比较两个特征值的关系：不能做减法

 

结论：

做减法是有问题的，那么我们如何来用λ1λ2来表示差值又能体现上面说的一些结论呢？

 定义一个新的函数R  

 det(M)=λ1λ2， 行列式

trace(M)=λ1+λ2  迹；

K是一个比较值，是一个我们给出的定值

 R取决于MM的特征值，

对于角点|R|很大，

平坦的区域|R|很小

边缘的R为负值;

• 特征值都比较大时，即窗口中含有角点；
• 特征值一个较大，一个较小，窗口中含有边缘；
• 特征值都比较小，窗口处在平坦区域；

四、Harris角点检测步骤

1、计算图像 x、y方向上的梯度Ix Iy

2、计算Ix* Iy

3、使用高斯函数对I2x、I2y、IxIy进行高斯加权(取σ=2,ksize=3)，计算中心点为(x,y)的窗口WW对应的矩阵M；

 4、计算R

5、过滤大于t的R值 

五、 cornerHarris

void cornerHarris(
 InputArray src, // 需要为单通道8位或者浮点型图像
 OutputArray dst,  // Harris角点检测的输出结果，和原图片有一样的尺寸和类型
 int block Size, //表示邻域的大小   cornerEigenValsAndVecs()中讲到
 int ksize,  //表示Sobel()算子的孔径的大小
 double k,  //计算响应公式中的k值，一般取0.04~0.06；
int borderType = BORDER_DEFAULT //图像像素的边界模式。注意它有默认值BORDER_DEFAULT
)

void cv::cornerHarris( InputArray _src,OutputArray _dst, int blockSize, int ksize, double k, int borderType )
{
   Mat src = _src.getMat();
   _dst.create( src.size(), CV_32F );
   Mat dst = _dst.getMat();
   cornerEigenValsVecs( src, dst, blockSize, ksize, HARRIS, k, borderType);
}

tatic void
cornerEigenValsVecs( const Mat& src,Mat& eigenv, int block_size,
                     int aperture_size, intop_type, double k=0.,
                     intborderType=BORDER_DEFAULT )
{
#ifdef HAVE_TEGRA_OPTIMIZATION
   if (tegra::cornerEigenValsVecs(src, eigenv, block_size, aperture_size,op_type, k, borderType))
       return;
#endif
 
   int depth = src.depth();
   double scale = (double)(1 << ((aperture_size > 0 ?aperture_size : 3) - 1)) * block_size;
   if( aperture_size < 0 )
       scale *= 2.;
   if( depth == CV_8U )
       scale *= 255.;
   scale = 1./scale;
 
   CV_Assert( src.type() == CV_8UC1 || src.type() == CV_32FC1 );
 
   Mat Dx, Dy;
   if( aperture_size > 0 )
    {
       Sobel( src, Dx, CV_32F, 1, 0, aperture_size, scale, 0, borderType );
       Sobel( src, Dy, CV_32F, 0, 1, aperture_size, scale, 0, borderType );
    }
   else
    {
       Scharr( src, Dx, CV_32F, 1, 0, scale, 0, borderType );
       Scharr( src, Dy, CV_32F, 0, 1, scale, 0, borderType );
    }
 
   Size size = src.size();
   Mat cov( size, CV_32FC3 );
   int i, j;
 
   for( i = 0; i < size.height; i++ )
    {
       float* cov_data = (float*)(cov.data + i*cov.step);
       const float* dxdata = (const float*)(Dx.data + i*Dx.step);
       const float* dydata = (const float*)(Dy.data + i*Dy.step);
 
        for( j = 0; j < size.width; j++ )
       {
           float dx = dxdata[j];
           float dy = dydata[j];
 
           cov_data[j*3] = dx*dx;
           cov_data[j*3+1] = dx*dy;
           cov_data[j*3+2] = dy*dy;
       }
    }
 
   boxFilter(cov, cov, cov.depth(), Size(block_size, block_size),
       Point(-1,-1), false, borderType );
 
   if( op_type == MINEIGENVAL )
       calcMinEigenVal( cov, eigenv );
   else if( op_type == HARRIS )
       calcHarris( cov, eigenv, k );
   else if( op_type == EIGENVALSVECS )
       calcEigenValsVecs( cov, eigenv );
}
 
}

六、代码演示

#if  1  // harris  角点检测
const  char* INPUT_TITLE = "原图";
const  char* OUT_TITLE = "原图显示角点";
const  char* OUT_TITLE2 = "灰度图显示角点";

Mat src, src1, gray;
int g_nThresh = 30; //当前阀值
int g_nMaxThresh = 175;

void callbackCornerHarris(int, void*)
{
	Mat dstImage; //目标图
	Mat normImage; //归一化后的图
	Mat scaledImage; //线性变换后的八位无符号整形的图

	// 置零当前需要显示的两幅图，即清除上一次调用此函数时他们的值
	dstImage = Mat::zeros(src.size(), CV_32FC1);
	src1 = src.clone();

	//角点检测
	cornerHarris(gray, dstImage, 2, 3, 0.04, BORDER_DEFAULT);
	normalize(dstImage, normImage, 0, 255, NORM_MINMAX, CV_32FC1, Mat());

	//将归一化后的图线性变换成8位无符号整形
	convertScaleAbs(normImage, scaledImage);

	//将检测到的，且符合阀值条件的角点绘制出来
	for (int i = 0; i < normImage.rows; i++)
	{
		for (int j = 0; j < normImage.cols; j++)
		{
			if ((int)normImage.at<float>(i, j) > g_nThresh + 80)
			{
				circle(src1, Point(j, i), 5, Scalar(10, 10, 255), 2, 8, 0);
				circle(scaledImage, Point(j, i), 5, Scalar(0, 10, 255), 2, 8, 0);
			}
		}
	}
	imshow(OUT_TITLE, src1);
	imshow(OUT_TITLE2, scaledImage);
}

int main()
{
	src = imread("C:\\Users\\19473\\Desktop\\opencv_images\\corners.png");
	if (!src.data)
	{
		printf("could not  load  image....\n");
	}
	imshow(INPUT_TITLE, src);
	src.copyTo(src1);
	cvtColor(src1, gray, COLOR_BGR2GRAY);

	namedWindow(OUT_TITLE, WINDOW_AUTOSIZE);
	createTrackbar("阀值：", OUT_TITLE, &g_nThresh, g_nMaxThresh, callbackCornerHarris);

	callbackCornerHarris(0, NULL);

	waitKey(0);
	return 0;
}

#endif