角点：最直观的印象就是在水平、竖直两个方向上变化均较大的点，即Ix、Iy都较大
边缘：仅在水平、或者仅在竖直方向有较大的变化量，即Ix和Iy只有其一较大
平坦地区：在水平、竖直方向的变化量均较小，即Ix、Iy都较小
2 strong eigenvalues======interest point
1 strong eigenvalues======contour/edge
0 eigenvalues ======uniform region
角点响应
R=det(M)-k*(trace(M)^2) (k=0.04~0.06)
det(M)=λ1*λ2 trace(M)=λ1+λ2
R取决于M的特征值，对于角点|R|很大，平坦的区域|R|很小。
编程步骤：
使用opencv进行测试：
#include "stdafx.h" #include "cv.h" #include "highgui.h" void drawcross(CvArr* img,CvPoint2D32f pt) { const int radius=3; int ptx=cvRound(pt.x); int pty=cvRound(pt.y); int ls=ptx-radius; int re=ptx+radius; int us=pty-radius; int de=pty+radius; cvLine(img,cvPoint(ls,pty),cvPoint(re,pty),CV_RGB(0,0,255),1,0); cvLine(img,cvPoint(ptx,us),cvPoint(ptx,de),CV_RGB(0,0,255),1,0); } int main(int argc, char* argv[]) { CvPoint2D32f pt[100]; int cornercount=30; IplImage* srcimg=cvLoadImage("2.bmp"); IplImage* grayimg=cvCreateImage(cvGetSize(srcimg),IPL_DEPTH_8U,1); IplImage* eigimg=cvCreateImage(cvGetSize(srcimg),IPL_DEPTH_32F,1); IplImage* tempimg=cvCloneImage(eigimg); //cvConvertImage(srcimg,grayimg,0); cvCvtColor(srcimg,grayimg,CV_BGR2GRAY); cvGoodFeaturesToTrack(grayimg,eigimg,tempimg,pt,&cornercount,0.1,10,NULL,3,0,0.04); for(int i=0;i<cornercount;i++) { drawcross(srcimg,pt[i]); } cvNamedWindow("corner detection",CV_WINDOW_AUTOSIZE); cvShowImage("corner detection",srcimg); cvWaitKey(0); return 0; }
不适用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大小的...求导 CvMat *mbys(CvMat *mat,int xwidth,int ywidth,double *a,int size1,int size2) { int i,j; int i1,j1; int px,py; int m; CvMat *mat1; mat1=cvCloneMat(mat); for(i=size1/2;i<ywidth-size1/2;i++) for(j=size2/2;j<xwidth-size2/2;j++) { m=0; for(i1=0;i1<size1;i1++) for(j1=0;j1<size2;j1++) { px=i-size1/2+i1; py=j-size2/2+j1; //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,int ywidth) { int i,j; CvMat *mat3; mat3=cvCloneMat(mat1); for(i=0;i<ywidth;i++) for (j=0;j<xwidth;j++) { 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,int ywidth) { int i,j; CvMat *mat; mat=cvCloneMat(mat1); for(i = 0; i <ywidth; i++) { for(j = 0; j < xwidth; 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,int ywidth,int size) { int i,j; double max=-1000; int i1,j1; int px,py; CvMat *mat; mat=cvCloneMat(mat1); for(i=size/2;i<ywidth-size/2;i++) for(j=size/2;j<xwidth-size/2;j++) { max=-10000; for(i1=0;i1<size;i1++) for(j1=0;j1<size;j1++) { px=i-size/2+i1; py=j-size/2+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,int ywidth,int size,double thresh) { CvMat *mat; int i,j; mat=cvCreateMat(ywidth,xwidth,CV_32FC1); for(i=size/2;i<ywidth-size/2;i++) for(j=size/2;j<xwidth-size/2;j++) { if(CV_MAT_ELEM(*mat1,double,i,j)==CV_MAT_ELEM(*mat2,double,i,j))//首先取得局部极大值 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* CHarris::harris_features(IplImage *src,int gausswidth,double sigma,int size,int 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图像 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); const int cxDIB=src->width ; // 图像宽度 const int cyDIB=src->height; // 图像高度 double *I=new double[cxDIB*cyDIB]; cvCvtColor(src,pImgGray,CV_RGB2GRAY);//灰度化 dst=cvCloneImage(src); int i,j; for(j=0;j<cyDIB;j++) for(i=0;i<cxDIB;i++) { 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<<CV_MAT_ELEM(*mat_I,double,200,200)<<endl; //-------------------------------------------------------------------------- // 第一步：利用差分算子对图像进行滤波 //-------------------------------------------------------------------------- //定义水平方向差分算子并求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<<CV_MAT_ELEM(*mat_Ix,double,200,200)<<endl; //定义垂直方向差分算子并求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<<CV_MAT_ELEM(*mat_Iy,double,200,200)<<endl; for(j=0;j<cyDIB;j++) for(i=0;i<cxDIB;i++) { 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<cyDIB;j++) for(i=0;i<cxDIB;i++) { 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的高斯模板 //计算模板参数 //int gausswidth=5; //double sigma=0.8; double *g=new double[gausswidth*gausswidth]; for(i=0;i<gausswidth;i++)//定义模板 for(j=0;j<gausswidth;j++) 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<gausswidth*gausswidth;i++) total+=g[i]; for(i=0;i<gausswidth;i++) for(j=0;j<gausswidth;j++) g[i*gausswidth+j]/=total; //进行高斯平滑 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<<CV_MAT_ELEM(*mat_cim,double,cyDIB-1,cxDIB-1)<<endl; //-------------------------------------------------------------------------- // 第四步：进行局部非极大值抑制 //-------------------------------------------------------------------------- CvMat *mat_locmax; //const int size=7; mat_locmax=mblocmax(mat_cim,cxDIB,cyDIB,size); // cout<<CV_MAT_ELEM(*mat_locmax,double,cyDIB-1,cxDIB-1)<<endl; //-------------------------------------------------------------------------- // 第五步：获得最终角点 //-------------------------------------------------------------------------- CvMat *mat_corner; //const double threshold=4500; //int cornernum=0; mat_corner=mbcorner(mat_cim,mat_locmax,cxDIB,cyDIB,gausswidth,threshold); //CCommon CommonClass; CvPoint pt[5000]; for(j=size/2;j<cyDIB-size/2;j++) for(i=size/2;i<cxDIB-size/2;i++) { 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<<i<<" "<<j<<endl; } } return pt; }
