角点:最直观的印象就是在水平、竖直两个方向上变化均较大的点,即Ix、Iy都较大
边缘:仅在水平、或者仅在竖直方向有较大的变化量,即Ix和Iy只有其一较大
平坦地区:在水平、竖直方向的变化量均较小,即Ix、Iy都较小
2 strong eigenvalues======interest point
1 strong eigenvalues======contour/edge
0 eigenvalues
角点响应
R=det(M)-k*(trace(M)^2)
det(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"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;ifor(j=size2/2;jfor(i1=0;i1for(j1=0;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;ifor (j=0;jdouble,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 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;ifor(j=size/2;jfor(i1=0;i1for(j1=0;j1if(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;ifor(j=size/2;jif(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;jfor(i=0;i//将灰度图像数值存入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;jfor(i=0;idouble,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;jfor(i=0;idouble,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//定义模板
for(j=0;jint(gausswidth/2))*(i-int(gausswidth/2))+(j-int(gausswidth/2))*(j-int(gausswidth/2)))/(2*sigma));
//归一化:使模板参数之和为1(其实此步可以省略)
double total=0;
for(i=0;ifor(i=0;ifor(j=0;j//进行高斯平滑
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;jfor(i=size/2;iif(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;
}