完成图像的角点检测实验

Harris 角点算法实现

根据上述讨论，可以将 Harris 图像角点检测算法归纳如下，共分以下五步：

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

Harris 角点的性质

参数α对角点检测的影响假设已经得到了矩阵 M 的特征值λ1≥λ2≥0，令λ2=kλ1,0≤k≤1。由特征值与矩阵 M 的直迹和行列式的关系可得： detM=∏iλi traceM=∑iλi 从而可以得到角点的响应 R=λ2λ2=α(λ2+λ2)2=λ2(k?α(1+k)2) 假设 R≥0，则有： 0≤α≤k(1+k)2≤0.25 对于较小的 k 值，R≈λ2(k?α),α<k。由此，可以得出这样的结论：增大α的值，将减小角点响应值 R，降低角点检测的灵性，减少被检测角点的数量；减小α值，将增大角点响应值 R，增加角点检测的灵敏性，增加被检测角点的数量。
Harris 角点检测算子对亮度和对比度的变化不敏感这是因为在进行 Harris 角点检测时，使用了微分算子对图像进行微分运算，而微分运算对图像密度的拉升或收缩和对亮度的抬高或下降不敏感。换言之，对亮度和对比度的仿射变换并不改变 Harris 响应的极值点出现的位置，但是，由于阈值的选择，可能会影响角点检测的数量。
Harris 角点检测算子具有旋转不变性 Harris 角点检测算子使用的是角点附近的区域灰度二阶矩矩阵。而二阶矩矩阵可以表示成一个椭圆，椭圆的长短轴正是二阶矩矩阵特征值平方根的倒数。当特征椭圆转动时，特征值并不发生变化，所以判断角点响应值 R 也不发生变化，由此说明 Harris 角点检测算子具有旋转不变性。
Harris 角点检测算子不具有尺度不变性如下图所示，当右图被缩小时，在检测窗口尺寸不变的前提下，在窗口内所包含图像的内容是完全不同的。左侧的图像可能被检测为边缘或曲线，而右侧的图像则可能被检测为一个角点。

close all;
clear all;
clc;
 
img=imread('testc.jpg');
imshow(img);
img = im2gray(img);
img =double(img);
[m n]=size(img);
 
tmp=zeros(m+2,n+2);
tmp(2:m+1,2:n+1)=img;
Ix=zeros(m+2,n+2);
Iy=zeros(m+2,n+2);
 
E=zeros(m+2,n+2);
 
Ix(:,2:n)=tmp(:,3:n+1)-tmp(:,1:n-1);
Iy(2:m,:)=tmp(3:m+1,:)-tmp(1:m-1,:);
 
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',[7 7],2);
Ix2=filter2(h,Ix2);
Iy2=filter2(h,Iy2);
Ixy=filter2(h,Ixy);
 
Rmax=0;
R=zeros(m,n);
for i=1:m
    for j=1:n
        M=[Ix2(i,j) Ixy(i,j);Ixy(i,j) Iy2(i,j)];
        R(i,j)=det(M)-0.06*(trace(M))^2;
        
        if R(i,j)>Rmax
            Rmax=R(i,j);
        end
    end
end
re=zeros(m+2,n+2);
 
tmp(2:m+1,2:n+1)=R;
img_re=zeros(m+2,n+2);
img_re(2:m+1,2:n+1)=img;
for i=2:m+1
    for j=2:n+1
        
        if tmp(i,j)>0.02*Rmax &&...
           tmp(i,j)>tmp(i-1,j-1) && tmp(i,j)>tmp(i-1,j) && tmp(i,j)>tmp(i-1,j+1) &&...
           tmp(i,j)>tmp(i,j-1) && tmp(i,j)>tmp(i,j+1) &&...
           tmp(i,j)>tmp(i+1,j-1) && tmp(i,j)>tmp(i+1,j) && tmp(i,j)>tmp(i+1,j+1)
                img_re(i,j)=255; 
        end   
    end
end
 
img_re=mat2gray(img_re(2:m+1,2:n+1));
figure,imshow(img_re);
复制代码

function re=get_coords(angle)       %angle是边缘法线角度，返回法线前后两点
    sigma=0.000000001;
    x1=ceil(cos(angle+pi/8)*sqrt(2)-0.5-sigma);
    y1=ceil(-sin(angle-pi/8)*sqrt(2)-0.5-sigma);
    x2=ceil(cos(angle-pi/8)*sqrt(2)-0.5-sigma);
    y2=ceil(-sin(angle-pi/8)*sqrt(2)-0.5-sigma);
    re=[x1 y1 x2 y2];

end
复制代码