特征检测器-Harris角点检测

最新推荐文章于 2023-06-10 10:12:05 发布

奋斗的_小强

最新推荐文章于 2023-06-10 10:12:05 发布

阅读量1.3k

点赞数

分类专栏：视觉算法与应用

本文链接：https://blog.csdn.net/xiaoqiangfighting/article/details/9982309

版权

视觉算法与应用专栏收录该内容

23 篇文章 0 订阅

订阅专栏

Harris角点检测

------------------------------------------------------------------------------------------------------------------------------------------------------------------------

原理：

在我们解决问题时，往往希望找到特征点，“特征”顾名思义，指能描述物体本质的东西，还有一种解释就是这个特征微小的变化都会对物体的某一属性产生重大的影响。而角点就是这样的特征。

角点：最直观的认识就是在水平、竖直两个方向上变化都比较大。理解为拐角。

观察日常生活中的角落就会发现，角落可以视为所有平面的交汇处，或者说是所有表面的发起处。假设我们要改变一个墙角的位置，那么由它而出发的平面势必都要有狠罚的变化。“如果某一点在任意方向的一个微小变动都会引起灰度很大的变化，那么我们就把它称之为角点”

harris是最典型的角点检测算子。角点经常被检测在边缘的交界处、被遮挡的边缘、纹理性很强的部分（不一定非得是拐点？）。满足这些条件一般都是稳定的、重复性比较高的点，所以实际上他们是不是角点并不重要（因为我们的目标就是找一些稳定、重复性高的点作为特征点）。

——————————————————————————————————————————————————————————

算法思路：去检测图像像素的灰度变化情况

求解： $E(u,v)=\sum [I(x+u,y+v)-I(x,y)]^{2}$

其中表示像素的灰度值。对于上式，我们希望找到使E得值尽量大的点（由之上描述可知），则将上式右边进行泰勒展开得到：

$E(u,v)=\sum [I(x,y)+uI_{x}+vI_{y}-I(x,y)]^{2}$ （泰勒展开？？？？？）

整理可得：

$E(u,v)=\sum [u^{2}I_{x}^{2}+v^{2}I_{y}^{2}+2uvI_{x}I_{y}]$

进而可以表示为

$E(u,v)=\sum [(u,v)\begin{pmatrix} I_{x}^{2}& I_{x}I_{y}\\I_{x}I_{y} & I_{y}^{2} \end{pmatrix}\begin{pmatrix} u\\v \end{pmatrix}]$

这里考虑进去窗函数，设

$M=\sum w(x,y) \begin{pmatrix} I_{x}^{2}& I_{x}I_{y}\\I_{x}I_{y} & I_{y}^{2} \end{pmatrix}$

上式得到：

$E(u,v)=(u,v)M\begin{pmatrix} u\\v \end{pmatrix}$

于是，整理出的Harris算子的公式为：
$R=det(M)-k*(trace(M)^2)$ ————————————————————————————————————————————————————————————————————

对于协方差矩阵：协方差矩阵的作用为什么比方差和均值要大？我们知道方差和均值只是一维随机变量的统计值，而协方差不同，它可以表示多维随机变量之间的相关性信息。协方差矩阵的一个很出色的应用就是在PCA中，选择主方向。协方差矩阵的对角线上的元素表示的是各个维度的方差，而非对角线上的元素表示的则是各个维度之间的相关性，因此，在PCA中，我们尽量将非对角线上的元素化为0，即将矩阵对角化。通过选特征值较大的维度，去掉特征值较小的维度来获得主方向，并且使主方向与其他方向的相关性尽量小。对于上面的矩阵M，通过以上对协方差的描述，我们完全可以把这个矩阵看作是一个二维随机分布的协方差矩阵，那么我们要做的就是将其对角化，求矩阵的两个特征值,矩阵的两个特征值可以用来描述两个主要方向上信号的变化，因此特征值可以用来判决是否为特征点，然后根据这两个特征值来判断是不是角点（两个特征值都大则代表角点）。

————————————————————————————————————————————————————————————————————

