Harris算子 Visual C++算法

Harris算法原理： 角点：最直观的印象就是在水平、竖直两个方向上变化均较大的点，即Ix、Iy都较大 边缘：仅在水平、或者仅在竖直方向有较大的变化量，即Ix和Iy只有其一较大 平坦地区：在水平、竖直方向的变化量均较小，即Ix、Iy都较小 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); }