openCV2马拉松第19圈——Harris角点检測（自己实现）
计算机视觉讨论群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是一个二阶矩矩阵，能够由我们原始图片的差分得到
假设a或者b非常接近0，也就是说他们非常小，那么这就不是个角点。角点的地方a,b肯定都非常大
我们发现。E是一个关于U,V的二次曲线，令E(u,v) = CONST,那么E就是一个椭圆
而M,由矩阵对角化能够写成左边的格式，我们能够得到两个特征值
长轴短轴由特征值决定。椭圆方向由矩阵R决定，以下没实用到R,所以你能够忽略它
在这里，我们的特征值就派上用场了，当两个特征值都非常大非常相近。说明我的椭圆非常小，整个E变化非常剧烈。那么我就找到了我的角点。
以下是 Harris给出的公式
- det(M) =
- trace(M) =
那么非常多人就在纠结要怎么计算了。
要怎么计算了。det = Ix2*Iy2 - Ixy^2 trace = 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 .
对每一个像素 
计算 
协方差矩阵 
over
a 
neighborhood. 用的是harris响应函数
荷枪实弹
以下是官方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;
}另外附加福利
