角点检测（1）Moravec's 算子

角点检测 流程

1.      使用角点操作子 : 这个步骤将图片和角点操作子需要的少数特定参数作为输入.对于输入图片的每个像素点将使用角点操作子为其生成一个角点幅值.角点幅值是一个简单的数字用于指明角点操作子将一个像素判断为角点的程度.兴趣点角点检测算法在角点操作子如何生成幅值上各有不同,但所有的算法在计算幅值时都只考虑一个以像素为中心点的小窗口中的像素.这个步骤将输出角点幅值表.因为对于输入图像的每个像素都会使用角点操作子来生成其角点幅值,因此角点幅值表是否和输入图片有同样的尺寸也可以作为输入图片的一种处理方式.

 

2.      角点幅值表门限 : 兴趣点角点检测子将角点定义为角点幅值表中的局部极大值点.从这点来看,角点幅值表中有很多局部极大值点但其中的一部分角点幅值相对较低,不是真正的角点.为了避免将这些点判断为角点,需要一个确定的门限值.角点幅值表中所有低于这个门限的的值都设为0.这个门限值的选择依赖于场景,以及测试和错误实验.门限值应设得足够高以移除不是真正角点的局部极大值,但也要足够低以保留描述真正角点的局部极大值.实际上只有极少的门限值可以移除假的角点并保留真的角点,所以要以应用场景的需求为基础来对门限值进行权衡平衡.

 

3.      非极大值抑制 : 经过门限值处理的角点幅值表中只有周围均为非零值的局部极大值点才考虑标记为角点.为了定位这样的局部极值点,需要使用非极大值抑制.对于角点幅值表中的每个点,如果一点的角点幅值小于以它为中心一定范围内其它的点,则将该点的幅值设为零.经过非极大值抑制,可以简单的从角点幅值表中剩下的非零点来确定角点.

Moravec's 算子算法

Moravec角点检测子的操作过程如下所示

图像中像素(x,y)的值由I(x,y)指出

输入:灰度图片,窗口尺寸,门限值T

输出:角点坐标表

1.      对于图片中的每个像素(x,y)移动(u,v)后计算其变化幅值,(u,v)使用如下参数(1,0),(1,1),(0,1),(-1,1),(-1,0),(-1,-1),(0,-1),(1,-1)

2.      对于每个像素(x,y)计算其角点幅值C(x,y)以建立角点幅值表:C(x,y)-min(V_(u,v)[x,y])

3.      对角点幅值表进行门限处理,将所有低于门限T的C(x,y)设为0

4.      使用非极大值抑制操作找出局部极大值点,角点幅值表中所有非0点即为角点.

限制

算子的旋转响应

Moravec算子是方向敏感的,因为变化幅值只在几个分离的方向集合上进行计算(确切的说,在8个主方向),Moravec算子不具有旋转不变性.这导致了角点检测子的可重复率很弱.

噪声响应

Moravec算子使用的窗口为正方形并且是二元整数.为了得到更精准的局部变化幅值可以使用圆形窗口,因为圆形窗口中心像素到窗口边沿的欧几里得距离在所有方向都相等.即使在不连续的坐标格子上,使用圆形窗口也能使局部变化幅值的精度提高.

边沿上的过激反应

边沿定义为灰阶图像中梯度变化较大的位置.也就是说,边沿是在某个方向上变化十分迅速而剧烈的位置.Moravec算子对于边沿的响应十分敏感,因为边沿上的任何瑕疵都会马上识别为噪声,导致生成的局部最小变化幅值变得异常大.如测试图片结果小结所展示的,边沿在确定的方向上会过激响应导致其变化幅值在8个主方向上进行计算.

结论

Moravec角点检测子是Moravec使用的一种非常简单的算法,现在基本可以认为是废弃的.它不具有旋转不变性(大多数现代角点检测子具有的特性),尽管它的响应各向异性,噪声敏感,易于在边沿上和独立像素点报告假角点.但是它的高效计算性是Moravec所需要的,因为Moravec的兴趣是在他的应用上拥有足够的实时性能和最少的计算消耗.过去数十年中随着计算能力的不断增强,现在大多数应用会选择计算量大但性能更好的角点检测子.这些角点检测子中使用极为广泛之一就是现在的Harris/Plessey角点检测子.它的计算十分费力但直接解决了Moravec角点检测子中的很多局限.

