本文共 16111 字，大约阅读时间需要 53 分钟。

参考博客：

今天想把Harris角点检测再捋一下，从公式到代码

 

在归一化特征图像上显示角点

在原图上显示角点 

opencv函数的官方文档

 具体代码：

#include 
   
    #include 
    
     using namespace cv;using namespace std;void call_back(int, void*);Mat Img_scr, Img_dst, Img_gray;int threshod_val = 30; //当前阈值 int max_threshod_val = 155;//最大阈值  int main(){	Img_scr = imread("img1.jpg");	imshow("original", Img_scr);	cvtColor(Img_scr, Img_gray, CV_BGR2GRAY);	namedWindow("work");	createTrackbar("thresh", "work", &threshod_val, max_threshod_val, call_back);	call_back(threshod_val, 0);	while (1)	{		waitKey(0);	}	return 0;}void call_back(int, void*){	Mat normImage, scaledImage;	//置零当前需要显示的两幅图，即清除上一次调用此函数时他们的值  	Mat Img_scr1 = Img_scr.clone();	Img_dst = Mat::zeros(Img_scr.size(), CV_32FC1);	//进行角点检测  	cornerHarris(Img_gray, Img_dst, 2, 3, 0.04, BORDER_DEFAULT);	// 归一化与转换  	normalize(Img_dst, normImage, 0, 255, NORM_MINMAX, CV_32FC1, Mat());	convertScaleAbs(normImage, scaledImage);//将归一化后的图线性变换成8位无符号整型   	// 将检测到的，且符合阈值条件的角点绘制出来 	for (int i = 0; i < normImage.rows; i++)	{		for (int j = 0; j< normImage.cols; j++)		{			if ((int)normImage.at
     
      (i, j) > threshod_val + 100)			{  //画圆				circle(Img_scr1, Point(i, j), 5, Scalar(10, 10, 255), 2, 8, 0);				circle(scaledImage, Point(i, j), 5, Scalar(0, 10, 255), 2, 8, 0);			}		}	}	imshow("work", Img_scr1);	imshow("work_1", scaledImage);}
     
    
   

源代码理解：

static voidcornerEigenValsVecs( const Mat& src, Mat& eigenv, int block_size,                     int aperture_size, int op_type, double k=0.,                     int borderType=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 );}

加上注释：

opencv 在Harris源代码里面并没有加上高斯滤波，是因为它希望这是一个单独的模块

我们进行工程的时候，图片先进行高斯滤波，然后对图片使用Harris角点检测

边缘检测的算法有很多种，这些算法通常容易受到图像本身的一些噪声的干扰，尤其当用偏微分方程获取图像边缘时候，如果边缘不连续，甚至导致函数水平集无法停止收敛。

高斯滤波通过平滑，可以去除噪声等干扰，是可以提高边缘检测的精准度。

Harris虽然是角点检测，但是内部却使用了索贝尔(Sobel)算子,计算差分，进行近似求导

static voidcornerEigenValsVecs( const Mat& src, Mat& eigenv, int block_size,                     int aperture_size, int op_type, double k=0.,                     int borderType=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)算子,实现任意阶导数和混合偏导数		//差分，即近似一阶导数        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滤波器运算符计算x或y方向的图像差分。		//其实它的参数变量和Sobel基本上是一样的，除了没有ksize核的大小。        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()    //当normalize = false时，为非归一化的盒式滤波，用于计算每个像素邻域内的积分特性，    //比如密集光流算法中用到的图像倒数的协方差矩阵    //方框滤波W(x,y)卷积, 也可用高斯核加权...    //W(Y,Y)与矩阵cov卷积运算得到 H 矩阵,后面通过H矩阵的特征值决定是否是角点    //这里默认采用了权值都为1的模板    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 );}

根据公式 R = det(H) - k*trace(H)*trace(H);

根据公式 R = a1*a2 - k*(a1+a2)*(a1+a2)

矩阵行列式  等于   其特征值乘积

矩阵的迹  等于   其矩阵特征值的和 

 补充：

矩阵行列式  等于   其特征值乘积

特征行列式：

|λI-A|=(λ-k1)(λ-k2)...(λ-kn)

其中k1,k2,...,kn是n个特征值

令上式中的λ=0，得到

|-A|=(0-k1)(0-k2)...(0-kn)

