首页 > 其他 > 详细

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

时间:2014-05-25 18:39:07      阅读:573      评论:0      收藏:0      [点我收藏+]
计算机视觉讨论群162501053
转载请注明:http://blog.csdn.net/abcd1992719g/article/details/26824529



收入囊中
  • 使用OpenCV的connerHarris实现角点检测
  • 自己实现Harris算法
下面是自己实现的一个效果图
bubuko.com,布布扣
因为阀值设置比较高,所以房屋周围没有找出来



葵花宝典
在此之前,我们讲过边缘的检测,边缘检测的基本原理就是x方向或者y方向梯度变化很大,角点,顾名思义,就是两个方向的梯度变化都很大。
bubuko.com,布布扣
左1,平滑区域,没有边缘和角点,窗口在任何方向移动都无变化
左2,边缘区域,在边缘方向移动没有变化
左3,角点区域,在任何方向移动都有显著的变化

下面我们定义

E(u,v)=(x,y)[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+ Iy

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

再代入下面的公式

E(u,v)=[I(x+u,y+v)?I(x,y)]2(x,y)

我们得到

bubuko.com,布布扣

因此,我们可以写成

bubuko.com,布布扣

其中,M是一个二阶矩矩阵,可以由我们原始图片的差分得到bubuko.com,布布扣

我们来想像一下,假使梯度是沿着x或者y方向的,也就是说Ix = 0或者是Iy = 0,那么我们的M就是
bubuko.com,布布扣
如果a或者b很接近0,也就是说他们很小,那么这就不是个角点,角点的地方a,b肯定都很大

我们发现,E是一个关于U,V的二次曲线,令E(u,v) = CONST,那么E就是一个椭圆
bubuko.com,布布扣而M,由矩阵对角化可以写成左边的格式,我们可以得到两个特征值
bubuko.com,布布扣长轴短轴由特征值决定,椭圆方向由矩阵R决定,下面没有用到R,所以你可以忽略它

bubuko.com,布布扣
在这里,我们的特征值就派上用场了,当两个特征值都很大很相近,说明我的椭圆很小,整个E变化很剧烈,那么我就找到了我的角点。
下面是 Harris给出的公式
bubuko.com,布布扣

  • det(M) = bubuko.com,布布扣
  • trace(M) = bubuko.com,布布扣

那么很多人就在纠结bubuko.com,布布扣要怎么计算了。
det = Ix2*Iy2 - Ixy^2
trace = Ix2 + Iy2
[还有另一种经验公式,bubuko.com,布布扣,我自己的实现就是采用了这一种]

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


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不具有,如下图

bubuko.com,布布扣


初识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  .

对每个像素 bubuko.com,布布扣计算 bubuko.com,布布扣 协方差矩阵 bubuko.com,布布扣 over a bubuko.com,布布扣 neighborhood. 用的是harris响应函数

bubuko.com,布布扣


荷枪实弹
下面是官方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 variables
Mat src, src_gray;
int thresh = 200;
int max_thresh = 255;

char* source_window = "Source image";
char* corners_window = "Corners detected";

/// Function header
void 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

openCV2马拉松第19圈——Harris角点检测(自己实现),布布扣,bubuko.com

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

原文:http://blog.csdn.net/abcd1992719g/article/details/26824529

(0)
(0)
   
举报
评论 一句话评论(0
关于我们 - 联系我们 - 留言反馈 - 联系我们:wmxa8@hotmail.com
© 2014 bubuko.com 版权所有
打开技术之扣,分享程序人生!