openCV2马拉松第19圈——Harris角点检测（自己实现）

来源：互联网发布：php用户管理系统代码编辑：程序博客网时间：2024/05/19 14:39

计算机视觉讨论群162501053

转载请注明：http://blog.csdn.net/abcd1992719g/article/details/26824529

收入囊中

使用OpenCV的connerHarris实现角点检测
自己实现Harris算法

下面是自己实现的一个效果图

因为阀值设置比较高，所以房屋周围没有找出来

葵花宝典

在此之前，我们讲过边缘的检测，边缘检测的基本原理就是x方向或者y方向梯度变化很大，角点，顾名思义，就是两个方向的梯度变化都很大。

左1，平滑区域，没有边缘和角点，窗口在任何方向移动都无变化

左2，边缘区域，在边缘方向移动没有变化

左3，角点区域，在任何方向移动都有显著的变化

下面我们定义

E(u,v)=∑(x,y)∈W [I(x+u,y+v)−I(x,y)]2（w是我们的窗口,[u,v]是我们的shift，也就是移动）

我们想了解微小移动对E到底有何影响

于是我们对I进行一阶泰勒展开

I(x+u,y+v)=I(x,y)+Ixu+Iyv+higherorder terms

≈I(x,y)+Ixu+Iyv

=I(x,y)+[Ix Iy ][u v]T

再代入下面的公式

E(u,v)=∑[I(x+u,y+v)−I(x,y)]2(x,y)∈W

我们得到

因此，我们可以写成

其中，M是一个二阶矩矩阵，可以由我们原始图片的差分得到

我们来想像一下，假使梯度是沿着x或者y方向的，也就是说Ix = 0或者是Iy = 0,那么我们的M就是

如果a或者b很接近0，也就是说他们很小，那么这就不是个角点，角点的地方a,b肯定都很大

我们发现，E是一个关于U,V的二次曲线，令E(u,v) = CONST,那么E就是一个椭圆

而M,由矩阵对角化可以写成左边的格式，我们可以得到两个特征值

长轴短轴由特征值决定，椭圆方向由矩阵R决定，下面没有用到R,所以你可以忽略它

在这里，我们的特征值就派上用场了，当两个特征值都很大很相近，说明我的椭圆很小，整个E变化很剧烈，那么我就找到了我的角点。

下面是 Harris给出的公式

det(M) = $\lambda_{1}\lambda_{2}$
trace(M) = $\lambda_{1}+\lambda_{2}$

那么很多人就在纠结 $\lambda_{1}\lambda_{2}$ 要怎么计算了。

det = Ix2*Iy2 - Ixy^2trace = Ix2 + Iy2

［还有另一种经验公式，

，我自己的实现就是采用了这一种］

说了这么多理论，实现起来其实没有那么困难，下面我们就来看看算法的步骤吧。

1. 用sobel算子分别计算出水平梯度和垂直梯度

-1 0 1

-2 0 2

-1 0 1

-1 -2 -1

0 0 0

1 2 1

2. 计算高斯二阶矩阵，就是一个窗口里所有差分和

3.计算响应函数R

4.设置阀值

5.非极大值抑制，判断这个点的响应函数R是不是周围最大的

强调一下，harris对仿射变换只有部分不变性质。平移和旋转具有covariant,但是scaling不具有，如下图

初识API

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

src – 8比特或者是浮点数矩阵.
dst – 存放harris响应函数，类型是CV_32FC1因为计算出的值比较大，size和src一样 .
blockSize – 窗口大小.
ksize – Sobel算子大小.
k – 一般是0.04-0.06.
borderType .

对每个像素 $(x, y)$ 计算 $2\times2$ 协方差矩阵 $M^{(x,y)}$ over a $\texttt{blockSize} \times \texttt{blockSize}$ neighborhood. 用的是harris响应函数

$\texttt{dst} (x,y) = \mathrm{det} M^{(x,y)} - k \cdot \left ( \mathrm{tr} M^{(x,y)} \right )^2$

荷枪实弹

下面是官方sample,我就不多解释了，纯粹是API的调用，要注意的就是得到响应函数后进行了归一化到［0，255］

