Harris角点检测算法详解

不适用opencv的代码（转）

////////////////////////////////////////////////////////////////////// // Construction/Destruction

////////////////////////////////////////////////////////////////////// #define B(image,x,y) ((uchar *)(image->imageData+image->widthStep*(y)))[(x)*3] #define G(image,x,y) ((uchar *)(image->imageData+image->widthStep*(y)))[(x)*3+1] #define R(image,x,y) ((uchar *)(image->imageData+image->widthStep*(y)))[(x)*3+2] #define S(image,x,y) ((uchar *)(image->imageData+image->widthStep*(y)))[(x)]

//卷积计算求Ix,Iy,以及滤波

//a指向的数组是size1*size2（滤波器大小）大小的...求导，a 指滤波器，xwidth,ywidth指图像大小 px=i-size1/2+i1;py=j-size2/2+j1;分别消掉了size1/2，size2/2,所以这样就更好理解了。 for(i=size1/2;i

CvMat *mbys(CvMat *mat,intxwidth,intywidth,double *a,int size1,int size2)//size { inti,j; int i1,j1; intpx,py; int m; CvMat *mat1;

mat1=cvCloneMat(mat);

for(i=size1/2;i

//CV_MAT_ELEM访问矩阵元素

m+=CV_MAT_ELEM(*mat,double,px,py)*a[i1*size1+j1]; }

CV_MAT_ELEM(*mat1,double,i,j)=m; } return mat1; }

//计算Ix2,Iy2,Ixy

CvMat *mbxy(CvMat *mat1,CvMat *mat2,int xwidth,intywidth) { inti,j; CvMat *mat3;

mat3=cvCloneMat(mat1); for(i=0;i

CV_MAT_ELEM(*mat3,double,i,j)=CV_MAT_ELEM(*mat1,double,i,j)*CV_MAT_ELEM(*mat2,double,i,j);

} return mat3; }

//用来求得响应度

CvMat *mbcim(CvMat *mat1,CvMat *mat2,CvMat *mat3,int xwidth,intywidth) { inti,j; CvMat *mat;

mat=cvCloneMat(mat1); for(i = 0; i

for(j = 0; j

//注意：要在分母中加入一个极小量以防止除数为零溢出

CV_MAT_ELEM(*mat,double,i,j)=(CV_MAT_ELEM(*mat1,double,i,j)*CV_MAT_ELEM(*mat2,double,i,j)- CV_MAT_ELEM(*mat3,double,i,j)*CV_MAT_ELEM(*mat3,double,i,j))/

(CV_MAT_ELEM(*mat1,double,i,j)+CV_MAT_ELEM(*mat2,double,i,j)+0.000001); } } return mat; }

//用来求得局部极大值

CvMat *mblocmax(CvMat *mat1,int xwidth,intywidth,int size) { inti,j;

double max=-1000; int i1,j1; intpx,py; CvMat *mat;

mat=cvCloneMat(mat1);

for(i=size/2;i

for(i1=0;i1

for(j1=0;j1

if(CV_MAT_ELEM(*mat1,double,px,py)>max) max=CV_MAT_ELEM(*mat1,double,px,py);

} if(max>0)

CV_MAT_ELEM(*mat,double,i,j)=max; else

CV_MAT_ELEM(*mat,double,i,j)=0; } return mat; }

//用来确认角点

CvMat *mbcorner(CvMat *mat1,CvMat *mat2,int xwidth,intywidth,intsize,double thresh) {

CvMat *mat; inti,j;

mat=cvCreateMat(ywidth,xwidth,CV_32FC1); for(i=size/2;i

if(CV_MAT_ELEM(*mat1,double,i,j)==CV_MAT_ELEM(*mat2,double,i,j))//首先取得局部极大值 //mat1是CRF相应的所有。Mat2是CRF局部极大值点.

if(CV_MAT_ELEM(*mat1,double,i,j)>thresh)//然后大于这个阈值 CV_MAT_ELEM(*mat,int,i,j)=255;//满足上两个条件，才是角点！ else

CV_MAT_ELEM(*mat,int,i,j)=0; } return mat; } CvPoint* threshold) {

CvMat *mat_I,*mat_Ix,*mat_Iy,*mat_Ixy,*mat_Ix2,*mat_Iy2;//相应的矩阵 IplImage *pImgGray=NULL; //灰度图像 IplImage *dst=NULL; //目标图像

IplImage *pImgDx=NULL; //水平梯度卷积后的图像 IplImage *pImgDy=NULL; //竖起梯度卷积后的图像 IplImage *pImgDx2=NULL;//Ix2图像 IplImage *pImgDy2=NULL;//Iy2图像

CHarris::harris_features(IplImage

*src,intgausswidth,doublesigma,intsize,int

IplImage *pImgDxy=NULL;//Ixy图像

pImgGray=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1); dst=cvCreateImage(cvGetSize(src),src->depth,3);

pImgDx=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);//创建图像 pImgDy=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1); pImgDx2=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1); pImgDy2=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1); pImgDxy=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1); constintcxDIB=src->width ; // 图像宽度 constintcyDIB=src->height; // 图像高度 double *I=new double[cxDIB*cyDIB];

cvCvtColor(src,pImgGray,CV_RGB2GRAY);//灰度化 dst=cvCloneImage(src); inti,j;