即(-1)^n|A|=(-1)^nk1k2...kn

则|A|=k1k2...kn

矩阵的迹  等于   其矩阵特征值的和 

static voidcalcHarris( const Mat& _cov, Mat& _dst, double k ){    int i, j;    Size size = _cov.size();#if CV_SSE    volatile bool simd = checkHardwareSupport(CV_CPU_SSE);#endif    if( _cov.isContinuous() && _dst.isContinuous() )    {        size.width *= size.height;        size.height = 1;    }    for( i = 0; i < size.height; i++ )    {        const float* cov = (const float*)(_cov.data + _cov.step*i);        float* dst = (float*)(_dst.data + _dst.step*i);        j = 0;    #if CV_SSE        if( simd )        {            __m128 k4 = _mm_set1_ps((float)k);            for( ; j <= size.width - 5; j += 4 )            {                __m128 t0 = _mm_loadu_ps(cov + j*3); // a0 b0 c0 x                __m128 t1 = _mm_loadu_ps(cov + j*3 + 3); // a1 b1 c1 x                __m128 t2 = _mm_loadu_ps(cov + j*3 + 6); // a2 b2 c2 x                __m128 t3 = _mm_loadu_ps(cov + j*3 + 9); // a3 b3 c3 x                __m128 a, b, c, t;                t = _mm_unpacklo_ps(t0, t1); // a0 a1 b0 b1                c = _mm_unpackhi_ps(t0, t1); // c0 c1 x x                b = _mm_unpacklo_ps(t2, t3); // a2 a3 b2 b3                c = _mm_movelh_ps(c, _mm_unpackhi_ps(t2, t3)); // c0 c1 c2 c3                a = _mm_movelh_ps(t, b);                b = _mm_movehl_ps(b, t);                t = _mm_add_ps(a, c);                a = _mm_sub_ps(_mm_mul_ps(a, c), _mm_mul_ps(b, b));                t = _mm_mul_ps(_mm_mul_ps(k4, t), t);                a = _mm_sub_ps(a, t);                _mm_storeu_ps(dst + j, a);            }        }    #endif        //根据公式 R = det(H) - k*trace(H)*trace(H);        //根据公式 R = a1*a2 - k*(a1+a2)*(a1+a2)        //矩阵行列式  等于   其特征值乘积        //矩阵的迹  等于   其矩阵特征值的和        for( ; j < size.width; j++ )        {            float a = cov[j*3];            float b = cov[j*3+1];            float c = cov[j*3+2];            dst[j] = (float)(a*c - b*b - k*(a + c)*(a + c));        }    }}

最后整理下上面的源代码，简化后

//--------------------【cornerHarris源码分析】------------------------------------//Harries角点实现函数,截取cornerHarries中的关键代码做了简化void myHarries(const Mat& src, Mat& eigenv, int block_size, int aperture_size, double k = 0.1){	eigenv.create(src.size(), CV_32F);	Mat eigenv_a(src.size(), CV_32FC1, Scalar(0));	eigenv_a.copyTo(eigenv);	Mat Dx, Dy;	//soble operation get Ix, Iy	Sobel(src, Dx, CV_32F, 1, 0, aperture_size);	Sobel(src, Dy, CV_32F, 0, 1, aperture_size);	//get covariance matrix	Size size = src.size();	Mat cov(size, CV_32FC3);     //创建一个三通道cov矩阵分别存储[Ix*Ix, Ix*Iy, Iy*Iy];	for (int i = 0; i < size.height; i++)	{		float* cov_data = cov.ptr
   
    (i);		const float* dxdata = Dx.ptr
    
     (i);		const float* dydata = Dy.ptr
     
      (i);		for (int j = 0; j < size.width; j++)		{			float dx = dxdata[j];			float dy = dydata[j];			cov_data[j * 3] = dx * dx;    //即  Ix*Ix			cov_data[j * 3 + 1] = dx*dy;  //即  Ix*Iy			cov_data[j * 3 + 2] = dy*dy;  //即  Iy*Iy		}	}	//方框滤波W(x,y)卷积, 也可用高斯核加权...	//W(Y,Y)与矩阵cov卷积运算得到 H 矩阵,后面通过H矩阵的特征值决定是否是角点	boxFilter(cov, cov, cov.depth(), Size(block_size, block_size), Point(-1, -1), false);	//cale Harris	size = cov.size();	//此处计算响应R= det(H) - k*trace(H)*trace(H);	for (int i = 0; i < size.height; i++)	{		const float* covPtr = cov.ptr
      
       (i);		float* dstPtr = eigenv.ptr
       
        (i);		int j = 0;		for (; j < size.width; j++)		{			float a = covPtr[j * 3];			float b = covPtr[j * 3 + 1];			float c = covPtr[j * 3 + 2];			//根据公式 R = det(H) - k*trace(H)*trace(H);			dstPtr[j] = (float)(a*c - b*b - k*(a + c)*(a + c));		}	}	double max, min;	minMaxLoc(eigenv, &min, &max);	//cout<< max << endl;	double threshold = 0.1*max;	cv::threshold(eigenv, eigenv, threshold, 1, CV_THRESH_BINARY);   //eigenv的类型是CV_32F,	imshow("eigenv", eigenv);}int main(){	Mat Img_scr, Img_gray;	Img_scr = imread("img1.jpg");	imshow("original", Img_scr);	cvtColor(Img_scr, Img_gray, CV_BGR2GRAY);	//进行Harris角点检测找出角点	Mat cornerStrength;	myHarries(Img_gray, cornerStrength, 2, 3, 0.01);	while (1)	{		waitKey(0);	}	return 0;}
       
      
     
    
   

最后的dst就是每一个点，Harris 角点检测准则的分值矩阵

对这些值进行非极大值抑制NMS

opencv 提供了函数，我们通过组合，可以自己写一个NMS

OpenCV中的dilate膨胀函数， 膨胀即是求局部最大值的操作，图像A与核B作卷积运算，计算核B覆盖区域的像素点的最大值

OpenCV中的compare比较函数，void compare(InputArray src1, InputArray src2, OutputArray dst, int cmpop);

enum { CMP_EQ=0,    //相等	CMP_GT=1,   //大于	CMP_GE=2,   //大于等于	CMP_LT=3,   //小于	CMP_LE=4,   //小于等于	CMP_NE=5 }; //不相等

OpenCV中的bitwise_and函数 是对二进制数据进行“与”操作，即对图像（灰度图像或彩色图像均可）每个像素值进行二进制“与”操作，1&1=1，1&0=0，0&1=0，0&0=0

---

刚刚得到的Harris 角点检测准则的分值矩阵，并不能可视化，为了可视化，把里面的值，转换到0-255，然后uint8

// 归一化与转换  

normalize(Img_dst, normImage, 0, 255, NORM_MINMAX, CV_32FC1, Mat());

//将归一化后的图线性变换成8位无符号整型 

convertScaleAbs(normImage, scaledImage);

我们前面得代码已经实现了非极大值抑制，你可以通过滑动条，更改最大值阈值，30-155区间

或者另外一种可视化的方式 直接显示哪些地方的 分值较高 白色部分为分值较高区域

#include 
   
    #include 
    
     using namespace cv;using namespace std;int main(){	Mat Img_scr, Img_gray;	Img_scr = imread("img1.jpg");	imshow("original", Img_scr);	cvtColor(Img_scr, Img_gray, CV_BGR2GRAY);	//进行Harris角点检测找出角点	Mat cornerStrength;	cornerHarris(Img_gray, cornerStrength, 2, 3, 0.01);	//对灰度图进行阈值操作，得到二值图并显示	Mat harrisCorner;	threshold(cornerStrength, harrisCorner, 0.00001, 255, THRESH_BINARY);	imshow("110", harrisCorner);	while (1)	{		waitKey(0);	}	return 0;}
    
   

---

突发奇想，想画出Harris的椭圆

参考文档：

 

 

 

1.矩阵特征值和椭圆长短轴的关系？

椭圆的长短轴分别沿着矩阵A的两个特征向量的方向，而两个与之对应的特征值分别是半长轴和半短轴的长度的平方的倒数。

2.harris的H矩阵和协方差矩阵的关系？

协方差矩阵：

harris的H矩阵 

协方差矩阵的作用为什么比方差和均值要大呢？显而易见方差和均值只是一维随机变量的统计值，而协方差就不一样了，它可以表示多维随机变量之间的相关性信息。协方差矩阵的一个很出色的应用就是在PCA中，选择主方向。协方差矩阵的对角线的元素表示的是各个维度的方差，而非对角线上的元素表示的是各个维度之间的相关性，因此，在PCA中，我们尽量将非对角线上的元素化为0，即将矩阵对角化，选特征值较大的维度，去掉特征值较小的维度，来获得主方向，并且使主方向与其他方向的相关性尽量小。

那现在看看这个矩阵H,通过上面对协方差的描述，我们完全可以把这个矩阵看做一个二维随机分布的协方差矩阵，那么我们要做的就是将其对角化，求矩阵的两个特征值，然后根据这两个特征值来判断是不是角点（两个特征值都大代表角点）。 

3.H矩阵的关键性

虽然我们利用E(u,v)来描述角点的基本思想，然而最终我们仅仅使用的是矩阵H。

让我们看看矩阵H形式，是不是跟协方差矩阵形式很像，像归像，但是还是有些不同，哪儿不同？

一般协方差矩阵对应维的随机变量需要减去该维随机变量的均值，但矩阵H中并没有这样做，所以在矩阵H里，我们先进行各维的均值化处理，那么各维所对应的随机变量的均值为0，协方差矩阵就大大简化了，简化的最终结果就是矩阵H，是否明白了？我们的目的是分析数据的主要成分，相信了解PCA原理的，应该都了解均值化的作用。

如果我们对协方差矩阵H进行对角化，很明显，特征值就是主分量上的方差，这点大家应该明白吧？不明白的话可以复习下PCA原理。如果存在两个主分量所对应的特征值都比较大，说明什么？

      像素点的梯度分布比较散，梯度变化程度比较大，符合角点在窗口区域的特点；

      如果是平坦区域，那么像素点的梯度所构成的点集比较集中在原点附近，因为窗口区域内的像素点的梯度幅值非常小，此时矩阵M的对角化的两个特征值比较小；

      如果是边缘区域，在计算像素点的x、y方向上的梯度时，边缘上的像素点的某个方向的梯度幅值变化比较明显，另一个方向上的梯度幅值变化较弱，其余部分的点都还是集中原点附近，这样M对角化后的两个特征值理论应该是一个比较大，一个比较小，当然对于边缘这种情况，可能是呈45°的边缘，致使计算出的特征值并不是都特别的大，总之跟含有角点的窗口的分布情况还是不同的。

4.opencv画出harris椭圆

AffineMatch功能使用HarrisLaplace检测器和仿射适应来提取仿射不变特征。同时使用SIFT描述符描述局部特征并使用flann和vfc进行匹配。

代码有一点问题，就是椭圆的长轴短轴，不知道哪里出问题了，太小了，无法画椭圆，就*30000扩大了，然后画出来了，意思达到了，我有点强迫症！

#include 
   
    #include 
    
     using namespace cv;using namespace std;void myHarries(const Mat& src, Mat& eigenv, int block_size, int aperture_size, double k = 0.1){	eigenv.create(src.size(), CV_32F);	Mat eigenv_a(src.size(), CV_32FC1, Scalar(0));	eigenv_a.copyTo(eigenv);	Mat Dx, Dy;	//soble operation get Ix, Iy	Sobel(src, Dx, CV_32F, 1, 0, aperture_size);	Sobel(src, Dy, CV_32F, 0, 1, aperture_size);	//get covariance matrix	Size size = src.size();	Mat cov(size, CV_32FC3);     //创建一个三通道cov矩阵分别存储[Ix*Ix, Ix*Iy, Iy*Iy];	for (int i = 0; i < size.height; i++)	{		float* cov_data = cov.ptr
     
      (i);		const float* dxdata = Dx.ptr
      
       (i);		const float* dydata = Dy.ptr
       
        (i);		for (int j = 0; j < size.width; j++)		{			float dx = dxdata[j];			float dy = dydata[j];			cov_data[j * 3] = dx * dx;    //即  Ix*Ix			cov_data[j * 3 + 1] = dx*dy;  //即  Ix*Iy			cov_data[j * 3 + 2] = dy*dy;  //即  Iy*Iy		}	}	//方框滤波W(x,y)卷积, 也可用高斯核加权...	//W(Y,Y)与矩阵cov卷积运算得到 H 矩阵,后面通过H矩阵的特征值决定是否是角点	boxFilter(cov, cov, cov.depth(), Size(block_size, block_size), Point(-1, -1), false);	//cale Harris	size = cov.size();	Mat harris_H(Size(2, 2), CV_32FC1, Scalar(0));	float ax1 = 0, ax2 = 0;	float phi = 0;	vector
        
          harris_H_ellipse; KeyPoint temp; //此处计算响应R= det(H) - k*trace(H)*trace(H); for (int i = 0; i < size.height; i++) { const float* covPtr = cov.ptr
         
          (i); float* dstPtr = eigenv.ptr
          
           (i); int j = 0; for (; j < size.width; j++) { float a = covPtr[j * 3]; float b = covPtr[j * 3 + 1]; float c = covPtr[j * 3 + 2]; //根据公式 R = det(H) - k*trace(H)*trace(H); //float det_H = a*c - b*b; // v1*v2 //float trace = a + c;// v1+v2 //det_H = trace*v2 - v2 * v2; // v2 * v2 - trace*v2 + det_H = 0 //A B C = 1 -trace det_H //先求 B*B - 4*A*C //float delta = (-trace)*(-trace) - 4 * det_H; //float x1 = (trace + sqrt(delta)) / (2 * a); //float x2 = (trace - sqrt(delta)) / (2 * a); harris_H.at
           
            (0, 0) = a; harris_H.at
            
             (0, 1) = b; harris_H.at
             
              (1, 0) = b; harris_H.at
              
               (1, 1) = c; Mat uVal, uV; eigen(harris_H, uVal, uV); //double sqrt( double x ) ax1 = 1.f / sqrt(std::abs(uVal.at
               
                (0, 0)));// ax2 = 1.f / sqrt(std::abs(uVal.at
                
                 (1, 0))); phi = float(atan(uV.at
                 
                  (1, 0) / uV.at
                  
                   (0, 0)) * (180) / CV_PI); temp.pt.x = ax1; temp.pt.y = ax2; temp.angle = phi; harris_H_ellipse.push_back(temp); dstPtr[j] = (float)(a*c - b*b - k*(a + c)*(a + c)); } } Mat normImage, scaledImage; // 归一化与转换 normalize(eigenv, normImage, 0, 255, NORM_MINMAX, CV_32FC1, Mat()); convertScaleAbs(normImage, scaledImage);//将归一化后的图线性变换成8位无符号整型 int threshod_val = 35; // 将检测到的，且符合阈值条件的角点绘制出来 Point center; // 中心点坐标 Size axes; // 椭圆长短轴 Mat ynh_d = src.clone(); cv::cvtColor(ynh_d, ynh_d, CV_GRAY2BGR); for (int i = 0; i < normImage.rows; i++) { for (int j = 0; j< normImage.cols; j++) { if ((int)normImage.at
                   
                    (i, j) > 20 + 100) { //画圆 center.x = j; center.y = i; if (harris_H_ellipse[i*normImage.cols + j].pt.x >harris_H_ellipse[i*normImage.cols + j].pt.y) { axes.width = (int)harris_H_ellipse[i*normImage.cols + j].pt.x * 30000; axes.height = harris_H_ellipse[i*normImage.cols + j].pt.y * 30000; } else{ axes.height = harris_H_ellipse[i*normImage.cols + j].pt.x * 30000; axes.width = harris_H_ellipse[i*normImage.cols + j].pt.y * 30000; } double angle = harris_H_ellipse[i].angle; // 角度 // 绘制椭圆图像 ellipse(ynh_d, center, axes, angle * 180 / CV_PI, 0, 360, Scalar(0, 0, 255), 2, 8); // 绘制中心点坐标 circle(ynh_d, center, 1, Scalar(255, 0, 0)); } } } cv::imshow("scaledImage", scaledImage); cv::imshow("ynh_d", ynh_d);}int main(){ Mat Img_scr, Img_gray; Img_scr = imread("img1.jpg"); cv::imshow("original", Img_scr); cvtColor(Img_scr, Img_gray, CV_BGR2GRAY); //进行Harris角点检测找出角点 Mat cornerStrength; myHarries(Img_gray, cornerStrength, 2, 3, 0.01); while (1) { waitKey(0); } return 0;}
                   
                  
                 
                
               
              
             
            
           
          
         
        
       
      
     
    
   

 

转载地址：http://mcwlf.baihongyu.com/