对Harris算子，我们可以写成： $R=\lambda _{1}\lambda_{2}-k*(\lambda_{1}+\lambda_{2})^2$ 或 $R=det(M)-k*(trace(M)^2)$

那么又为什么上式取值较大时能保证λ1和λ2的取值都很大呢？

a) λ1和λ2一个大而另一个小时，det小而trace大，‘-’号就能使cornerness小（而‘+’号却使cornerness依然很大，所以必须是减号而不是加号)；

b) λ1和λ2都很小时，显然cornerness很小；

c) λ1和λ2都很大时（比参数λ更大），此时det会更大于trace从而使cornerness很大。

可以参考这样一个图：描述了不同纹理下λ1和λ2的取值情况（其中λ1和λ2是矩阵M的两个特征值）：

- 没有什么纹理的情况下，两个值都很小（很小的正值）
- 边缘的点，一个值大，另外一个值小（由于k取了很小的值，所以3.4的结果为一个小负值）
- 角点：两个值都比较大（比较大的正值）

这样，当我们把目标函数定义为3.4式的时候，得到的结果就会尽力满足两个特征值都比较大了。当然，除此之外，还有Harmonic mean等方式实现更理想的组合方式达到检测出的两个特征值都尽可能大。

——————————————————————————————————————————————————————————————————

代码实现：

？？？？？？？？？？如何由R推导到cim???????????????????????????//

一：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

/
#include "cv.h";
#include "highgui.h"

IplImage *image;
#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++)//考虑到滤波模板为3x3的
	  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;

}