for(j=0;j

I[j*cxDIB+i]=S(pImgGray,i,j);//将灰度图像数值存入I中 }

mat_I=cvCreateMat(cyDIB,cxDIB,CV_64FC1);

cvInitMatHeader(mat_I,cyDIB,cxDIB,CV_64FC1,I);//用I来初始化相应的矩阵 // cout<

//-------------------------------------------------------------------------- // 第一步：利用差分算子对图像进行滤波

//-------------------------------------------------------------------------- //定义水平方向差分算子并求Ix double dx[9]={-1,0,1,-1,0,1,-1,0,1};

mat_Ix=mbys(mat_I,cxDIB,cyDIB,dx,3,3); //求Ix矩阵 // cout<

//定义垂直方向差分算子并求Iy double dy[9]={-1,-1,-1,0,0,0,1,1,1};

mat_Iy=mbys(mat_I,cxDIB,cyDIB,dy,3,3);//求Iy矩阵 // cout<

S(pImgDx,i,j)=CV_MAT_ELEM(*mat_Ix,double,j,i);//为相应图像赋值 S(pImgDy,i,j)=CV_MAT_ELEM(*mat_Iy,double,j,i); }

mat_Ix2=mbxy(mat_Ix,mat_Ix,cxDIB,cyDIB);//计算Ix2,Iy2,Ixy矩阵 mat_Iy2=mbxy(mat_Iy,mat_Iy,cxDIB,cyDIB);

mat_Ixy=mbxy(mat_Ix,mat_Iy,cxDIB,cyDIB); for(j=0;j

S(pImgDxy,i,j)=CV_MAT_ELEM(*mat_Ixy,double,j,i);//为相应图像赋值 S(pImgDx2,i,j)=CV_MAT_ELEM(*mat_Ix2,double,j,i); S(pImgDy2,i,j)=CV_MAT_ELEM(*mat_Iy2,double,j,i); }

//-------------------------------------------------------------------------- // 第二步：对Ix2/Iy2/Ixy进行高斯平滑，以去除噪声

//-------------------------------------------------------------------------- //本例中使用5×5的高斯模板 //计算模板参数 //intgausswidth=5; //double sigma=0.8;

double *g=new double[gausswidth*gausswidth]; for(i=0;i

g[i*gausswidth+j]=exp(-((i-int(gausswidth/2))*(i-int(gausswidth/2))+(j-int(gausswidth/2))*(j-int(gausswidth/2)))/(2*sigma));

//归一化：使模板参数之和为1（其实此步可以省略） double total=0;

for(i=0;i

for(i=0;i

//进行高斯平滑

mat_Ix2=mbys(mat_Ix2,cxDIB,cyDIB,g,gausswidth,gausswidth); mat_Iy2=mbys(mat_Iy2,cxDIB,cyDIB,g,gausswidth,gausswidth); mat_Ixy=mbys(mat_Ixy,cxDIB,cyDIB,g,gausswidth,gausswidth);

//-------------------------------------------------------------------------- // 第三步：计算角点量

//--------------------------------------------------------------------------

//计算cim：即cornerness of image，我们把它称做‘角点量’ CvMat *mat_cim;

mat_cim=mbcim(mat_Ix2,mat_Iy2,mat_Ixy,cxDIB,cyDIB);

// cout<

//-------------------------------------------------------------------------- // 第四步：进行局部非极大值抑制

//-------------------------------------------------------------------------- CvMat *mat_locmax; //constint size=7;

mat_locmax=mblocmax(mat_cim,cxDIB,cyDIB,size);

// cout<

//-------------------------------------------------------------------------- // 第五步：获得最终角点

//-------------------------------------------------------------------------- CvMat *mat_corner;

//const double threshold=4500; //intcornernum=0;

mat_corner=mbcorner(mat_cim,mat_locmax,cxDIB,cyDIB,gausswidth,threshold); //CCommonCommonClass; CvPointpt[5000];

for(j=size/2;j

if(CV_MAT_ELEM(*mat_corner,int,j,i)==255) { pt[cornerno].x=i; pt[cornerno].y=j; cornerno++;

// CommonClass.DrawCross(showImg2,pt,CV_RGB(0,0,255),1,4); // cvCircle(dst,pt,2,CV_RGB(255,0,0),1,8,0); // cout<

//-------------------------------------------------------------------------- // 第四步：进行局部非极大值抑制

//-------------------------------------------------------------------------- CvMat *mat_locmax; //constint size=7;

mat_locmax=mblocmax(mat_cim,cxDIB,cyDIB,size);

// cout<

//-------------------------------------------------------------------------- // 第五步：获得最终角点

//-------------------------------------------------------------------------- CvMat *mat_corner;

//const double threshold=4500; //intcornernum=0;

mat_corner=mbcorner(mat_cim,mat_locmax,cxDIB,cyDIB,gausswidth,threshold); //CCommonCommonClass; CvPointpt[5000];

for(j=size/2;j

if(CV_MAT_ELEM(*mat_corner,int,j,i)==255) { pt[cornerno].x=i; pt[cornerno].y=j; cornerno++;

// CommonClass.DrawCross(showImg2,pt,CV_RGB(0,0,255),1,4); // cvCircle(dst,pt,2,CV_RGB(255,0,0),1,8,0); // cout<

本文来源：https://www.bwwdw.com/article/poi7.html