Harris角點檢測

http://www.cnblogs.com/doucontorl/archive/2011/01/02/1924157.html

角點：最直觀的印象就是在水平、豎直兩個方向上變化均較大的點，即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;}