CV——图像角点检测（Moravec Harris SUSAN）

一、介绍

1. CV中常常使用角点检测完成两张图像的匹配工作
2. 本人使用了Matlab实现了Moravec Harris SUSAN 三种角点检测算法

二、Moravec

原理

1. 设置一个小窗口，它不断在图像上移动
2. 公式含义为 小窗口上一个位置包囊像素之和 与 下一个位置包囊像素之和 的 平方差
3. 小窗口会遇到的三种情况，一种是平滑区域，一种是直线边缘区域，一种是角点
   我们可以想象，当小窗口移动到角点时的前后方差是最大的
4. 小窗口的移动方向设置为四个：右 、下 、左下 、右下
5. 设置一个阈值，当平方值为多大，才将其视为角点

步骤

1. 确定移动窗口的大小
2. 确定窗口中心移动范围
3. 计算窗口中心偏移到四个方向上时，像素之和引起的移动平方差，取最小值
4. 设置平方差阈值作为是否为角点的判断标准

代码

% 1、设置圆周模板 并 提取圆周模板的逻辑地址 并将其转换为 局部地址形式
radius=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;

% 2、读取图像 H  W 确定模板原心的水平和垂直运动范围
I = imread('lena.png');
[H,W,~]=size(I);
Xstep=1;Ystep=1;
nucleas_X=radius+1:Xstep:W-radius;
nucleas_Y=radius+1:Ystep:H-radius;

% 3、模板在图像上滑动
threshold=200;%——自定义阈值参数
result=zeros(H,W);
for y=nucleas_Y %行
    for x=nucleas_X %列
       delta=zeros(4,1);
       for e=1:length(tem_x)%不同中心的模板像素差值累加
           if y+1<length(nucleas_Y) && x+1<length(nucleas_X) && x-1>1 && y-1>1%出界判断
                delta(1)=delta(1)+power(I(y+1+tem_y(e),x+tem_x(e))-I(y+tem_y(e),x+tem_x(e)),2);%右
                delta(2)=delta(2)+power(I(y+tem_y(e),x+1+tem_x(e))-I(y+tem_y(e),x+tem_x(e)),2);%下
                delta(3)=delta(3)+power(I(y+1+tem_y(e),x+1+tem_x(e))-I(y+tem_y(e),x+tem_x(e)),2);%右下
                delta(4)=delta(4)+power(I(y-1+tem_y(e),x-1+tem_x(e))-I(y+tem_y(e),x+tem_x(e)),2);%左下
           end
       end
       o_min=min(delta(:));%获取四个方向上平方差的最小值
       if o_min>threshold
            result(y,x)=1;
       end
    end
end

% 4、显示边缘检测结果——白色为边缘
subplot(1,2,1),imshow(result,[]),title('Moravec角点检测结果');

% 5、角点的获取及绘制
[posr, posc] = find(result == 1);
subplot(1,2,2),imshow(I),hold on,plot(posc,posr,'r+'),title('角点在原图上绘制');

结果

三、Harris

原理

相较于Moravec有三点改进

1. 窗口移动方向不再是4个
2. 窗口数值大小不再是简单的0 1
3. 角点判断不再是 简单的一个平方差阈值

1. 对小窗口移动公式利用泰勒展开式进行化简计算——窗口移动方向数量增加
   

2. 化简得到实对称矩阵——包含图像的许多边缘信息
   

3. 对角化处理，即利用Ix，Iy进行坐标系的构建，散点图绘制后进行椭圆曲线拟合
   

4. 从绘制的椭圆曲线特征得到，长短半轴均比较大时，为角点
   

5. 将以上步骤绘制成下列公式
   

步骤

1. 计算自适应矩阵需要的参数 并对其进行高斯平滑
2. 构建自适应矩阵
3. 计算上述提到的公式 设置阈值 提取角点
4. 非极大抑制获取角点

代码

% 1、读取图像 H W
ori_im = imread('lena.png');    
[H,W,~]=size(ori_im);

% 2、获取自适应矩阵参数 并对其进行高斯平滑
fx = [-2 -1 0 1 2];   Ix = filter2(fx,ori_im);
fy = [-2;-1;0;1;2];   Iy = filter2(fy,ori_im);     
Ix2 = Ix.^2;    Iy2 = Iy.^2;   Ixy = Ix.*Iy;
h= fspecial('gaussian',[7 7],2);% 产生7*7的高斯窗函数，size为7*7，sigma为2
Ix2 = filter2(h,Ix2);
Iy2 = filter2(h,Iy2);	
Ixy = filter2(h,Ixy);	

% 3、计算自适应矩形 和 指定公式
R = zeros(H,W);
for i = 1:H
    for j = 1:W
        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

% 4、根据阈值挑选角点
BIN=zeros(H,W);
BIN(R>0)=1;
subplot(1,3,1),imshow(~BIN,[]),title('Harris角点检测结果');

% 5、非极大值抑制获取角点
result = zeros(H,W);
for i = 2:H-1
    for j = 2:W-1
        %当前角点为8领域的最大值
        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))])
            result(i,j) = 1;
        end
    end
end
subplot(1,3,2),imshow(result),title('非极大值抑制后');

% 6、角点的获取及绘制
[posr, posc] = find(result == 1);
subplot(1,3,3),imshow(ori_im),hold on,plot(posc,posr,'r+'),title('角点在原图上绘制');

结果

四、SUSAN

思路

1. 通过7*7圆形模板实现图像的滑动检测，相似灰度值的区域称为USAN；
   

2. 当模板遇到角点，USAN为整体面积的1/4；当模板遇到直线边缘，USAN为整体面积的1/2；当模板遇到平滑区域，USAN几乎等于整体面积；
   

3. 具体定义 其他像素与中心像素 差值为多少时，才算相似
   

4. 具体定义USAN区域占整体区域多少时，被定义为角点
   

步骤

1. 设置圆周模板 并 提取圆周模板的逻辑地址 并将其转换为 局部地址形式
2. 确定模板核的运动范围
3. 确定差值小于多少 周围像素与中心像素 才为相似
4. 用模板对图像滑动进行相似度累计
5. 确定多大的USAN区域将被视为角点
6. 非极大抑制检测角点

代码

% 1、设置圆周模板 并 提取圆周模板的逻辑地址 并将其转换为 局部地址形式
radius=3;
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;

% 2、读取图像 H  W 确定模板原心的水平和垂直运动范围
I = imread('lena.png');
[H,W,~]=size(I);
Xstep=1;Ystep=1;
nucleas_X=radius+1:Xstep:W-radius;
nucleas_Y=radius+1:Ystep:H-radius;

% 3、使用模板扫描图像确定 USAN大小
t=15;%差值为多少被判断为相似
USAN=zeros(H,W);
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  
    end
end

% 4、角点挑选 小于USAN最大值的多少视为角点
g=2/3*max(USAN(:));%——自定义阈值
R=zeros(H,W);
for i=1:H
    for j=1:W
        if USAN(i,j)<g
            R(i,j)=1;
        else
            R(i,j)=0;
        end
    end
end
subplot(1,3,1),imshow(R),title('SUSAN角点检测结果');

% 5、非极大值抑制获取角点
corners=zeros(H,W);
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))]) %当前角点为8领域的最大值
            corners(i,j)=1;
        end
   end
end
subplot(1,3,2),imshow(corners),title('非极大值抑制后');

% 6、角点的获取及绘制
[posr, posc] = find(corners == 1);
subplot(1,3,3),imshow(I),hold on,plot(posc,posr,'r+'),title('角点在原图上绘制');

结果