【算法】Harris角点算法matlab程序

Harris角点算法实现

基本原理

人眼对角点的识别通常是在一个局部的小区域或小窗口完成的。如果在各个方向上移动这个特征的小窗口，窗口内区域的灰度发生了较大的变化，那么就认为在窗口内遇到了角点。如果这个特定的窗口在图像各个方向上移动时，窗口内图像的灰度没有发生变化，那么窗口内就不存在角点；如果窗口在某一个方向移动时，窗口内图像的灰度发生了较大的变化，而在另一些方向上没有发生变化，那么，窗口内的图像可能就是一条直线的线段。

Harris角点算法实现

1. 计算图像I(x,y)在X和Y两个方向的梯度Ix、Iy
2. 计算图像两个方向梯度的乘积
3. 使用高斯函数对Ix^2、Iy^2和Ixy进行高斯加权（取σ=1），生成矩阵M的元素A、B和C
4. 计算每个像素的Harris响应值R，并对小于某一阈值t的R置为零
5. 在3×3或5×5的邻域内进行非最大值抑制，局部最大值点即为图像中的角点

Harris角点的性质

参数α对角点检测的影响:增大α的值，将减小角点响应值R，降低角点检测的灵性，减少被检测角点的数量；减小α值，将增大角点响应值R，增加角点检测的灵敏性，增加被检测角点的数量。
Harris角点检测算子对亮度和对比度的变化不敏感:这是因为在进行Harris角点检测时，使用了微分算子对图像进行微分运算，而微分运算对图像密度的拉升或收缩和对亮度的抬高或下降不敏感。换言之，对亮度和对比度的仿射变换并不改变Harris响应的极值点出现的位置，但是，由于阈值的选择，可能会影响角点检测的数量。
**Harris角点检测算子具有旋转不变性**Harris角点检测算子使用的是角点附近的区域灰度二阶矩矩阵。而二阶矩矩阵可以表示成一个椭圆，椭圆的长短轴正是二阶矩矩阵特征值平方根的倒数。当特征椭圆转动时，特征值并不发生变化，所以判断角点响应值R也不发生变化，由此说明Harris角点检测算子具有旋转不变性。
Harris角点检测算子不具有尺度不变性:如下图所示，当右图被缩小时，在检测窗口尺寸不变的前提下，在窗口内所包含图像的内容是完全不同的。左侧的图像可能被检测为边缘或曲线，而右侧的图像则可能被检测为一个角点。
Matlab代码实现(也可以通过调用matlab函数cornermetric函数实现):

function [posr,posc]=Harris1(in_image,a)% 功能：检测图像harris角点% in_image-待检测的rgb图像数组% a--角点参数响应，取值范围：0.04~0.06% [posr，posc]-角点坐标in_image=rgb2gray(in_image);I=double(in_image);%%%%计算xy方向梯度%%%%%fx=[-1,0,1];%x方向梯度模板Ix=filter2(fx,I);%x方向滤波fy=[-1;0;1];%y方向梯度模板(注意是分号)Iy=filter2(fy,I);%%%%计算两个方向梯度的乘积%%%%%Ix2=Ix.^2;Iy2=Iy.^2;Ixy=Ix.*Iy;%%%%使用高斯加权函数对梯度乘积进行加权%%%%%产生一个7*7的高斯窗函数，sigma值为2h=fspecial('gaussian',[7,7],2);IX2=filter2(h,Ix2);IY2=filter2(h,Iy2);IXY=filter2(h,Ixy);%%%%%计算每个像元的Harris响应值%%%%%[height,width]=size(I);R=zeros(height,width);%像素(i,j)处的Harris响应值for i=1:height    for j=1:width        M=[IX2(i,j) IXY(i,j);IXY(i,j) IY2(i,j)];        R(i,j)=det(M)-a*(trace(M))^2;    endend%%%%%去掉小阈值的Harris值%%%%%Rmax=max(max(R));%阈值t=0.01*Rmax;for i=1:height    for j=1:width        if R(i,j)<t            R(i,j)=0;        end    endend%%%%%进行3*3领域非极大值抑制%%%%%%%%%corner_peaks=imregionalmax(R);%imregionalmax对二维图片，采用8领域（默认，也可指定）查找极值，三维图片采用26领域%极值置为1，其余置为0num=sum(sum(corner_peaks));%%%%%%显示所提取的Harris角点%%%%[posr,posc]=find(corner_peaks==1);figureimshow(in_image);hold onfor i=1:length(posr)    plot(posc(i),posr(i),'r+');endend