Harris角点检测及实现

角点：最直观的印象就是在水平、竖直两个方向上变化均较大的点，即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,IxyCvMat *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;}