vc实现harris角点检测算法

   Harris算法是最常用的图像角点检测算法之一。本文作者将讲述Harris算法的直观意义及检测方法，并给出VC实现源代码。

源代码可以去我的网络硬盘下载:

http://contact.ys168.com/

Harris算法原理：

 

 

 

 

 

 

    角点：最直观的印象就是在水平、竖直两个方向上变化均较大的点，即Ix、Iy都较大

 

 

 

    边缘：仅在水平、或者仅在竖直方向有较大的变化量，即Ix和Iy只有其一较大

 

 

 

    平坦地区：在水平、竖直方向的变化量均较小，即Ix、Iy都较小

 

 

 

    接下来我们来分析更一般的情况，得出Harris算法的判定公式。因为本算法公式不好打，我直接把我的文章作了一个抓图。

---

VC编程实现：

void CImagetestDoc::OnImgHarris()
{
 //gausswidth：二维高斯窗口宽度
 //sigma：高斯函数的方差
 //size：非极大值抑制的邻域宽度
 //thresh：最终确定角点所需的阈值
 int i,j,m,n,size,thresh,gausswidth;
 double sigma;

 //输入四个参数
 CInput2 input;
 input.m_gausswidth =5;
 input.m_sigma =0.8;
 input.m_size =5;
 input.m_thresh =5000;
 input.DoModal ();
 gausswidth=input.m_gausswidth ;
 sigma=input.m_sigma ;
 size=input.m_size ;
 thresh=input.m_thresh ;

    unsigned char *lpSrc;//一个指向源、目的像素的移动指针
 LPSTR lpDIB = (LPSTR) ::GlobalLock((HGLOBAL)m_hDIB);
 int cxDIB = (int) ::DIBWidth(lpDIB);         // 图像宽度
 int cyDIB = (int) ::DIBHeight(lpDIB);        // 图像高度
 LPSTR lpDIBBits=::FindDIBBits (lpDIB);
 long lLineBytes = WIDTHBYTES(cxDIB * 8);     // 计算灰度图像每行的字节数

 //创建I、Ix、Ix2、Iy、Iy2、Ixy、cim、mx、corner数组
 double *I=new double[cxDIB*cyDIB];
 double *Ix=new double[cxDIB*cyDIB];
 double *Ix2=new double[cxDIB*cyDIB];
 double *Iy=new double[cxDIB*cyDIB];
 double *Iy2=new double[cxDIB*cyDIB];
 double *Ixy=new double[cxDIB*cyDIB];
 double *cim=new double[cxDIB*cyDIB];
 double *mx=new double[cxDIB*cyDIB];

 corner=new bool[cxDIB*cyDIB];
 memset(corner, 0, cxDIB*cyDIB*sizeof(bool));

 //定义宏以方便访问元素
 #define I(ROW,COL) I[cxDIB*(ROW)+(COL)]
 #define Ix(ROW,COL) Ix[cxDIB*(ROW)+(COL)]
 #define Ix2(ROW,COL) Ix2[cxDIB*(ROW)+(COL)]
 #define Iy(ROW,COL) Iy[cxDIB*(ROW)+(COL)]
 #define Iy2(ROW,COL) Iy2[cxDIB*(ROW)+(COL)]
 #define Ixy(ROW,COL) Ixy[cxDIB*(ROW)+(COL)]
 #define cim(ROW,COL) cim[cxDIB*(ROW)+(COL)]
 #define mx(ROW,COL) mx[cxDIB*(ROW)+(COL)]
 #define corner(ROW,COL) corner[cxDIB*(ROW)+(COL)]

 //将图像灰度值复制到I中，这步很重要！想想为什么？
 for(i = 0; i < cyDIB; i++)
 {
  for(j = 0; j < cxDIB; j++)
  {   
   lpSrc = (unsigned char*)lpDIBBits + lLineBytes * (cyDIB - 1 - i) + j;
   //将256级灰度图像转化为double型
   I(i,j)=double(*lpSrc);
  }
 }

 //--------------------------------------------------------------------------
 //                     第一步：利用差分算子对图像进行滤波
 //--------------------------------------------------------------------------
 
 //定义水平方向差分算子并求Ix
 double dx[9]={-1,0,1,-1,0,1,-1,0,1};
 Ix=mbys(I,cxDIB,cyDIB,dx,3,3);

 //定义垂直方向差分算子并求Iy
 double dy[9]={-1,-1,-1,0,0,0,1,1,1};
 Iy=mbys(I,cxDIB,cyDIB,dy,3,3);

 
 //将中间结果Ix写入到文本文件以便后续分析
 FILE *fp;
 fp=fopen("Ix.txt","w+");
 for(i = 0; i < cyDIB; i++)
 {
  for(j = 0; j < cxDIB; j++)
   fprintf(fp,"%f ",Ix(i,j));
  fprintf(fp,"/n");
 }

 fp=fopen("Iy.txt","w+");
 for(i = 0; i < cyDIB; i++)
 {
  for(j = 0; j < cxDIB; j++)
   fprintf(fp,"%f ",Iy(i,j));
  fprintf(fp,"/n");
 }
 //计算Ix2、Iy2、Ixy
 for(i = 0; i < cyDIB; i++)
 {
  for(j = 0; j < cxDIB; j++)
  { Ix2(i,j)=Ix(i,j)*Ix(i,j);
   Iy2(i,j)=Iy(i,j)*Iy(i,j);
   Ixy(i,j)=Ix(i,j)*Iy(i,j);
  }
 }

 //--------------------------------------------------------------------------
 //                  第二步：对Ix2/Iy2/Ixy进行高斯平滑，以去除噪声
 //--------------------------------------------------------------------------
 
 //本例中使用5×5的高斯模板
 //计算模板参数
 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;
 
 fp=fopen("g.txt","w+");
 for(i = 0; i < gausswidth; i++)
 {
  for(j = 0; j < gausswidth; j++)
   fprintf(fp,"%f ",g[i*gausswidth+j]);
  fprintf(fp,"/n");
 }
 //进行高斯平滑
 Ix2=mbys(Ix2,cxDIB,cyDIB,g,gausswidth,gausswidth);
 Iy2=mbys(Iy2,cxDIB,cyDIB,g,gausswidth,gausswidth);
 Ixy=mbys(Ixy,cxDIB,cyDIB,g,gausswidth,gausswidth);
 

 //--------------------------------------------------------------------------
 //                        第三步：计算角点量
 //--------------------------------------------------------------------------
 
 //计算cim：即cornerness of image，我们把它称做‘角点量’
 for(i = 0; i < cyDIB; i++)
 {
  for(j = 0; j < cxDIB; j++)
  {
   //注意：要在分母中加入一个极小量以防止除数为零溢出
   cim(i,j) = (Ix2(i,j)*Iy2(i,j) - Ixy(i,j)*Ixy(i,j))/(Ix2(i,j) + Iy2(i,j) + 0.000001);
  }
 }

 fp=fopen("cim.txt","w+");
 for(i = 0; i < cyDIB; i++)
 {
  for(j = 0; j < cxDIB; j++)
   fprintf(fp,"%f ",cim(i,j));
  fprintf(fp,"/n");
 }
 //--------------------------------------------------------------------------
 //                 第四步：进行局部非极大值抑制以获得最终角点
 //--------------------------------------------------------------------------
 
 //注意进行局部极大值抑制的思路

 //const double size=7;
 double max;
 //对每个点在邻域内做极大值滤波：即将该点的值设为邻域中最大的那个值（跟中值滤波有点类似）
 for(i = 0; i < cyDIB; i++)
 {
  for(j = 0; j < cxDIB; j++)
  {
   max=-1000000;
   if(i>int(size/2) && i<cyDIB-int(size/2) && j>int(size/2) && j<cxDIB-int(size/2))
    for(m=0;m<size;m++)
    {
     for(n=0;n<size;n++)
      {      
       if(cim(i+m-int(size/2),j+n-int(size/2))>max)
        max=cim(i+m-int(size/2),j+n-int(size/2));
      }
    }
   if(max>0)
    mx(i,j)=max;
   else
    mx(i,j)=0;
  }
 }
 
 fp=fopen("mx.txt","w+");
 for(i = 0; i < cyDIB; i++)
 {
  for(j = 0; j < cxDIB; j++)
   fprintf(fp,"%f ",mx(i,j));
  fprintf(fp,"/n");
 }
 //最终确定角点
 //const double thresh=4500;
 for(i = 0; i < cyDIB; i++)
 {
  for(j = 0; j < cxDIB; j++)
  { 
   if(cim(i,j)==mx(i,j))  //首先取得局部极大值
     if(mx(i,j)>thresh)  //然后大于这个阈值
      corner(i,j)=1;   //满足上两个条件，才是角点！
  }
 }
 
 fp=fopen("corner.txt","w+");
 for(i = 0; i < cyDIB; i++)
 {
  for(j = 0; j < cxDIB; j++)
   fprintf(fp,"%d ",corner(i,j));
  fprintf(fp,"/n");
 }

 ::GlobalUnlock((HGLOBAL) m_hDIB);  
    UpdateAllViews(NULL, 0, NULL); 
}