harris corner c的实现代码和一些问题

[mohyu]

先简述一下算法内容
构造一个函数

（x,y）是窗口偏移，I为灰度，（u，v）是像素点。
在点角处，图像窗口的便宜将造成自相关函数E（x，y）的显著变化。

E的泰勒展开



可以吧E表示成


这里的w(x,y) 是高斯平滑滤波函数



特征值可以表示为 R=Det(M)一kTr^2(M) 为什么要表示成如此？

BOOL CDIBDisplayView::Harris(LPSTR lpDIBBits,LONG lWidth, LONG lHeight) 
{
	//gausswidth：二维高斯窗口宽度
	//sigma：高斯函数的方差
	//size：非极大值抑制的邻域宽度
	//thresh：最终确定角点所需的阈值
	int i,j,m,n,size,thresh,gausswidth;
	double sigma;
    int totalcorner = 0;   //角点计数
	//设置四个参数
	 
	 gausswidth =5;
	 sigma =0.8;
	 size =5;
	 thresh =5000;
 
    unsigned char *lpSrc;//一个指向源、目的像素的移动指针

	int cxDIB = (int) lWidth;         // 图像宽度
	int cyDIB = (int) lHeight;        // 图像高度

	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=Template(I,cxDIB,cyDIB,dx,3,3);  //Template做卷积运算

	//定义垂直方向差分算子并求Iy
	double dy[9]={-1,-1,-1,0,0,0,1,1,1};
	Iy=Template(I,cxDIB,cyDIB,dy,3,3);
 
	//计算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));
	
 

	//进行高斯平滑//Template用于卷积
	Ix2=Template(Ix2,cxDIB,cyDIB,g,gausswidth,gausswidth);//A
	Iy2=Template(Iy2,cxDIB,cyDIB,g,gausswidth,gausswidth);//B
	Ixy=Template(Ixy,cxDIB,cyDIB,g,gausswidth,gausswidth);//C
	

  
   /*计算角点量*/
 
	
	//计算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); 
		}
	}

       //注意到这里是特征值的计算,可是程序给出的公式跟论文不符合啊

  /*进行局部非极大值抑制以获得最终角点*/
 
	
	//注意进行局部极大值抑制的思路

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

	//最终确定角点
	//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;   //满足上两个条件，才是角点！
					    totalcorner++;
					}
		}
	}
	//画点
	toal =totalcorner;
	OnImgShowcorner() ;
	delete corner ;
    return true;
}

我的问题是1.特征值可以表示为 R=Det(M)一kTr^2(M) 是如何给出的？

2.程序中的计算方法与论文中看起来不一样。为何？
论文中为： R=Det(M)一kTr^2(M)
程序中为： cim(i,j) = (Ix2(i,j)*Iy2(i,j) - Ixy(i,j)*Ixy(i,j))/(Ix2(i,j) + Iy2(i,j) + 0.000001); （要在分母中加入一个极小量以防止除数为零溢出）

3.为何要做局部极大值抑制，意义何在？

[天鹅梦]

这个比较专业，所以我帮你顶一下

[moonskypxj]

问题1其实论文中已经给出了呀，关于问题2这篇文章给出的：J.A.Noble.Finding corners.Image and Vision Co mputing,1988。问题3，Harris角点最后就是检测局部最大值并大于给定阈值才确认是否是角点，一般这一步是最大值滤波

[fengbingchun]

Harris角点大体意思是一样的，但是实现起来可以根据需要做相应的变动

[mohyu]

[Quote=引用 3 楼 moonskypxj 的回复:]

问题1其实论文中已经给出了呀，关于问题2这篇文章给出的：J.A.Noble.Finding corners.Image and Vision Co mputing,1988。问题3，Harris角点最后就是检测局部最大值并大于给定阈值才确认是否是角点，一般这一步是最大值滤波
[/Quote]

1. 呵呵，我知道论文中有1，我想知道为什么能表示极值，我已经找到

2.那篇文章我找到了，但是。。。看不懂，囧。能不能简单的认为就是（AB-c^2）/(A+B)也是求极值的一种方法

3.这个还在看。

谢谢您的解惑，非常感谢！