harris角点检测算法实现VC

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);  }