Harris角点检测算法编程步骤及示例演示

1、Harris角点检测算法编程步骤及示例演示也不说那么多废话了，要介绍啥背景意义之类的，角点检测，顾名思义，就是检测角点，最简单的就是两条线的交点了，还有比如下国际象棋的棋盘格子的交点之类的，反正就是检测这些点。简单将Harris角点检测算法的思想说下，就是拿一个小窗在图像中移动，通过考察这个小窗口内图像灰度的平均变换值来确定角点。(1)如果窗口内区域图像的灰度值恒定，那么所有不同方向的偏移几乎不发生变化；(2)如果窗口跨越一条边，那么沿着这条边的偏移几乎不发生变化，但是与边垂直的偏移会发生很大的变化；(3)如果窗口包含一个孤立的点或者角点，那么所有不同方向的偏移会发生很大的变化。下面给出具体数

2、学推导：设图像窗口平移量为(u,v),产生的灰度变化为E(u,v),有E(u,v尸sumw(x,y)I(x+u,y+v)-I(x,y)八2/中w(x,y)为窗口函数，I(x+u,y+v为平移后的灰度值，I(x,y)为平移前的灰度值。有泰勒公式展开可得：I(x+u,y+v)=I(x,y)+Ix*u+Iy*v+O(uA2,vA2);Ix,Iy分别为偏微分，在图像中为图像的方向导数.因此E(u,v)=sumw(x,y)Ix*u+Iy*v+O(uA2,vA2)A2,可以近似得到E(u,v)=sumw(x,y)Ix*u+Iy*vA2,即E(u,v)=u,vIxA2,Ix*Iy;Ix*Iy,IyA2u,v

3、T令M=IxA2,Ix*Iy;Ix*Iy,IyA2,因此最后对角点的检测成了对矩阵M的特征值的分析了，令M其特征值为x1,x2;当x1>>x2或者x2>>x1,则检测到的是边缘部分；当x1,x2都很小，图像窗口在所有移动的方向上移动灰度级都无明显变化.当X1,X2者B很大时且相当，检测到的是角点。编程时用x1,x2不方便，因此定义角点响应函数；R=det(M)-k(trace(M)八2;其中det(M)为矩阵M的行列式，trace(M)为矩阵M的迹。下面给出更具数学公式实际编程的步骤：1 .利用水平，竖直差分算子对图像的每个像素进行滤波以求得Ix,Iy进而求得M中的四个

4、元素的值。M=IxA2,Ix*Iy;Ix*Iy,IyA22 .对M的四个元素进行高斯平滑滤波，为的是消除一些不必要的孤立点和凸起，得到新的矩阵M。3 .接下来利用M计算对应每个像素的角点响应函数R,即：R=det(M)-k(trace(M)A2;也可以使用改进的R：R=IxA2*IyA2-(Ix*Iy)A2/(IxA2+IyA2);里面没有随意给定的参数k，取值应当比第一个令人满意。4 .在矩阵R中，同时满足R(i,j庆于一定阈值threshold和R(i,j)是某领域内的局部极大值，则被认为是角点。下面给出程序代码：还不错，自己一步步查资料写的，还算详细。functionmainclc;cl

5、ear;closeall;%TheconerpointsextractionsusingHarrismethods;frame=imread('*.jpg');figure(1);imshow(frame);%调用harris角点检测子函数;%输入参数解释：frame为输入图像,7为高斯滤波窗口大小，5 为均方差sigma的值，%0,04为推荐的k值，winsize为极大抑制窗口的大小且给参数时为奇数;%输出参数解释：posX为检测到角点X坐标，posY为检测到角点Y坐标，%cnt为检测到角点的个数，Out_Image为输出图像；posX,posY,cnt,Out_Image=

6、conerdetection(frame,7,2,0.04,7);%输出的图像已经是2值化的figure(2);imshow(Out_Image);holdon;plot(posX,posY,'ro','MarkerSize',15);disp(cnt);endfunctionposX,posY,cnt,Out_Image=conerdetection(frame,GaussWindow,sigma,k,winsize)%ImageData:gracyscaleimageofinput%GaussWindow:ThesizesofGausswindow%sigm

7、a:Thevariance%defaultvalue%winsize为极大抑制窗口大小Out_Image=frame;ImageData=frame;ImageData=double(ImageData(:,:,2);%ImageData数据%或者将 3 通矩阵式一个三通道的，我们角点标记只需选择一个可以;%ImageData=im2bw(ImageData,0.5);道的的图像转换成2值化的图像，完成提取%算法解释：%1：利用水平，竖直差分算子对图像的每个像素进行滤波以求得Ix,Iy,进而求得M矩阵中四个元素的值；%M=Ix*Ix,Ix*Iy;Ix*Iy,Iy*Iyorig_image=Im

8、ageData;fx=-2,-1,0,1,2;Ix=filter2(fx,orig_image);fy=-2;-1;0;1;2;Iy=filter2(fy,orig_image);Ix2=Ix.*Ix;Iy2=Iy.*Iy;Ixy=Ix.*Iy;%2:对M的四个元素进行高斯平滑滤波，得到新的矩阵M;%滤波平滑，消除突出点,得到新的矩阵M;h=fspecial('gaussian',GaussWindow,GaussWindow,sigma);%建立滤波算子Ix2=filter2(h,Ix2);%filter2是用h滤波器放在Ix2移动进行模板滤波Iy2=filter2(h,Iy

9、2);%消除y 方向上的突兀点Ixy=filter2(h,Ixy);%提取前的图像矩阵的预处理height=size(orig_image,1);高width=size(orig_image,2);给宽result=zeros(height,width);值为1R=zeros(height,width);相同的零矩阵Rmax=0;值%返回图像矩阵的行数给%返回图像矩阵的列数% 纪录角点位置,角点处%创建与图像矩阵大小% 图像中最大的 R%3:接下来利用M计算对应于每个像素的角点响应函数Cim(即R);%计算公式为：R=det(M)-k*(trace(M)八2,其中k为一个任意数，经验选取0.0

10、4就可以;%由于k的取值有些太随意，因此，改用此公式，用新的公式定义R：R=det(M)/Tr(M);%即Cim=R=Ix*Ix*Iy*Iy-(Ix*Iy)*(Ix*Iy)/Ix*Ix+Iy*Iy;fori=1:heightforj=1:widthM=Ix2(i,j),Ixy(i,j);Ixy(i,j),Iy2(i,j);%自相关矩阵% 计算 RR(i,j尸det(M)-0.04*(trace(M)八2;值,det()求一个方阵的行列式(Determinant);trace()求方阵的迹,即该方阵对角线上元素之和;ifR(i,j)>RmaxRmax=R(i,j);endendend%wi

11、nsize为非极大抑制窗口neighborhoodwinr=(winsize-1)/2;%the radius of theistart=winr+1;jstart=winr+1;iend=height-winr;jend=width-winr;cnt=0;fori=istart:iendforj=jstart:jendsubr=R(i-winr):(i+winr),(j-winr):(j+winr);%取出winr*winr这块区域里面的的矩阵;subrmax=max(max(subr);if(R(i,j)>k*Rmax)&&(R(i,j)=subrmax)result(i,j)=1;cnt=cnt+1;endendendposY,posX=find(result=1);%cnt为检测出来的角点的个数；%