计算机视觉与模式识别实验2-1 角点检测算法（Harris，SUSAN，Moravec）

🧡🧡实验流程🧡🧡

Harris算法

Harris算法实现步骤：

• w表示窗函数，在这个算法中，用于对梯度图像的平方进行滤波，以减少图像中的噪声，从而提高角点检测的准确性。
• R(i,j) 表示图像中像素点 (i,j) 处的角点响应值。这个值描述了在该像素点周围区域的灰度变化情况，从而表征了该像素是否可能是角点。通常情况下， R(i,j) 的值越大，表示该像素点越可能是角点。因此对于每个像素点，非极大值抑制会检查其周围邻域内的角点响应值，当它是该邻域内的最大值时，才将该像素点的角点响应值保留下来。

SUSAN算法

SUSAN角点检测算法实现步骤：
①　定义局部邻域：对于图像中的每个像素，定义一个局部邻域，通常是以该像素为中心的圆形区域。
②　计算灰度差异度：对于局部邻域内的每个像素，计算它与中心像素的灰度差异度。通常使用的是像素间的灰度差的绝对值作为度量。
③　确定相似像素：根据计算得到的灰度差异度，判断每个像素是否属于同一类别。具体而言，如果像素与中心像素的灰度差异度小于一定阈值，则认为它们属于同一类别；否则，它们属于不同的类别。
④　统计局部邻域中相似像素的数量：统计局部邻域中与中心像素属于同一类别的像素数量。
⑤　角点判断：如果局部邻域中相似像素的数量小于一定阈值，则认为该中心像素是一个角点

原始图像

Harris、SUSAN角点检测效果如下：

分析：直观上看，SUSAN标记的角点数更多，由于它是基于像素灰度值的局部统计特征来判断角点的，因此在一些局部特征明显的图像部位中效果较好。它们在理论实现方面有一些相似之处（都需要非极大值抑制），不同的是Harris角点检测算法基于图像局部区域的灰度变化来检测角点，通过计算图像的梯度和自相关矩阵的特征值来判断是否存在角点；而SUSAN 角点检测算法基于像素点的灰度值与周围像素的灰度值的差异来检测角点，通过比较像素与其周围邻域的灰度值之间的差异来判断是否为角点。另外从计算效率感受，SUSAN角点检测算法的计算相对较简单和高效。

Moravec算法

实现步骤：
①　计算窗口中心偏移到四个方向上时的灰度平方差，取最小值，作为兴趣值IV。
②　设置阈值T，将兴趣值大于阈值的点作为候选点（IV>T）。
③　局部区域选取候选点中的极值点作为特征点。 (局部非最大值抑制)
其中，四个方向灰度差的计算举例如下：

效果如下：

分析：可以看到，在参数设置合理的情况下，检测效果不比上述两种算法差，并且检测出了上述方法没有测出的角点，如嘴角处、肩膀处等等。另外，由于算法的局限性，在灰度值很相近的情况下，可能产生误检测，比如这里头发处其实较多误检测。

综上所述，从实验直观感受来说，三种算法总结：

• 精度上：Harris 算法更优。
• 鲁棒性：SUSAN 对于噪声有较好的鲁棒性。
• 计算效率和理论理解上：Moravec 实现简单，对精度要求不是特别高，可以选择 Moravec 算法。

🧡🧡全部代码🧡🧡

%%% Harris


clear; clc ;tic;
ori_im = imread('img/test1_Lena.tif');     % 读取图像
if(size(ori_im,3)==3)
    ori_im = rgb2gray(uint8(ori_im));  %转为灰度图像
end
imshow(ori_im);
% ===================计算图像的梯度===================
% fx = [5 0 -5;8 0 -8;5 0 -5];          % 高斯函数一阶微分，x方向(用于改进的Harris角点提取算法)
fx = [-2 -1 0 1 2];                 % x方向梯度算子(用于Harris角点提取算法)
Ix = filter2(fx,ori_im);              % x方向滤波
% fy = [5 8 5;0 0 0;-5 -8 -5];          % 高斯函数一阶微分，y方向(用于改进的Harris角点提取算法)
fy = [-2;-1;0;1;2];                 % y方向梯度算子(用于Harris角点提取算法)
Iy = filter2(fy,ori_im);              % y方向滤波

% ===================计算自相关矩阵===================
Ix2 = Ix.^2;
Iy2 = Iy.^2;
Ixy = Ix.*Iy;
clear Ix;
clear Iy;
h= fspecial('gaussian',[7 7],2);      % 产生7*7的高斯窗函数，sigma=2
Ix2 = filter2(h,Ix2);
Iy2 = filter2(h,Iy2);
Ixy = filter2(h,Ixy);

% ===================计算角点响应值 R===================
height = size(ori_im,1);
width = size(ori_im,2);
result = zeros(height,width);  % 纪录角点位置，角点处值为1
R = zeros(height,width);
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)-0.06*(trace(M))^2;   
    end
end

% ===================进行非极大值抑制===================
cnt = 0;
for i = 2:height-1
    for j = 2:width-1
        % 进行非极大抑制，窗口大小3*3
        if  R(i,j) > R(i-1,j-1) && R(i,j) > R(i-1,j) && R(i,j) > R(i-1,j+1) && R(i,j) > R(i,j-1) && R(i,j) > R(i,j+1) && R(i,j) > R(i+1,j-1) && R(i,j) > R(i+1,j) && R(i,j) > R(i+1,j+1)
            result(i,j) = 1;
            cnt = cnt+1;
        end
    end
