【转】&&Harris角点检测算子
原文地址：
&Harris角点检测算子是于1988年由CHris Harris & Mike
Stephens提出来的。在具体展开之前，不得不提一下Moravec早在1981就提出来的Moravec角点检测算子。
1.Moravec角点检测算子
Moravec角点检测算子的思想其实特别简单，在图像上取一个W*W的“滑动窗口”，不断的移动这个窗口并检测窗口中的像素变化情况E。像素变化情况E可简单分为以下三种：A&
如果在窗口中的图像是什么平坦的，那么E的变化不大。B&
如果在窗口中的图像是一条边，那么在沿这条边滑动时E变化不大，而在沿垂直于这条边的方向滑动窗口时，E的变化会很大。
C& 如果在窗口中的图像是一个角点时，窗口沿任何方向移动E的值都会发生很大变化。
<img WIDTH="268" HEIGHT="199" ALT="" src="/blog7style/images/common/sg_trans.gif" real_src ="http://img.my.csdn.net/uploads//_4764.jpg"
TITLE="【转】&&Harris角点检测算子" />
上图就是对Moravec算子的形象描述。用数学语言来表示的话就是：
<img ALT="" src="/blog7style/images/common/sg_trans.gif" real_src ="http://img.my.csdn.net/uploads//_6755.jpg"
TITLE="【转】&&Harris角点检测算子" />
其中（x，y）就表示四个移动方向（1，0）（1，1）（0，1）（-1，1），E就是像素的变化值。Moravec算子对四个方向进行加权求和来确定变化的大小，然和设定阈值，来确定到底是边还是角点。
2.Harris角点检测算子
Harris角点检测算子实质上就是对Moravec算子的改良和优化。在原文中，作者提出了三点Moravec算子的缺陷并且给出了改良方法：
Moravec算子对方向的依赖性太强，在上文中我们可以看到，Moravec算子实际上只是移动了四个45度角的离散方向，真正优秀的检测算子应该能考虑到各个现象的移动变化情况。为此，作者采用微分的思想（这里不清楚的话可以复习下高数中的全微分）：
<img WIDTH="287" HEIGHT="113" ALT="" src="/blog7style/images/common/sg_trans.gif" real_src ="http://img.my.csdn.net/uploads//_9194.jpg"
TITLE="【转】&&Harris角点检测算子" />
<img WIDTH="221" HEIGHT="81" ALT="" src="/blog7style/images/common/sg_trans.gif" real_src ="http://img.my.csdn.net/uploads//_5324.jpg"
TITLE="【转】&&Harris角点检测算子" />
所以E就可以表示为：
<img ALT="" src="/blog7style/images/common/sg_trans.gif" real_src ="http://img.my.csdn.net/uploads//_7663.jpg"
TITLE="【转】&&Harris角点检测算子" />
2.由于Moravec算子采用的是方形的windows，因此的E的响应比较容易受到干扰，Harris采用了一个较为平滑的窗口——高斯函数：
<img ALT="" src="/blog7style/images/common/sg_trans.gif" real_src ="http://img.my.csdn.net/uploads//_5797.jpg"
TITLE="【转】&&Harris角点检测算子" />
3.Moravec算子对边缘响应过于灵敏。为此，Harris提出了对E进行变形：
<img ALT="" src="/blog7style/images/common/sg_trans.gif" real_src ="http://img.my.csdn.net/uploads//_3245.jpg"
TITLE="【转】&&Harris角点检测算子" />
对，没错，变成了二次型，果然是大牛，更牛的还在后面！其中，
<img ALT="" src="/blog7style/images/common/sg_trans.gif" real_src ="http://img.my.csdn.net/uploads//_9407.jpg"
TITLE="【转】&&Harris角点检测算子" />
用α，β表示矩阵M的特征值，这样会产生三种情况：A&
如果α，β都很小，说明高斯windows中的图像接近平坦。 B& 如果一个大一个小，则表示检测到边。
C& 如果α，β都很大，那么表示检测到了角点。
α，β是什么？α，β就是椭圆的长短轴的度量，什么？椭圆哪里来？椭圆就是那个二次型函数来的！下图的意思我就不详细讲解了，相信大家能懂。
<img WIDTH="416" HEIGHT="302" ALT="" src="/blog7style/images/common/sg_trans.gif" real_src ="http://img.my.csdn.net/uploads//_8342.jpg"
TITLE="【转】&&Harris角点检测算子" />
有人又要问了，你怎么知道我检测到α，β算大还是算小？对此天才哈里斯定义了一个角点响应函数：
<img ALT="" src="/blog7style/images/common/sg_trans.gif" real_src ="http://img.my.csdn.net/uploads//_9598.jpg"
TITLE="【转】&&Harris角点检测算子" />
其中（这些都是线性代数里的知识）：
<img ALT="" src="/blog7style/images/common/sg_trans.gif" real_src ="http://img.my.csdn.net/uploads//_7675.jpg"
TITLE="【转】&&Harris角点检测算子" />
我们惊喜的发现,R为正值是，检测到的是角点，R为负时检测到的是边，R很小时检测到的是平坦区域。至于他怎么想出来的，我就不得而知了......
<img ALT="" src="/blog7style/images/common/sg_trans.gif" real_src ="http://img.my.csdn.net/uploads//_6039.jpg"
TITLE="【转】&&Harris角点检测算子" />
&Harris角点检测算法有诸多优点：A&
旋转不变性，椭圆转过一定角度但是其形状保持不变（特征值保持不变）
<img ALT="" src="/blog7style/images/common/sg_trans.gif" real_src ="http://img.my.csdn.net/uploads//_8431.jpg"
TITLE="【转】&&Harris角点检测算子" />
对于图像灰度的仿射变化具有部分的不变性，由于仅仅使用了图像的一介导数，对于图像灰度平移变化不变；对于图像灰度尺度变化不变
当然Harris也有许多不完善的地方：A& 它对尺度很敏感，不具备几何尺度不变性。
&<img ALT="" src="/blog7style/images/common/sg_trans.gif" real_src ="http://img.my.csdn.net/uploads//_6925.jpg"
TITLE="【转】&&Harris角点检测算子" />
提取的角点是像素级的。以至于后来又有许多牛人提出了更多更完善的检测算子，且听下回分解！
已投稿到：
以上网友发言只代表其个人观点，不代表新浪网的观点或立场。&#xe621; 上传我的文档
&#xe602; 下载
&#xe60c; 收藏
关注毕业设计资料的整理，部分来源与网络。
&#xe602; 下载此文档
正在努力加载中...
改进的Harris算子用于点特征的精确定位
下载积分：300
内容提示：改进的Harris算子用于点特征的精确定位
文档格式：PDF|
浏览次数：0|
上传日期： 12:32:20|
文档星级：&#xe60b;&#xe612;&#xe612;&#xe612;&#xe612;
该用户还上传了这些文档
改进的Harris算子用于点特征的精确定位
官方公共微信/*--------------------CSS部分-------------------*/
/*--------------------JS部分-------------------*/
特征点检测广泛应用到目标匹配、目标跟踪、三维重建等应用中，在进行目标建模时会对图像进行目标特征的提取，常用的有颜色、角点、特征点、轮廓、纹理等特征。现在开始讲解常用的特征点检测，其中Harris角点检测是特征点检测的基础，提出了应用邻近像素点灰度差值概念，从而进行判断是否为角点、边缘、平滑区域。Harris角点检测原理是利用移动的窗口在图像中计算灰度变化值，其中关键流程包括转化为灰度图像、计算差分图像、高斯平滑、计算局部极值、确认角点。
一、基础知识
图像的变化类型：
在特征点检测中经常提出尺度不变、旋转不变、抗噪声影响等，这些是判断特征点是否稳定的指标。
性能较好的角点：
检测出图像中&真实&的角点
准确的定位性能
很高的重复检测率
噪声的鲁棒性
较高的计算效率
角点的类型：
基于图像灰度的方法通过计算点的曲率及梯度来检测角点，避免了第一类方法存在的缺陷，此类方法主要有Moravec算子、Forstner算子、Harris算子、SUSAN算子等。
二、Harris角点原理
1、算法思想
角点原理来源于人对角点的感性判断，即图像在各个方向灰度有明显变化。算法的核心是利用局部窗口在图像上进行移动判断灰度发生较大的变化，所以此窗口用于计算图像的灰度变化为：[-1,0,1;-1,0,1;-1,0,1][-1,-1,-1;0,0,0;1,1,1]。人各个方向上移动这个特征的小窗口，如图3中窗口内区域的灰度发生了较大的变化，那么就认为在窗口内遇到了角点。如图1中，窗口内图像的灰度没有发生变化，那么窗口内就不存在角点；如果窗口在某一个方向移动时，窗口内图像的灰度发生了较大的变化，而在另一些方向上没有发生变化，那么，窗口内的图像可能就是一条直线的线段。
2、数学模型
根据算法思想，构建数学模型，计算移动窗口的的灰度差值。
为了减小计算量，利用泰勒级数进行简化公式：
上图中W函数表示窗口函数，M矩阵为偏导数矩阵。对于矩阵可以进行对称矩阵的变化，假设利用两个特征值进行替代，其几何含义类似下图中的表达。在几何
模型中通过判断两个特征值的大小，来判定像素的属性。
M为梯度的协方差矩阵 ，在实际应用中为了能够应用更好的编程，定义了角点响应函数R，通过判定R大小来判断像素是否为角点。
R取决于M的特征值，对于角点|R|很大，平坦的区域|R|很小，边缘的R为负值。
3、算法流程
1.利用水平，竖直差分算子对图像的每个像素进行滤波以求得Ix,Iy,进而求得M中的四个元素的值。
算法源码：
代码中如果array为-1,0,1,-1,0,1,-1,0,1}则是求解X方向的，如果为{-1,-1,-1,0,0,0,1,1,1}为Y方向的，则Ix和Iy求解结束
//这里的size/2是为了不把图像边界算进去。
//array也就是3*3的窗口函数为：double dx[9]={-1,0,1,-1,0,1,-1,0,1};或者 //定义垂直方向差分算子并求Iy
double dy[9]={-1,-1,-1,0,0,0,1,1,1};
CvMat *mbys(CvMat *mat,int xwidth,int ywidth,double *array,int size1,int size2)//size
int i1,j1;
CvMat *mat1;
mat1=cvCloneMat(mat);
//为了去除边界，从框体一半开始遍历
for(i=size1/2;i&ywidth-size1/2;i++)
for(j=size2/2;j&xwidth-size2/2;j++)
for(i1=0;i1&size1;i1++)
for(j1=0;j1&size2;j1++)
px=i-size1/2+i1;
py=j-size2/2+j1;
//CV_MAT_ELEM访问矩阵元素
//每个元素和窗体函数遍历相加
m+=CV_MAT_ELEM(*mat,double,px,py)*array[i1*size1+j1];
CV_MAT_ELEM(*mat1,double,i,j)=m;
return mat1;
求解IX2相对比较简单，像素相乘即可。下面为基于opencv的C++版本，很是简单
//构建模板
Mat xKernel = (Mat_&double&(1,3) && -1, 0, 1);
Mat yKernel = xKernel.t();
//计算IX和IY
Mat Ix,Iy;
//可参考filter2d的定义
filter2D(gray, Ix, CV_64F, xKernel);
filter2D(gray, Iy, CV_64F, yKernel);
//计算其平方
Mat Ix2,Iy2,I
Ix2 = Ix.mul(Ix);
Iy2 = Iy.mul(Iy);
Ixy = Ix.mul(Iy);
2.对M的四个元素进行高斯平滑滤波，为的是消除一些不必要的孤立点和凸起，得到新的矩阵M。
//本例中使用5&5的高斯模板
//计算模板参数
//int gausswidth=5;
//double sigma=0.8;
double *g=new double[gausswidth*gausswidth];
for(i=0;i&i++)//定义模板
for(j=0;j&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*i++)
total+=g[i];
for(i=0;i&i++)
for(j=0;j&j++)
g[i*gausswidth+j]/=
//进行高斯平滑
mat_Ix2=mbys(mat_Ix2,cxDIB,cyDIB,g,gausswidth,gausswidth);
mat_Iy2=mbys(mat_Iy2,cxDIB,cyDIB,g,gausswidth,gausswidth);
mat_Ixy=mbys(mat_Ixy,cxDIB,cyDIB,g,gausswidth,gausswidth);
利用opencv接口则是很简单
//构建高斯函数
Mat gaussKernel = getGaussianKernel(7, 1);
//卷积计算
filter2D(Ix2, Ix2, CV_64F, gaussKernel);
filter2D(Iy2, Iy2, CV_64F, gaussKernel);
filter2D(Ixy, Ixy, CV_64F, gaussKernel);
3.接下来利用M计算对应每个像素的角点响应函数R，即：
也可以使用改进的R：R=[Ix^2*Iy^2-(Ix*Iy)^2]/(Ix^2+Iy^2);里面没有随意给定的参数k，取值应当比第一个令人满意。
//计算cim：即cornerness of image，我们把它称做&角点量&
CvMat *mat_
mat_cim=mbcim(mat_Ix2,mat_Iy2,mat_Ixy,cxDIB,cyDIB);
//用来求得响应度
CvMat *mbcim(CvMat *mat1,CvMat *mat2,CvMat *mat3,int xwidth,int ywidth)
mat=cvCloneMat(mat1);
for(i = 0; i & i++)
for(j = 0; j & j++)
//注意：要在分母中加入一个极小量以防止除数为零溢出
CV_MAT_ELEM(*mat,double,i,j)=(CV_MAT_ELEM(*mat1,double,i,j)*CV_MAT_ELEM(*mat2,double,i,j)-
CV_MAT_ELEM(*mat3,double,i,j)*CV_MAT_ELEM(*mat3,double,i,j))/
(CV_MAT_ELEM(*mat1,double,i,j)+CV_MAT_ELEM(*mat2,double,i,j)+0.000001);
//opencv代码
Mat cornerStrength(gray.size(), gray.type());
for (int i = 0; i & gray. i++)
for (int j = 0; j & gray. j++)
double det_m = Ix2.at&double&(i,j) * Iy2.at&double&(i,j) - Ixy.at&double&(i,j) * Ixy.at&double&(i,j);
double trace_m = Ix2.at&double&(i,j) + Iy2.at&double&(i,j);
cornerStrength.at&double&(i,j) = det_m - alpha * trace_m *trace_m;
4、局部极大值抑制，同时选取其极大值
//--------------------------------------------------------------------------
第四步：进行局部非极大值抑制
//--------------------------------------------------------------------------
CvMat *mat_
//const int size=7;
mat_locmax=mblocmax(mat_cim,cxDIB,cyDIB,size);
cout&&CV_MAT_ELEM(*mat_locmax,double,cyDIB-1,cxDIB-1)&&
//用来求得局部极大值
CvMat *mblocmax(CvMat *mat1,int xwidth,int ywidth,int size)
double max=-1000;
int i1,j1;
mat=cvCloneMat(mat1);
for(i=size/2;i&ywidth-size/2;i++)
for(j=size/2;j&xwidth-size/2;j++)
max=-10000;
for(i1=0;i1&i1++)
for(j1=0;j1&j1++)
px=i-size/2+i1;
py=j-size/2+j1;
if(CV_MAT_ELEM(*mat1,double,px,py)&max)
max=CV_MAT_ELEM(*mat1,double,px,py);
CV_MAT_ELEM(*mat,double,i,j)=
CV_MAT_ELEM(*mat,double,i,j)=0;
在opencv中应用maxminloc函数
// threshold
double maxS
minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);
Mat localM
//默认3*3核膨胀，膨胀之后，除了局部最大值点和原来相同，其它非局部最大值点被
//3*3邻域内的最大值点取代
dilate(cornerStrength, dilated, Mat());
//与原图相比，只剩下和原图值相同的点，这些点都是局部最大值点，保存到localMax
compare(cornerStrength, dilated, localMax, CMP_EQ);
5.在矩阵R中，同时满足R(i,j)大于一定阈值threshold和R(i,j)是某领域内的局部极大值，则被认为是角点。
cv::Mat cornerM
// 根据角点响应最大值计算阈值
threshold= qualityLevel*maxS
cv::threshold(cornerStrength,cornerTh,
threshold,255,cv::THRESH_BINARY);
// 转为8-bit图
cornerTh.convertTo(cornerMap,CV_8U);// 和局部最大值图与，剩下角点局部最大值图，即：完成非最大值抑制
cv::bitwise_and(cornerMap,localMax,cornerMap);
imgDst = cornerMap.clone();
for( int y = 0; y & cornerMap. y++ )
const uchar* rowPtr = cornerMap.ptr&uchar&(y);
for( int x = 0; x & cornerMap. x++ )
// 非零点就是角点
if (rowPtr[x])
points.push_back(cv::Point(x,y));
三、算法源码
里面函数最基础的代码解释和注释
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
#define B(image,x,y) ((uchar *)(image-&imageData+image-&widthStep*(y)))[(x)*3]
#define G(image,x,y) ((uchar *)(image-&imageData+image-&widthStep*(y)))[(x)*3+1]
#define R(image,x,y) ((uchar *)(image-&imageData+image-&widthStep*(y)))[(x)*3+2]
#define S(image,x,y) ((uchar *)(image-&imageData+image-&widthStep*(y)))[(x)]
//卷积计算求Ix,Iy,以及滤波
//a指向的数组是size1*size2大小的...求导
CvMat *mbys(CvMat *mat,int xwidth,int ywidth,double *a,int size1,int size2)
int i1,j1;
CvMat *mat1;
mat1=cvCloneMat(mat);
for(i=size1/2;i&ywidth-size1/2;i++)
for(j=size2/2;j&xwidth-size2/2;j++)
for(i1=0;i1&size1;i1++)
for(j1=0;j1&size2;j1++)
px=i-size1/2+i1;
py=j-size2/2+j1;
//CV_MAT_ELEM访问矩阵元素
m+=CV_MAT_ELEM(*mat,double,px,py)*a[i1*size1+j1];
CV_MAT_ELEM(*mat1,double,i,j)=m;
return mat1;
//计算Ix2,Iy2,Ixy
CvMat *mbxy(CvMat *mat1,CvMat *mat2,int xwidth,int ywidth)
CvMat *mat3;
mat3=cvCloneMat(mat1);
for(i=0;i&i++)
for (j=0;j&j++)
CV_MAT_ELEM(*mat3,double,i,j)=CV_MAT_ELEM(*mat1,double,i,j)*CV_MAT_ELEM(*mat2,double,i,j);
return mat3;
//用来求得响应度
CvMat *mbcim(CvMat *mat1,CvMat *mat2,CvMat *mat3,int xwidth,int ywidth)
mat=cvCloneMat(mat1);
for(i = 0; i & i++)
for(j = 0; j & j++)
//注意：要在分母中加入一个极小量以防止除数为零溢出
CV_MAT_ELEM(*mat,double,i,j)=(CV_MAT_ELEM(*mat1,double,i,j)*CV_MAT_ELEM(*mat2,double,i,j)-
CV_MAT_ELEM(*mat3,double,i,j)*CV_MAT_ELEM(*mat3,double,i,j))/
(CV_MAT_ELEM(*mat1,double,i,j)+CV_MAT_ELEM(*mat2,double,i,j)+0.000001);
//用来求得局部极大值
CvMat *mblocmax(CvMat *mat1,int xwidth,int ywidth,int size)
double max=-1000;
int i1,j1;
mat=cvCloneMat(mat1);
for(i=size/2;i&ywidth-size/2;i++)
for(j=size/2;j&xwidth-size/2;j++)
max=-10000;
for(i1=0;i1&i1++)
for(j1=0;j1&j1++)
px=i-size/2+i1;
py=j-size/2+j1;
if(CV_MAT_ELEM(*mat1,double,px,py)&max)
max=CV_MAT_ELEM(*mat1,double,px,py);
CV_MAT_ELEM(*mat,double,i,j)=
CV_MAT_ELEM(*mat,double,i,j)=0;
//用来确认角点
CvMat *mbcorner(CvMat *mat1,CvMat *mat2,int xwidth,int ywidth,int size,double thresh)
mat=cvCreateMat(ywidth,xwidth,CV_32FC1);
for(i=size/2;i&ywidth-size/2;i++)
for(j=size/2;j&xwidth-size/2;j++)
if(CV_MAT_ELEM(*mat1,double,i,j)==CV_MAT_ELEM(*mat2,double,i,j))//首先取得局部极大值
if(CV_MAT_ELEM(*mat1,double,i,j)&thresh)//然后大于这个阈值
CV_MAT_ELEM(*mat,int,i,j)=255;//满足上两个条件，才是角点！
CV_MAT_ELEM(*mat,int,i,j)=0;
CvPoint* CHarris::harris_features(IplImage *src,int gausswidth,double sigma,int size,int threshold)
CvMat *mat_I,*mat_Ix,*mat_Iy,*mat_Ixy,*mat_Ix2,*mat_Iy2;//相应的矩阵
IplImage *pImgGray=NULL;
//灰度图像
IplImage *dst=NULL;
//目标图像
IplImage *pImgDx=NULL; //水平梯度卷积后的图像
IplImage *pImgDy=NULL; //竖起梯度卷积后的图像
IplImage *pImgDx2=NULL;//Ix2图像
IplImage *pImgDy2=NULL;//Iy2图像
IplImage *pImgDxy=NULL;//Ixy图像
pImgGray=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);
dst=cvCreateImage(cvGetSize(src),src-&depth,3);
pImgDx=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);//创建图像
pImgDy=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);
pImgDx2=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);
pImgDy2=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);
pImgDxy=cvCreateImage(cvGetSize(src),IPL_DEPTH_8U,1);
const int cxDIB=src-&
// 图像宽度
const int cyDIB=src-&
// 图像高度
double *I=new double[cxDIB*cyDIB];
cvCvtColor(src,pImgGray,CV_RGB2GRAY);//灰度化
dst=cvCloneImage(src);
for(j=0;j&cyDIB;j++)
for(i=0;i&cxDIB;i++)
I[j*cxDIB+i]=S(pImgGray,i,j);//将灰度图像数值存入I中
mat_I=cvCreateMat(cyDIB,cxDIB,CV_64FC1);
cvInitMatHeader(mat_I,cyDIB,cxDIB,CV_64FC1,I);//用I来初始化相应的矩阵
cout&&CV_MAT_ELEM(*mat_I,double,200,200)&&
//--------------------------------------------------------------------------
第一步：利用差分算子对图像进行滤波
//--------------------------------------------------------------------------
//定义水平方向差分算子并求Ix
double dx[9]={-1,0,1,-1,0,1,-1,0,1};
mat_Ix=mbys(mat_I,cxDIB,cyDIB,dx,3,3); //求Ix矩阵
cout&&CV_MAT_ELEM(*mat_Ix,double,200,200)&&
//定义垂直方向差分算子并求Iy
double dy[9]={-1,-1,-1,0,0,0,1,1,1};
mat_Iy=mbys(mat_I,cxDIB,cyDIB,dy,3,3);//求Iy矩阵
cout&&CV_MAT_ELEM(*mat_Iy,double,200,200)&&
for(j=0;j&cyDIB;j++)
for(i=0;i&cxDIB;i++)
S(pImgDx,i,j)=CV_MAT_ELEM(*mat_Ix,double,j,i);//为相应图像赋值
S(pImgDy,i,j)=CV_MAT_ELEM(*mat_Iy,double,j,i);
mat_Ix2=mbxy(mat_Ix,mat_Ix,cxDIB,cyDIB);//计算Ix2,Iy2,Ixy矩阵
mat_Iy2=mbxy(mat_Iy,mat_Iy,cxDIB,cyDIB);
mat_Ixy=mbxy(mat_Ix,mat_Iy,cxDIB,cyDIB);
for(j=0;j&cyDIB;j++)
for(i=0;i&cxDIB;i++)
S(pImgDxy,i,j)=CV_MAT_ELEM(*mat_Ixy,double,j,i);//为相应图像赋值
S(pImgDx2,i,j)=CV_MAT_ELEM(*mat_Ix2,double,j,i);
S(pImgDy2,i,j)=CV_MAT_ELEM(*mat_Iy2,double,j,i);
//--------------------------------------------------------------------------
第二步：对Ix2/Iy2/Ixy进行高斯平滑，以去除噪声
//--------------------------------------------------------------------------
//本例中使用5&5的高斯模板
//计算模板参数
//int gausswidth=5;
//double sigma=0.8;
double *g=new double[gausswidth*gausswidth];
for(i=0;i&i++)//定义模板
for(j=0;j&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*i++)
total+=g[i];
for(i=0;i&i++)
for(j=0;j&j++)
g[i*gausswidth+j]/=
//进行高斯平滑
mat_Ix2=mbys(mat_Ix2,cxDIB,cyDIB,g,gausswidth,gausswidth);
mat_Iy2=mbys(mat_Iy2,cxDIB,cyDIB,g,gausswidth,gausswidth);
mat_Ixy=mbys(mat_Ixy,cxDIB,cyDIB,g,gausswidth,gausswidth);
//--------------------------------------------------------------------------
第三步：计算角点量
//--------------------------------------------------------------------------
//计算cim：即cornerness of image，我们把它称做&角点量&
CvMat *mat_
mat_cim=mbcim(mat_Ix2,mat_Iy2,mat_Ixy,cxDIB,cyDIB);
cout&&CV_MAT_ELEM(*mat_cim,double,cyDIB-1,cxDIB-1)&&
//--------------------------------------------------------------------------
第四步：进行局部非极大值抑制
//--------------------------------------------------------------------------
CvMat *mat_
//const int size=7;
mat_locmax=mblocmax(mat_cim,cxDIB,cyDIB,size);
cout&&CV_MAT_ELEM(*mat_locmax,double,cyDIB-1,cxDIB-1)&&
//--------------------------------------------------------------------------
第五步：获得最终角点
//--------------------------------------------------------------------------
CvMat *mat_
//const double threshold=4500;
//int cornernum=0;
mat_corner=mbcorner(mat_cim,mat_locmax,cxDIB,cyDIB,gausswidth,threshold);
//CCommon CommonC
CvPoint pt[5000];
for(j=size/2;j&cyDIB-size/2;j++)
for(i=size/2;i&cxDIB-size/2;i++)
if(CV_MAT_ELEM(*mat_corner,int,j,i)==255)
pt[cornerno].x=i;
pt[cornerno].y=j;
cornerno++;
CommonClass.DrawCross(showImg2,pt,CV_RGB(0,0,255),1,4);
cvCircle(dst,pt,2,CV_RGB(255,0,0),1,8,0);
cout&&i&&" "&&j&&
2、基于opencv源码
1 #include "imgFeat.h"
3 void feat::detectHarrisCorners(const Mat& imgSrc, Mat& imgDst, double alpha)
if (imgSrc.channels() == 3)
cvtColor(imgSrc, gray, CV_BGR2GRAY);
gray = imgSrc.clone();
gray.convertTo(gray, CV_64F);
Mat xKernel = (Mat_&double&(1,3) && -1, 0, 1);
Mat yKernel = xKernel.t();
Mat Ix,Iy;
filter2D(gray, Ix, CV_64F, xKernel);
filter2D(gray, Iy, CV_64F, yKernel);
Mat Ix2,Iy2,I
Ix2 = Ix.mul(Ix);
Iy2 = Iy.mul(Iy);
Ixy = Ix.mul(Iy);
Mat gaussKernel = getGaussianKernel(7, 1);
filter2D(Ix2, Ix2, CV_64F, gaussKernel);
filter2D(Iy2, Iy2, CV_64F, gaussKernel);
filter2D(Ixy, Ixy, CV_64F, gaussKernel);
Mat cornerStrength(gray.size(), gray.type());
for (int i = 0; i & gray. i++)
for (int j = 0; j & gray. j++)
double det_m = Ix2.at&double&(i,j) * Iy2.at&double&(i,j) - Ixy.at&double&(i,j) * Ixy.at&double&(i,j);
double trace_m = Ix2.at&double&(i,j) + Iy2.at&double&(i,j);
cornerStrength.at&double&(i,j) = det_m - alpha * trace_m *trace_m;
// threshold
double maxS
minMaxLoc(cornerStrength, NULL, &maxStrength, NULL, NULL);
Mat localM
dilate(cornerStrength, dilated, Mat());
compare(cornerStrength, dilated, localMax, CMP_EQ);
Mat cornerM
double qualityLevel = 0.01;
double thresh = qualityLevel * maxS
cornerMap = cornerStrength &
bitwise_and(cornerMap, localMax, cornerMap);
imgDst = cornerMap.clone();
63 void feat::drawCornerOnImage(Mat& image, const Mat&binary)
Mat_&uchar&::const_iterator it = binary.begin&uchar&();
Mat_&uchar&::const_iterator itd = binary.end&uchar&();
for (int i = 0; it != it++, i++)
circle(image, Point(i%image.cols, i / image.cols), 3, Scalar(0, 255, 0), 1);
opencv+C++推荐
3、opencv中Harris接口
OpenCV的Hairrs角点检测的函数为cornerHairrs()，但是它的输出是一幅浮点值图像，浮点值越高，表明越可能是特征角点，我们需要对图像进行阈值化。
C++: void cornerHarris(InputArray src, OutputArray dst, int blockSize, int apertureSize, double k, int borderType = BORDER_DEFAULT);
src & 输入的单通道8-bit或浮点图像。
dst & 存储着Harris角点响应的图像矩阵，大小与输入图像大小相同，是一个浮点型矩阵。
blockSize & 邻域大小。
apertureSize & 扩展的微分算子大。
k & 响应公式中的，参数$\alpha$。
boderType & 边界处理的类型。
int main()
Mat image = imread("../buliding.png");
cvtColor(image, gray, CV_BGR2GRAY);
Mat cornerS
cornerHarris(gray, cornerStrength, 3, 3, 0.01);
threshold(cornerStrength, cornerStrength, 0.001, 255, THRESH_BINARY);
从上面上间一幅图像我们可以看到，有很多角点都是粘连在一起的，我们下面通过加入非极大值抑制来进一步去除一些粘在一起的角点。
非极大值抑制原理是，在一个窗口内，如果有多个角点则用值最大的那个角点，其他的角点都删除，窗口大小这里我们用3*3，程序中通过图像的膨胀运算来达到检测极大值的目的，因为默认参数的膨胀运算就是用窗口内的最大值替代当前的灰度值。
int main()
Mat image = imread("buliding.png");
cvtColor(image, gray, CV_BGR2GRAY);
Mat cornerS
cornerHarris(gray, cornerStrength, 3, 3, 0.01);
double maxS
double minS
// 找到图像中的最大、最小值
minMaxLoc(cornerStrength, &minStrength, &maxStrength);
// 膨胀图像，最找出图像中全部的局部最大值点
dilate(cornerStrength, dilated, Mat());
// compare是一个逻辑比较函数，返回两幅图像中对应点相同的二值图像
compare(cornerStrength, dilated, locaMax, CMP_EQ);
Mat cornerM
double qualityLevel = 0.01;
double th = qualityLevel*maxS // 阈值计算
threshold(cornerStrength, cornerMap, th, 255, THRESH_BINARY);
cornerMap.convertTo(cornerMap, CV_8U);
// 逐点的位运算黑色图片显示的结果
bitwise_and(cornerMap, locaMax, cornerMap);
drawCornerOnImage(image, cornerMap);
namedWindow("result");
imshow("result", image);
waitKey();
void drawCornerOnImage(Mat& image, const Mat&binary)
Mat_&uchar&::const_iterator it = binary.begin&uchar&();
Mat_&uchar&::const_iterator itd = binary.end&uchar&();
for (int i = 0; it != it++, i++)
circle(image, Point(i%image.cols, i / image.cols), 3, Scalar(0, 255, 0), 1);
四、Harris角点性质
1、阈值决定检测点数量
增大$\alpha$的值，将减小角点响应值$R$，降低角点检测的灵性，减少被检测角点的数量；减小$\alpha$值，将增大角点响应值$R$，增加角点检测的灵敏性，增加被检测角点的数量
2、Harris角点检测算子对亮度和对比度的变化不敏感
这是因为在进行Harris角点检测时，使用了微分算子对图像进行微分运算，而微分运算对图像密度的拉升或收缩和对亮度的抬高或下降不敏感。换言之，对亮度和对比度的仿射变换并不改变Harris响应的极值点出现的位置，但是，由于阈值的选择，可能会影响角点检测的数量。
3. Harris角点检测算子具有旋转不变性
Harris角点检测算子使用的是角点附近的区域灰度二阶矩矩阵。而二阶矩矩阵可以表示成一个椭圆，椭圆的长短轴正是二阶矩矩阵特征值平方根的倒数。当特征椭圆转动时，特征值并不发生变化，所以判断角点响应值也不发生变化，由此说明Harris角点检测算子具有旋转不变性。
4. Harris角点检测算子不具有尺度不变性
如下图所示，当右图被缩小时，在检测窗口尺寸不变的前提下，在窗口内所包含图像的内容是完全不同的。左侧的图像可能被检测为边缘或曲线，而右侧的图像则可能被检测为一个角点。
五、函数备注
1、compare
功能：两个数组之间或者一个数组和一个常数之间的比较
void compare(InputArray src1, InputArray src2, OutputArray dst, int cmpop)
：第一个数组或者标量，如果是数组，必须是单通道数组。src2 ：第二个数组或者标量，如果是数组，必须是单通道数组。dst
：输出数组，和输入数组有同样的size和type=CV_8UC1cmpop ：
标志指明了元素之间的对比关系CMP_EQ src1&相等&src2.
CMP_GT src1 大于 src2.CMP_GE src1 大于或等于 src2.CMP_LT src1 小于
src2.CMP_LE src1 小于或等于 src2.CMP_NE src1 不等于 src2.
如果对比结果为true，那么输出数组对应元素的值为255，否则为0
//获取角点图
Mat getCornerMap(double qualityLevel) {
Mat cornerM
// 根据角点响应最大值计算阈值
thresholdvalue= qualityLevel*maxS
threshold(cornerStrength,cornerTh,
thresholdvalue,255,cv::THRESH_BINARY);
// 转为8-bit图
cornerTh.convertTo(cornerMap,CV_8U);
// 和局部最大值图与，剩下角点局部最大值图，即：完成非最大值抑制
bitwise_and(cornerMap,localMax,cornerMap);
return cornerM
2、bitwise_and（位运算函数）
功能：计算两个数组或数组和常量之间与的关系
void bitwise_and(InputArray src1, InputArray src2, OutputArray dst, InputArray mask=noArray())
：第一个输入的数组或常量src2 ：第二个输入的数组或常量dst ：输出数组，和输入数组有同样的size和typemask
：可选择的操作掩码，为8位单通道数组，指定了输出数组哪些元素可以被改变，哪些不可以
操作过程为：
如果为多通道数组，每个通道单独处理
3、filter2D
OpenCV中的内部函数filter2D可以直接用来做图像卷积操作
void filter2D(InputArray src, OutputArray dst, int ddepth, InputArray kernel, Point anchor=Point(-1,-1), double delta=0, int borderType=BORDER_DEFAULT )
六、参考文章&
阅读(...) 评论()}