moravec.m

利用moravec兴趣算子提取图像的点特征

clear all
close all
clc

tic
%load lvdongdata
% imshow(I)

I=imread('girl.bmp');
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('图像的行数')