harris算子也称为Plessey算子,是为了改善Moravec算子性能,由Harris和Stephens提出的。经典的Harris角点检测是基于亮度变化的角点检测算法,具有较高的稳定性和稳健性,能够在图像旋转,灰度变化以及噪声干扰等情况下准确的检测到角点。
它的计算步骤如下:
1.首先计算矩阵M
其中Ix是图像该点对x方向的偏导数,Iy是对y方向的偏导数,w是你所选择的邻域,可以是3*3或者是5*5等等都可以。这个矩阵是一个实对称矩阵,它有两个特征值,如果这两个特征值都比较大,那么该点就是角点;如果一个大一个小,那么该点就是边缘点;如果都比较小,就说明该点的灰度变化比较平缓,是一个平面。
2.在角点判别中,一般不需要计算两个特征值,而是采用角点响应函数R作为检测角点特征的依据:
根据响应值R的大小可以区分出图像中的边缘和角点,所以可以在这里设置一个阈值来判断该点是否为角点,阈值提高,可以减少检测到的角点数目。邻域大小的改变会影响提取角点的数目和容忍度,通常在3*3或5*5的邻域内进行非极大抑制,局部极大值即为图像中的角点。
我用opencv2.3.1+vs2008实现的代码如下:
#include "stdafx.h"
#include <opencv2/opencv.hpp>
#include "highgui.h"
#include <math.h>
typedef unsigned long uint32;
typedef unsigned int uint16;
typedef unsigned char uint8;
#define THRESHOLD 5
IplImage *src_gray1, *src_gray2, *src_gray3;
IplImage* src_img, *dst_img;
void AllocateImage(IplImage* I) //给图像分配大小
{
CvSize sz = cvGetSize(I);
dst_img = cvCreateImage( sz, IPL_DEPTH_8U, 1);
cvSetZero(dst_img);
src_gray1 = cvCreateImage( sz, IPL_DEPTH_8U, 1); //原图的三个通道
src_gray2 = cvCreateImage( sz, IPL_DEPTH_8U, 1);
src_gray3 = cvCreateImage( sz, IPL_DEPTH_8U, 1);
}
void DeallocateImage()
{
cvReleaseImage(&src_img);
cvReleaseImage(&src_gray1);
cvReleaseImage(&src_gray2);
cvReleaseImage(&src_gray3);
}
void Harris_check(IplImage* I, IplImage* dst, int n) //harris角点检测函数
{
int i,j,x,y,num;
int Value,ValueX0,ValueY0;
int IX,IY;
long M00,M01,M11;
double DetM,TraceM;
double *R = new double[I->width*I->height];
for( i=1+n/2; i<(I->height-1-n/2); i++ ) //计算角点响应函数R
{
for( j=1+n/2; j<(I->width-1-n/2); j++ )
{
M00=0;
M01=0;
M11=0;
for( y=0; y<n; y++ )
{
for( x=0; x<n; x++ )
{
Value = *( I->imageData + (i-n/2+y)*I->widthStep + (j-n/2+x) );
ValueX0 = *( I->imageData + (i-n/2+y)*I->widthStep + (j-n/2+x-1) );
ValueY0 = *( I->imageData + (i-n/2+y-1)*I->widthStep + (j-n/2+x) );
IX = Value - ValueX0;
IY = Value - ValueY0;
M00 += IX*IX;
M01 += IX*IY;
M11 += IY*IY;
}
}
DetM = M00*M11 - M01*M01;
TraceM = M00 + M11;
R[i*I->width+j] = DetM - 0.05*TraceM*TraceM;
}
}
for( i=1+n/2; i<(I->height-1-n/2); i++ ) //非极大抑制
{
for( j=1+n/2; j<(I->width-1-n/2); j++ )
{
if( R[i*I->width+j]>THRESHOLD)
{
num=0;
for(y=0; y<n; y++)
{
for(x=0; x<n; x++)
{
if(R[i*I->width+j]>R[(i-n/2+y)*I->width+j-n/2+x])
num++;
}
}
if(num==n*n-1)
{
*(dst->imageData+i*dst->widthStep+j) =255;
}
}
}
}
}
int _tmain(int argc, _TCHAR* argv[])
{
src_img = cvLoadImage("角点.bmp");
AllocateImage(src_img);
cvSplit( src_img, src_gray1, src_gray2, src_gray3, 0);
Harris_check(src_gray1, dst_img, 5);
cvNamedWindow("my picture",CV_WINDOW_AUTOSIZE);
cvNamedWindow("my dst",CV_WINDOW_AUTOSIZE);
cvShowImage("my picture",src_img);
cvShowImage("my dst",dst_img);
cvWaitKey(0);
DeallocateImage();
cvDestroyWindow("my picture");
cvDestroyWindow("my dst");
return 0;
}