＜图像处理＞ Harris角点检测

Harris角点检测原理

Harris角点检测是一种计算机视觉中常用的角点检测算法，用于在图像中检测出角点特征。角点通常被定义为两条边的交点，或者说，角点的局部邻域应该具有两个不同区域的不同方向的边界。Harris角点检测算法是最常用且最基础的角点检测算法之一。

Harris角点检测算法的原理是通过计算图像中每个像素点的响应函数来判断其是否为角点。该响应函数基于图像的灰度变化和局部窗口的协方差矩阵，通过计算特征值来确定角点的位置。当特征值较大时，表示该像素点处存在角点。

优点： 对图像旋转、尺度变化和亮度变化具有一定的不变性，且计算简单快速。

当一幅图像$f$，检测窗口在图像上滑动，当检测窗口在点$(x,y)$处位移$(\Delta x,\Delta y)$后，它的自相关函数为：
$$C(x,y)=\sum _{x_i\in w}\sum _{y_i\in w}(f(x_i,y_i)-f(x_i-\Delta x,y_i-\Delta y){)}^2$$
角点不会受光圈问题的影响，对于所有$\Delta x,\Delta y，S_w(\Delta x,\Delta y)$都是高响应，如果平移图像用一阶泰勒展开近似，则可表示为
$$f(x_i-\Delta x,y_i-\Delta y)\approx (f(x_i,y_i)+\left[\begin{array}{cc}\frac{\partial f(x_i,y_i)}{\partial x}& \frac{\partial f(x_i,y_i)}{\partial y}\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta y\end{array}\right]$$
带入：
$$C(x,y)=\sum _{x_i\in w}\sum _{y_i\in w}(f(x_i,y_i)-f(x_i,y_i)-\left[\begin{array}{cc}\frac{\partial f(x_i,y_i)}{\partial x}& \frac{\partial f(x_i,y_i)}{\partial y}\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta y\end{array}\right]{)}^2$$
$$=\sum _{x_i\in w}\sum _{y_i\in w}(-\left[\begin{array}{cc}\frac{\partial f(x_i,y_i)}{\partial x}& \frac{\partial f(x_i,y_i)}{\partial y}\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta y\end{array}\right]{)}^2$$
$$=\sum _{x_i\in w}\sum _{y_i\in w}(\left[\begin{array}{cc}\frac{\partial f(x_i,y_i)}{\partial x}& \frac{\partial f(x_i,y_i)}{\partial y}\end{array}\right]\left[\begin{array}{c}\Delta x\\ \Delta y\end{array}\right]{)}^2$$
因为$u^2=u^{T}u$
$$=\sum _{x_i\in w}\sum _{y_i\in w}[\Delta x,\Delta y]\left(\left[\begin{array}{c}\frac{\partial f}{\partial x}\\ \frac{\partial f}{\partial y}\end{array}\right]\left[\begin{array}{c}\frac{\partial f}{\partial x}\frac{\partial f}{\partial y}\end{array}\right]\right)\left[\begin{array}{c}\Delta x\\ \Delta y\end{array}\right]$$
$$=[\Delta x,\Delta y]\left(\sum _{x_i\in w}\sum _{y_i\in w}\left[\begin{array}{c}\frac{\partial f}{\partial x}\\ \frac{\partial f}{\partial y}\end{array}\right]\left[\begin{array}{c}\frac{\partial f}{\partial x}\frac{\partial f}{\partial y}\end{array}\right]\right)\left[\begin{array}{c}\Delta x\\ \Delta y\end{array}\right]$$
$$=[\Delta x,\Delta y]A_w(x,y)\left[\begin{array}{c}\Delta x\\ \Delta y\end{array}\right]$$
$A_w(x,y)$表示图像区域$w$在点$(x,y)=(0,0)$处的二阶导数的一半。$A$为：
$$A(x,y)=2\bullet \left[\begin{array}{cc}{\sum }_{x_i\in w}{\sum }_{y_i\in w}(\frac{\partial f(x_i,y_i)}{\partial x}{)}^2& {\sum }_{x_i\in w}{\sum }_{y_i\in w}\frac{\partial f(x_i,y_i)}{\partial x}\frac{\partial f(x_i,y_i)}{\partial y}\\ {\sum }_{x_i\in w}{\sum }_{y_i\in w}\frac{\partial f(x_i,y_i)}{\partial x}\frac{\partial f(x_i,y_i)}{\partial y}& {\sum }_{x_i\in w}{\sum }_{y_i\in w}(\frac{\partial f(x_i,y_i)}{\partial y}{)}^2\end{array}\right]$$

局部结构矩阵$A$代表邻域的变化，其有两个特征值${\lambda }_1，{\lambda }_2$，有以下三种情形出现：
1、两个特征值都很小，这意味着，图像f在检测点处平坦——在该点没有边缘或角点；
2、一个特征值很小，另一个很大，局部邻域呈脊状，如有垂直于脊部的微小移动，图像f将发生显著的变化；
3、两个特征值都很大，任何方向的移动，都会造成图像f的显著变化，此时，检测到一个角点。

对于求解特征值需要一元二次方程，这必然会降低运算效率。因此，通过计算响应函数$R(A)=det(A)-k(trace(A){)}^2$，即$R={\lambda }_1{\lambda }_2-k({\lambda }_1+{\lambda }_2{)}^2$，
当$A(x,y)=\left[\begin{array}{cc}a& b\\ b& c\end{array}\right]$
$$det(A)={\lambda }_1{\lambda }_2=ac-b^2$$
$$trace(A)={\lambda }_1+{\lambda }_2=a+c$$
所以，
$$R=ac-b^2-k(a+c{)}^2$$

检测步骤：

1. 计算图像$I$的梯度，获取水平和垂直方向上的梯度$I_x$和$I_y$；
2. 构建协方差矩阵，计算每个像素点的协方差矩阵，即计算两个方向梯度的乘积$I_x^2=I_x\ast I_x，I_y^2=I_y\ast I_y，I_{xy}=I_x\ast I_y$；
3. 使用加权窗口$w$（可以是高斯加权，可以是平均加权，OpenCV中用的平均加权），对$I_x^2,I_y^2,I_{xy}$进行加权处理，得到自相关元素$a=I_x^2\ast w，b=I_{xy}^2\ast w，c=I_y^2\ast w$；
4. 计算协方差矩阵的特征值$R$，判断每个像素点是否为角点；
5. 根据设定的阈值，筛选出最终的角点。

OpenCV函数

void cv::cornerHarris(InputArray src, OutputArray dst, int blockSize, int ksize, double k, int borderType=BORDER_DEFAULT)

Parameters
src				输入图像，单通道8位或浮点型；
dst				输出Harris算子响应图像，该变量存储的R值，类型为单通道32位浮点型；
blockSize		步骤3中的加权窗口尺寸；
ksize			Sobel算子求梯度的尺寸；
k				响应函数中的k；
borderType		填充边缘方法，见边缘填充博客；

OpenCV源码分析

cornerHarris函数在sources\modules\imgproc\src\corner.cpp文件内定义

void cv::cornerHarris( InputArray _src, OutputArray _dst, int blockSize, int ksize, double k, int borderType )
{
    CV_INSTRUMENT_REGION();

    CV_OCL_RUN(_src.dims() <= 2 && _dst.isUMat(),
               ocl_cornerMinEigenValVecs(_src, _dst, blockSize, ksize, k, borderType, HARRIS))

    Mat src = _src.getMat();
    _dst.create( src.size(), CV_32FC1 );
    Mat dst = _dst.getMat();

#ifdef HAVE_IPP
    int borderTypeNI = borderType & ~BORDER_ISOLATED;
    bool isolated = (borderType & BORDER_ISOLATED) != 0;
#endif
    CV_IPP_RUN(((ksize == 3 || ksize == 5) && (_src.type() == CV_8UC1 || _src.type() == CV_32FC1) &&
        (borderTypeNI == BORDER_CONSTANT || borderTypeNI == BORDER_REPLICATE) && CV_MAT_CN(_src.type()) == 1 &&
        (!_src.isSubmatrix() || isolated)) && IPP_VERSION_X100 >= 810, ipp_cornerHarris( src, dst, blockSize, ksize, k, borderType ));

	//调用计算角点特征值的函数，其中参数HARRIS表示该函数用于HARRIS算子
    cornerEigenValsVecs( src, dst, blockSize, ksize, HARRIS, k, borderType );
}

cornerEigenValsVecs函数：

static void
cornerEigenValsVecs( const Mat& src, Mat& eigenv, int block_size,
                     int aperture_size, int op_type, double k=0.,
                     int borderType=BORDER_DEFAULT )
{
#if CV_TRY_AVX
    bool haveAvx = CV_CPU_HAS_SUPPORT_AVX;
#endif

    int depth = src.depth();
    //计算尺度，该值在梯度计算的时候使用，主要目的是通过减小尺度，来提高平滑处理的速度
    double scale = (double)(1 << ((aperture_size > 0 ? aperture_size : 3) - 1)) * block_size;
    if( aperture_size < 0 )
        scale *= 2.0;
    if( depth == CV_8U )
        scale *= 255.0;
    scale = 1.0/scale;

    CV_Assert( src.type() == CV_8UC1 || src.type() == CV_32FC1 );
	
	//Dx：水平方向梯度；Dy：垂直方向梯度
    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;

	//遍历每个像素，计算Ix^2,Iy^2,Ixy
    for( i = 0; i < size.height; i++ )
    {
        float* cov_data = cov.ptr<float>(i);
        const float* dxdata = Dx.ptr<float>(i);
        const float* dydata = Dy.ptr<float>(i);

#if CV_TRY_AVX
        if( haveAvx )
            j = cornerEigenValsVecsLine_AVX(dxdata, dydata, cov_data, size.width);
        else
#endif // CV_TRY_AVX
            j = 0;

//指令集加速
#if CV_SIMD128
        {
            for( ; j <= size.width - v_float32x4::nlanes; j += v_float32x4::nlanes )
            {
                v_float32x4 v_dx = v_load(dxdata + j);
                v_float32x4 v_dy = v_load(dydata + j);

                v_float32x4 v_dst0, v_dst1, v_dst2;
                v_dst0 = v_dx * v_dx;
                v_dst1 = v_dx * v_dy;
                v_dst2 = v_dy * v_dy;

                v_store_interleave(cov_data + j * 3, v_dst0, v_dst1, v_dst2);
            }
        }
#endif // CV_SIMD128

        for( ; j < size.width; j++ )
        {
            float dx = dxdata[j];
            float dy = dydata[j];

            cov_data[j*3] = dx*dx;	//Ix^2
            cov_data[j*3+1] = dx*dy;	//Ixy
            cov_data[j*3+2] = dy*dy;	//Iy^2
        }
    }

	//步骤3，均值滤波
    boxFilter(cov, cov, cov.depth(), Size(block_size, block_size),
        Point(-1,-1), false, borderType );

	//步骤4，计算特征值
    if( op_type == MINEIGENVAL )
        calcMinEigenVal( cov, eigenv );
    else if( op_type == HARRIS )
        calcHarris( cov, eigenv, k );
    else if( op_type == EIGENVALSVECS )
        calcEigenValsVecs( cov, eigenv );
}

calcHarris函数：

static void calcHarris( const Mat& _cov, Mat& _dst, double k )
{
    int i, j;
    Size size = _cov.size();
#if CV_TRY_AVX
    bool haveAvx = CV_CPU_HAS_SUPPORT_AVX;
#endif

    if( _cov.isContinuous() && _dst.isContinuous() )
    {
        size.width *= size.height;
        size.height = 1;
    }

	//遍历每个像素
    for( i = 0; i < size.height; i++ )
    {
        const float* cov = _cov.ptr<float>(i);
        float* dst = _dst.ptr<float>(i);

#if CV_TRY_AVX
        if( haveAvx )
            j = calcHarrisLine_AVX(cov, dst, k, size.width);
        else
#endif // CV_TRY_AVX
            j = 0;

//指令集加速
#if CV_SIMD128
        {
            v_float32x4 v_k = v_setall_f32((float)k);

            for( ; j <= size.width - v_float32x4::nlanes; j += v_float32x4::nlanes )
            {
                v_float32x4 v_a, v_b, v_c;
                v_load_deinterleave(cov + j * 3, v_a, v_b, v_c);

                v_float32x4 v_ac_bb = v_a * v_c - v_b * v_b;
                v_float32x4 v_ac = v_a + v_c;
                v_float32x4 v_dst = v_ac_bb - v_k * v_ac * v_ac;
                v_store(dst + j, v_dst);
            }
        }
#endif // CV_SIMD128

        for( ; j < size.width; j++ )
        {
            float a = cov[j*3];		//自相关元素a
            float b = cov[j*3+1];	//自相关元素b
            float c = cov[j*3+2];	//自相关元素c
            dst[j] = (float)(a*c - b*b - k*(a + c)*(a + c));	//计算特征值R=ac-b*b-k*(a+c)*(a+c)
        }
    }
}

实战

#include <iostream>
#include "opencv2/core/core.hpp"
#include "opencv2/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"

using namespace cv;
using namespace std;

int main()
{
	const char* imageName = "SiliconChip.bmp";
	cv::Mat src = imread(imageName, IMREAD_GRAYSCALE);

	int blocksize = 3;
	int apertureSize = 7;
	double k = 0.001;
	int thresh = 80; //步骤5中阈值

	Mat dst;
	cv::cornerHarris(src, dst, blocksize, apertureSize, k, BORDER_DEFAULT);

	//归一化R图像，使其在0-255之间
	Mat dst_norm;
	cv::normalize(dst, dst_norm, 0, 255, NORM_MINMAX, CV_32FC1);
	//将R图像转换成CV_8U类型，以便能够正确显示
	Mat dst_norm_scaled;
	cv::convertScaleAbs(dst_norm, dst_norm_scaled);

	Mat srcColor;
	cvtColor(src, srcColor, COLOR_GRAY2RGB);//转rgb
	//步骤5，在角点画圈
	for (int i = 0;i < dst.rows;i++)
	{
		for (int j = 0;j < dst.cols;j++)
		{
			if (dst_norm_scaled.at<uchar>(i, j) > thresh)
			{
				cv::circle(srcColor, Point(j, i), 15, Scalar(0, 0, 255), 1, 8, 0);
			}
		}
	}

	system("pause");
	return 0;
}

检测结果显示：