#include "opencv2/highgui/highgui.hpp"#include "opencv2/imgproc/imgproc.hpp"#include <iostream>#include <stdio.h>#include <stdlib.h>using namespace cv;using namespace std;/// Global variablesMat src, src_gray;int thresh = 200;int max_thresh = 255;char* source_window = "Source image";char* corners_window = "Corners detected";/// Function headervoid cornerHarris_demo( int, void* );/** @function main */int main( int argc, char** argv ){  /// Load source image and convert it to gray  src = imread( argv[1], 1 );  cvtColor( src, src_gray, CV_BGR2GRAY );  /// Create a window and a trackbar  namedWindow( source_window, CV_WINDOW_AUTOSIZE );  createTrackbar( "Threshold: ", source_window, &thresh, max_thresh, cornerHarris_demo );  imshow( source_window, src );  cornerHarris_demo( 0, 0 );  waitKey(0);  return(0);}/** @function cornerHarris_demo */void cornerHarris_demo( int, void* ){  Mat dst, dst_norm, dst_norm_scaled;  dst = Mat::zeros( src.size(), CV_32FC1 );  /// Detector parameters  int blockSize = 2;  int apertureSize = 3;  double k = 0.04;  /// Detecting corners  cornerHarris( src_gray, dst, blockSize, apertureSize, k, BORDER_DEFAULT );  /// Normalizing  normalize( dst, dst_norm, 0, 255, NORM_MINMAX, CV_32FC1, Mat() );  convertScaleAbs( dst_norm, dst_norm_scaled );  /// Drawing a circle around corners  for( int j = 0; j < dst_norm.rows ; j++ )     { for( int i = 0; i < dst_norm.cols; i++ )          {            if( (int) dst_norm.at<float>(j,i) > thresh )              {               circle( dst_norm_scaled, Point( i, j ), 5,  Scalar(0), 2, 8, 0 );              }          }     }  /// Showing the result  namedWindow( corners_window, CV_WINDOW_AUTOSIZE );  imshow( corners_window, dst_norm_scaled );}

举一反三

下面是我自己的实现，我找了好多地方都没找到别人的源代码...

#include "opencv2/highgui/highgui.hpp"  #include "opencv2/imgproc/imgproc.hpp"  #include <cmath>#include <iostream>using namespace cv;   Mat harris(Mat &im, double sigma, int thresh, int radius){    Mat dx,dy,Ix,Iy,Ix2,Iy2,Ixy,cim;    Sobel( im, Ix, CV_64F, 1, 0, 3);   //算法第一步，计算水平垂直差分    Sobel( im, Iy, CV_64F, 0, 1, 3);    int ksize = max(1, (int)(6*sigma));    if(ksize % 2 == 0)    ksize++;GaussianBlur(Ix.mul(Ix), Ix2, Size(ksize, ksize), sigma);   //算法第二步，计算二阶高斯差分矩阵   GaussianBlur(Iy.mul(Iy), Iy2, Size(ksize, ksize), sigma);GaussianBlur(Ix.mul(Iy), Ixy, Size(ksize, ksize), sigma);//Harris corner measure//cim = (Ix2.*Iy2 - Ixy.^2)./(Ix2 + Iy2);cim = (Ix2.mul(Iy2) - Ixy.mul(Ixy)) / (Ix2+Iy2);            //算法第三步，计算响应函数，我使用了另外一种    Mat structedElement(radius, radius, CV_8U, Scalar(1));    Mat mx,norm_cim;    normalize( cim, norm_cim, 0, 255, NORM_MINMAX, CV_8U, Mat() );dilate(norm_cim, mx, structedElement);norm_cim = ( norm_cim == mx) & (norm_cim>thresh);           //算法第4第5步融合，非极大值抑制和阀值检测return norm_cim;}int main( int, char** argv )  {  Mat src,gray;    src = imread( argv[1] );    cvtColor( src, gray, CV_RGB2GRAY );    Mat corners = harris(gray, 1.5, 30, 2);for( int j = 0; j < corners.rows ; j++ ) { for( int i = 0; i < corners.cols; i++ ) {            if( corners.at<unsigned char>(j,i) > 0){circle( gray, Point( i, j ), 3,  Scalar(0), 2, 8, 0 );}}}        namedWindow("result", 1);     imshow("result", gray);    waitKey();      return 0;  }

另外附加福利

harris matlab版本 http://www.cs.illinois.edu/~slazebni/spring13/harris.m

harris python版本 http://www.janeriksolem.net/2009/01/harris-corner-detector-in-python.html

11 1