end

% ===================对检测到的角点按照角点响应值进行排序，然后取前 ps 个角点（设定阈值）===================
Rsort=zeros(cnt,1);
[posr, posc] = find(result == 1);
for i=1:cnt
    Rsort(i)=R(posr(i),posc(i));
end
[Rsort,ix]=sort(Rsort,1);
Rsort=flipud(Rsort);
ix=flipud(ix);
ps=100;
posr2=zeros(ps,1);
posc2=zeros(ps,1);
for i=1:ps
    posr2(i)=posr(ix(i));
    posc2(i)=posc(ix(i));
end

% ===================显示原始图像，并在图像上标记检测到的角点===================
imshow(ori_im);
hold on;
plot(posc2,posr2,'g+');
title('Harris');

%%% SUSAN


%设置圆周模板半径和滑动窗口的步长
radius=3;Xstep=1;Ystep=1;
template=fspecial('disk',radius);
template(template>0.01)=1;  %模板二值化
template(template<0.01)=0;
 %提取圆周模板的逻辑地址
[tem_x,tem_y]=find(template==1);
tem_x=tem_x-radius-1;
tem_y=tem_y-radius-1;

t=45; %USAN判定阈值
I = imread('img/test1_Lena.tif');
W=size(I,2);H=size(I,1); %图像大小
nucleas_X=radius+1:Xstep:W-radius; %模板圆心即nucleas运动范围
nucleas_Y=radius+1:Ystep:H-radius;
USAN=zeros(size(I,1),size(I,2));  %初始化USAN累加器
%网格遍历
tic;
for y=nucleas_Y
    for x=nucleas_X
       for e=1:length(tem_x)  %圆周模板上进行判定
           delta=I(y+tem_y(e),x+tem_x(e))-I(y,x);
           if delta<t
               USAN(y,x)=USAN(y,x)+1;  %低于阈值则收纳
           end
       end  
        fprintf(strcat('已处理第','(',num2str(y),',',num2str(x),')','像素点..\n'));
    end
end
toc;
%边缘检测
%
%formula： R(r0) =   g-USAN(r0)  , if USAN(r0)<g
%                    0           ,  otherwise   
%          where g=3/4*max(USAN(:)), USAN越小边缘响应越强,其中角点极小，对应R极大
g=1/2*max(USAN(:));
R=zeros(H,W);
for i=1:size(USAN,1)
    for j=1:size(USAN,2)
        if USAN(i,j)<g
            R(i,j)=g-USAN(i,j);
        else
            R(i,j)=0;
        end
    end
end
BIN=zeros(H,W);
BIN(R>0)=1;
figure(1),imshow(BIN,[]),title('SUSAN  边缘检测');

%角点检测（非极大值抑制）
corners=[];
for i=2:H-1
   for j=2:W-1 
        if R(i,j)>max([max(R(i-1,j-1:j+1)),R(i,j-1) ,R(i,j+1), max(R(i+1,j-1:j+1))])
            corners=[corners;[i,j]];
        end
   end
end
% I = insertMarker(I,corners);
figure(2);imshow(I);hold on;
set(gca,'xaxislocation','top','yaxislocation','left','ydir','reverse');
scatter(corners(:,2),corners(:,1),'x','g'),title('SUSAN');
hold off;

%%% Moravec

close all;
clear all;
clc;
img=(imread('img/test1_LenaRGB.tif'));
img = double(rgb2gray(img));
[X,Y]=size(img);
img2 = imread ('img/test1_LenaRGB.tif');
img2 = rgb2gray(img2);

imgn=zeros(X,Y);
n=2;%控制窗口大小。窗口大小为 2*n+1

% ====计算四个方向的灰度变化平方====
for x=1+n:X-n
   for y=1+n:Y-n
       window=img(x-n:x+n,y-n:y+n);%窗口为5*5
       V=zeros(1,4);
       for i=2:2*n+1        %垂直，水平，对角，反对角四个方向领域灰度差的平方和
            V(1)=V(1)+(window(i,n+1)-window(i-1,n+1))^2;%垂直
            V(2)=V(2)+(window(n+1,i)-window(n+1,i-1))^2;%水平
            V(3)=V(3)+(window(i,i)-window(i-1,i-1))^2;%正对角线
            V(4)=V(4)+(window(i,(2*n+1)-(i-1))-window(i-1,(2*n+1)-(i-2)))^2;%反对角线
       end

       IV=min(V);          %四个方向中选最小值
       imgn(x,y)=IV;      
   end
end

% =====设定阈值T，>=T为兴趣点（潜在角点）====
T=2000; % 阈值
ind=find(imgn<T);
imgn(ind)=0;

imshow(img2);hold on;title('Moravec算子提取');
m=3;
count = 0;
% ====非极大值抑制====
for x=1+m:X-m           %选局部最大（窗口为11*11）且非零值作为特征点
    for y=1+m:Y-m
        window2=imgn(x-m:x+m,y-m:y+m);
        if max(window2(:))==imgn(x,y) && imgn(x,y)~=0
            plot(y,x,'+','color','g');
            count=count+1; 
        end
    end
end