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