Harris角点检测算法原理及其MATLAB编程实现

Harris角点检测算法原理及其MATLAB编程实现对于图像I(x,y),当在点(x,y)处平移(,)xy后的自相似性可以通过自相关函数给出：2(,)(,)(,,,)(,)((,)(,))uvWxycxyxywuvIuvIuxvy(1)其中，W(x,y)是以点(x,y)为中心的窗口，w(x,y)为加权函数，既可以用常数表示，又可以用高斯加权函数表示。根据泰勒展开，对图像I(x,y)在平移(,)xy后进行一阶近似：(,)(,)(,)(,)(,)(,),(,)xyxyIuxvyIuvuvxuvyxIuvuvuvyIIII(2)在公式(2)中，xI、yI是图像I(x,y)的偏导数。这样，公式(1)就能近似表示为：22(,,,)(,)((,)(,))(,)(,)wwcxyxyxxyMxyxyyIuvIuxvyxuvuvIIy(3)其中2222(,)(,)(,)(,)(,)(,)(,)(,)(,)(,)(,)(,)(,)xxywxyyxxywwxyywwuvuvMxyuvuvuvuvACCBuvuvuvIIIuvIIIuvIIIuvIII(4)也就是说图像I(x,y)在点(x,y)处平移(,)xy后的自相关函数可以近似为二次项函数：(,,,)(,)xcxyxyxyMxyy(5)根据二次项函数特征值的计算方法，可以求式(4)的特征值。Harris算法中给出的角点判别方法并不需要计算具体的特征值，而是计算一个角点响应值R来判断角点[1][2]。R的计算公式为2det()RMtraceM(6)式中，detM为矩阵(,)ACMxyCB的行列式；traceM为矩阵M的直迹；为经验常数，取值范围为0.04~0.06.事实上，特征值是隐含在detM和traceM中的，因为：21212detMABtraceMABC(7)Harris角点检测实现步骤如下：1)求出I(x,y)在x、y方向上的的梯度xI、yI;2)分别求出在x、y方向上的梯度乘积,2xxxIII2yyyIIIxyxyIII;3)对2xI、2yI和xyI进行高斯加权，从而产生A、B、C三个元素，如下：22()xxAgwII22()yyBgwII22()xyxyCgwII4)求出每个像素的Harris响应值R，令小于阀值的响应值R为零；2:det()RRMttraceM5)进行3×3邻域非极大值抑制，图像中的角点用局部极大值的点表示；6)记录下角点在原图像中的位置，即图像角点所在位置。MATLAB运行结果如下：附MATLAB代码：[filename,pathname,~]=uigetfile('*.jpg','选择JPG格式图片');if~ischar(filename)returnendstr=[pathnamefilename];pic=imread(str);iflength(size(pic))==3img=rgb2gray(pic);end[m,n]=size(img);tmp=zeros(m+2,n+2);tmp(2:m+1,2:n+1)=img;%扩展图像边缘1个像素Ix=zeros(m+2,n+2);Iy=zeros(m+2,n+2);Ix(:,2:n+1)=tmp(:,3:n+2)-tmp(:,1:n);%x方向上的差分Iy(2:m+1,:)=tmp(3:m+2,:)-tmp(1:m,:);%y方向上的差分Ix2=Ix(2:m+1,2:n+1).^2;Iy2=Iy(2:m+1,2:n+1).^2;Ixy=Ix(2:m+1,2:n+1).*Iy(2:m+1,2:n+1);h=fspecial('gaussian',[77],2);%生成高斯滤波器，消弱噪声的影响Ix2=filter2(h,Ix2);%滤波Iy2=filter2(h,Iy2);Ixy=filter2(h,Ixy);R=zeros(m,n);fori=1:mforj=1:nM=[Ix2(i,j)Ixy(i,j);Ixy(i,j)Iy2(i,j)];%构造Hessian矩阵R(i,j)=det(M)-0.06*(trace(M))^2;%角点判定公式endendRmax=max(max(R));loc=[];%记录角点位置tmp(2:m+1,2:n+1)=R;%扩展图像边缘1个像素fori=2:m+1forj=2:n+1iftmp(i,j)0.01*Rmax%要求矩阵R中的每个元素值大于0.01倍的最大值sq=tmp(i-1:i+1,j-1:j+1);sq=reshape(sq,1,9);sq=[sq(1:4),sq(6:9)];iftmp(i,j)sq%局部非极大值抑制loc=[loc;[j-1,i-1]];endendendend%以下代码显示获取的Harris角点X=loc(:,1);Y=loc(:,2);subplot(1,2,1);imshow(pic);subplot(1,2,2);imshow(pic);holdonplot(X,Y,'*');holdoff