clear all 
close all 
clc 
 
tic  
I=imread('girl_.jpg'); 
originalmap=I; 
%计算各像元的兴趣值IV(Interest value) 
 
%以w*w为像元对图像I进行分割 
w=5; 
%分割为m*n个像元 
m=floor(size(I,1)/w); %行数，round 
n=floor(size(I,2)/w); %列数 
%得到各中心点的坐标值 
c=ceil(w/2)+(0:n-1)*w; %所有中心点x坐标的值 
r=ceil(w/2)+(0:m-1)*w; %所有中心点y坐标的值 
 
%计算各个像元的兴趣值 
step=floor(w/2); 
for y_unit=1:m 
    for x_unit=1:n 
        %计算横向相邻像素灰度差的平方和 
        v1=0;v2=0;v3=0;v4=0; 
        for i=-step:step-1 
            %计算横向相邻像素灰度差的平方和 
            v1=v1+(I(r(y_unit),c(x_unit)+i)-I(r(y_unit),c(x_unit)+i+1))^2; 
            %计算右斜线方向相邻像素灰度差的平方和 
            v2=v2+(I(r(y_unit)+i,c(x_unit)+i)-I(r(y_unit)+i+1,c(x_unit)+i+1))^2; 
            %计算纵向相邻像素灰度差的平方和 
            v3=v3+(I(r(y_unit)+i,c(x_unit))-I(r(y_unit)+i+1,c(x_unit)))^2; 
            %计算左斜线方向相邻像素灰度差的平方和 
            v4=v4+(I(r(y_unit)-i,c(x_unit)+i)-I(r(y_unit)-i-1,c(x_unit)+i+1))^2; 
        end 
        %取v1,v2,v3,v4中最小者作为像素(c,r)的兴趣值 
        IV_cr(y_unit,x_unit)=min([v1,v2,v3,v4]); 
    end 
end 
 
%给定一个经验阈值。阈值的选取应以候选点包含所需要的特征点，而又不包含过多的非特征点为原则. 
exper_thr=20; 
IV_cr(IV_cr<exper_thr)=NaN; %将候选点中小于阈值的点全部去除 
 
%选择候选点的极值点为特征点 
 
%选择计算窗口大小 
wf=9; 
%对候选点进行分割,分割为mf*nf个区域 
mf=floor(m/wf); %行数 
nf=floor(n/wf); %列数 
 
%得到特征点的坐标 
xc=[]; 
yc=[]; 
for y_unit=1:mf 
    for x_unit=1:nf 
        %计算分割区域中的最大值 
        [C,I]=max(IV_cr((y_unit-1)*wf+1:y_unit*wf,(x_unit-1)*wf+1:x_unit*wf));%得到行 
        [C1,I1]=max(C);%得到列 
        IV_crch(y_unit,x_unit)=C1; %得到分割区域中的最大值 
        row=I(I1); 
        col=I1; 
        crch_row(y_unit,x_unit)=(y_unit-1)*wf+row; %得到最大值在候选区域中的行数 
        crch_col(y_unit,x_unit)=(x_unit-1)*wf+col; %得到最大值在候选区域中的列数 
        yc=[yc,r((y_unit-1)*wf+row)];xc=[xc,c((x_unit-1)*wf+col)]; 
        IV_cr((y_unit-1)*wf+1:y_unit*wf,(x_unit-1)*wf+1:x_unit*wf)=NaN; %先去掉所有的点 
        IV_cr((y_unit-1)*wf+row,(x_unit-1)*wf+col)=C1; %加上符合要求的候选点 
    end 
end 
 
%先给出原来的图像 
figure(1) 
imshow(originalmap) 
title('原图像') 
xlabel('横向') 
ylabel('纵向') 
axis on 
 
%描绘特征点分布 
figure(2) 
plot(xc,yc,'*') 
view(0,-90) 
title('特征点的分布') 
xlabel('图像的列数') 
ylabel('图像的行数') 
toc; 
t=toc; 
disp(['本程序的运行时间为',num2str(t),'秒。']) 
 
%图像中显示特征点 
figure(3) 
imshow(originalmap) 
hold on 
plot(xc,yc,'R*') 
axis on 
title('图像中显示特征点') 
xlabel('图像的列数') 
ylabel('图像的行数')