//计算Ix^2,Iy^2,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++)
   {
	   //在分母中加入一个极小的量是为了防止除数为0时溢出
    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;
}
int main()
{
	IplImage *src;
	src=cvLoadImage("e:\\lena.bmp",1);
  CvMat *mat_I,*mat_Ix,*mat_Iy,*mat_IxIy,*mat_Ix2,*mat_Iy2;//相应的矩阵
  IplImage *Gray=NULL;
  IplImage *dst=NULL;
  IplImage *imgDx=NULL;//水平梯度卷积后的图像
  IplImage *imgDy=NULL;
  IplImage *imgDx2=NULL;//Ix2图像
  IplImage *imgDy2=NULL;
  IplImage *imgDxy=NULL;

  Gray=cvCreateImage(cvGetSize(src),8,1);
  dst=cvCreateImage(cvGetSize(src),src->depth,3);
  imgDx=cvCreateImage(cvGetSize(src),8,1);
  imgDy=cvCreateImage(cvGetSize(src),8,1);
  imgDx2=cvCreateImage(cvGetSize(src),8,1);
  imgDy2=cvCreateImage(cvGetSize(src),8,1);
  imgDxy=cvCreateImage(cvGetSize(src),8,1);

  const int cxDIB=src->width;//图像的宽度
  const int cyDIB=src->height;//图像的高度
  double *I=new double[cxDIB*cyDIB];
  cvCvtColor(src,Gray,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(Gray,i,j);//将灰度图像的数值存入I中
	  }
	  mat_I=cvCreateMat(cyDIB,cxDIB,CV_64FC1);
	  cvInitMatHeader(mat_I,cyDIB,cxDIB,CV_64FC1,I);//用I来初始化相应的矩阵

	  //——————————————————————————————————————
      //                     第一步：利用差分算子对图像进行滤波
	  //——————————————————————————————————————
      //定义水平方向差分算子并求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矩阵
      
	  for(j=0;j<cyDIB;j++)
	  {
		  for(i=0;i<cxDIB;i++)
		  {
		   //S(imgDx,i,j)=CV_MAT_ELEM(*mat_Ix,double,j,i);//???
		   //S(imgDy,i,j)=CV_MAT_ELEM(*mat_Iy,double,j,i);//为相应图像赋值
			  ((uchar *)(imgDx->imageData+imgDx->widthStep*j))[i]=CV_MAT_ELEM(*mat_Ix,double,j,i);
			  ((uchar *)(imgDy->imageData+imgDy->widthStep*j))[i]=CV_MAT_ELEM(*mat_Iy,double,j,i);
		  }
	  }
	  mat_Ix2=mbxy(mat_Ix,mat_Ix,cxDIB,cyDIB);
	  mat_Iy2=mbxy(mat_Iy,mat_Iy,cxDIB,cyDIB);
	  mat_IxIy=mbxy(mat_Ix,mat_Iy,cxDIB,cyDIB);
	  for(j=0;j<cyDIB;j++)
	  {
		  for(i=0;i<cxDIB;i++)
		  {
		   //S(imgDxy,i,j)=CV_MAT_ELEM(*mat_IxIy,double,j,i);
		   //S(imgDx2,i,j)=CV_MAT_ELEM(*mat_Ix2,double,j,i);
		   //S(imgDy2,i,j)=CV_MAT_ELEM(*mat_Iy2,double,j,i);
			  ((uchar *)(imgDxy->imageData+imgDxy->widthStep*j))[i]=CV_MAT_ELEM(*mat_IxIy,double,j,i);
			  ((uchar *)(imgDx2->imageData+imgDx2->widthStep*j))[i]=CV_MAT_ELEM(*mat_Ix2,double,j,i);
			  ((uchar *)(imgDy2->imageData+imgDy2->widthStep*j))[i]=CV_MAT_ELEM(*mat_Iy2,double,j,i);
		  }
	  }
	  //————————————————————————————————————————
	  //                第二步：对Ix2 Iy2 Ixy 进行高斯平滑，以去除噪声
	  //————————————————————————————————————————
      //本例使用5x5的高斯模板
		  //计算模板参数
		  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++)
				  //考虑到5x5的模板
				  g[i*gausswidth+j]=exp(-((i-int(gausswidth/2))*(i-int(gausswidth/2))+(j-int(gausswidth/2))*(j-int(gausswidth/2)))/(2*sigma));

		  //归一化(可省略)
		  
		  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_IxIy=mbys(mat_IxIy,cxDIB,cyDIB,g,gausswidth,gausswidth);
     //————————————————————————————————————————
		  //     第三步：计算角点量
     //————————————————————————————————————————
     //计算cim:即cornerness of image,称之为角点量
		  CvMat *mat_cim;
		  mat_cim=mbcim(mat_Ix2,mat_Iy2,mat_IxIy,cxDIB,cyDIB);

	  //————————————————————————————————————————
		  //     第四步：进行局部非极大值抑制
	  //————————————————————————————————————————
      CvMat *mat_locmax;
	  const int localsize=7;
	  mat_locmax=mblocmax(mat_cim,cxDIB,cyDIB,localsize);

	  //————————————————————————————————————————
	  //          第五步：获得最终角点
	  //————————————————————————————————————————
      CvMat *mat_corner;
	  const double threshold=4500;
	  int cornernum=0;
	  mat_corner=mbcorner(mat_cim,mat_locmax,cxDIB,cyDIB,gausswidth,threshold);
	  //CvPoint pt[5000];
	  int size=7;
	  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[cornernum].x=i;
			  //pt[cornernum].y=j;
			  cornernum++;
			  //cvCircle(dst,pt[cornernum],10,CV_RGB(255,255,255),5,8,0);??
			  cvCircle(dst,cvPoint(i,j),2,CV_RGB(255,0,0),2,8,0);

			}
		  }
		  std::cout<<"cornernum:"<<cornernum<<std::endl;
		  cvNamedWindow("src",1);
		  cvShowImage("src",src);
		  cvNamedWindow("dst",1);
		  cvShowImage("dst",dst);
		  cvWaitKey(0);
     return 0;

}

三：Matlab

整理自：1、http://www.cnblogs.com/ztfei/archive/2012/05/08/2489900.html

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

奋斗的_小强

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
特征检测器-Harris角点检测

Harris角点检测原理：角点：最直观的认识就是在水平、竖直两个方向上变化都比较大。理解为拐角。 harris是最典型的角点检测算子。角点经常被检测在边缘的交界处、被遮挡的边缘、纹理性很强的部分（不一定非得是拐点？）。满足这些条件一般都是稳定的、重复性比较高的点，所以实际上他们是不是角点并不重要（因为我们的目标就是找一些稳定、重复性高的点作为特征点）。
复制链接

扫一扫