摄影测量学点特征提取模拟实践.doc

辽宁工程技术大学
模 拟 实 践

教学单位 测绘学院
专 业 测绘工程
名 称 点特征提取模拟实践
班 级 测绘13-3
组 长 秦小茜
组 员 柳少冬 倪贺星 刘建琪
指导教师 徐辛超

一、目的
理解影像中每个像素灰度值的概念；理解点特征在灰度方面的特点；掌握常用的点特征提取算子，及其对应的点特征提取方法，包括：Moravec算子、Forstner算子。
二、要求
运用VC++、C#、MATLAB其中一种自己擅长的语言，实现两种点特征提取算法。
三、方法与步骤
(一) Moravec算子提取过程
a. 读取灰度影像，得到影像的灰度矩阵；
b. 计算各像元的兴趣值 IV：

c. 给定一经验阈值，将兴趣值()大于阈值的点作为候选点；
d. 选取候选点中的极值点作为特征点；
e. 将特征点在原始图像上进行标记、展示。
(二) Forstner算子提取过程
a. 计算各像素的Robert’s梯度 ；

b. 计算l´l（如5´5或更大）窗口中灰度的协方差矩阵 ；

其中：
c. 计算兴趣值q与w ；

式中：DetN代表矩阵N之行列式，trN代表矩阵N之迹。
d. 确定待选点 ；

当同时，该像元为待选点。
e. 选取极值点。
即在一个适当窗口中选择最大的待选点。
f. 将特征点在原始图像上进行标记、展示。

（三） 、流程图

四、 实验内容及成果
1、Moravec算子程序相关内容
（1）Moravec算子程序源代码：
clc
clear
% I=imread('C:\Users\LBS\Desktop\Moravec MATLAB\girl.bmp');
I=imread('C:\Users\Administrator\Desktop\点特征提取\Moravec MATLAB\left缩小.bmp');
originalmap=I;
% if isrgb(I)
% I=rgb2gray(I); %将彩色图像转换为灰度图像
% end
I=double(I); %转换为双精度图像
[row,column]=size(I);
w1 = floor(5/2);%确定窗口大小
w2 = 5;
step=w1;
ir=10000;%设置阈值
num=0;
%得到特征点的坐标
xc=[];

yc=[];
tic;
for i=w1+1:1:row-w1-1
if i+w2>row-w1-1
continue;
end
for j=w1+1:1:column-w1-1
if j+w2>column-w1-1
continue;
end
interest=0; %兴趣值
v1=0;v2=0;v3=0;v4=0;
for k=-step:step-1
%计算纵向相邻像素灰度差的平方和
v1=v1+(I(i+k,j)-I(i+k+1,j))^2;
%计算125度斜线方向相邻像素灰度差的平方和
v2=v2+(I(i+k,j+k)-I(i+k+1,j+k+1))^2;
%计算横相邻像素灰度差的平方和
v3=v3+(I(i,j+k)-I(i,j+k+1))^2;
%计算45度斜线方向相邻像素灰度差的平方和
v4=v4+(I(i+k,j-k)-I(i+k+1,j-k-1))^2;
end
%取v1,v2,v3,v4中最小者作为像素(c,r)的兴趣值
interest=min([v1,v2,v3,v4]);
if interest>ir %//最大兴趣值与阈值比较

num=num+1;
pointx(num)=j;
pointy(num)=i;
yc=[yc,i];
xc=[xc,j];
end
end
end

%图像中显示特征点
figure(4);
imshow(originalmap);
hold on;
plot(xc,yc,'R*');
axis on;
title('图像中显示特征点');
xlabel('图像的列数');
ylabel('图像的行数');

toc;
t=toc;
disp(['本程序的运行时间为',num2str(t),'秒。']);

图一Moravec算子程序的相关截图
（2） 原始灰度影像

图二原始灰度影像

（3）程序运行结果

图三Moravec算子程序运行结果
2、 Forstner算子程序的相关内容
（1） Forstner算子程序源代码：
clc
clear
I=imread('C:\Users\Administrator\Desktop\点特征提取\Moravec MATLAB\left缩小.bmp');
% I=imread('C:\Users\LBS\Desktop\Lenna.bmp');
originalmap=I;
% if isrgb(I)

% I=rgb2gray(I); %将彩色图像转换为灰度图像
% end
I=double(I); %转换为双精度图像
[row,column]=size(I);
w1 = 5;%确定窗口大小
w2 = 5;
k=floor(w1/2);
step=w1;
iq=0.65;%设置阈值
ir_w=2.5;
num=0;
%得到特征点的坐标
xc=[];
yc=[];
tic;
sum_w=0;
count_w=0;
for i=k+1:row-k
for j=k+1:column-k
%协方差矩阵
G0=0;G2=0;G1=0;
for m=i-k:i+k-1
for n=j-k:j+k-1
G0=G0+(I(m+1,n+1)-I(m,n))^2;
G2=G2+(I(m,n+1)-I(m+1,n))^2;

G1=G1+(I(m+1,n+1)-I(m,n))*(I(m,n+1)-I(m+1,n));
end
end
DetN=G0*G2-G1*G1;
trN=G0+G2;
if trN==0
ww(i,j)=0;
q(i,j)=0;
else
ww(i,j)=DetN/trN;
q(i,j)=DetN*4/trN/trN;
end
sum_w=sum_w+ww(i,j);
count_w=count_w+1;
end
end
avg_w=sum_w/count_w;
for i=k:w2:row-w1
if i+w2>row-w1
continue;
end
for j=k:w2:column-w1
if j+w2>column-w1
continue;
end

max_w=0
max_x=0;max_y=0;
for mh=1:w2
for mw=1:w2
if (max_w<ww(i+mh,j+mw))
max_w=ww(i+mh,j+mw);
max_x=j+mw;
max_y=i+mh;
end
end
end
if (max_w>ir_w*avg_w&&q(i+mh,j+mw)>iq) %//最大兴趣值与阈值比较
num = num+1;
pointx(num)=max_x;
pointy(num)=max_y;
yc=[yc,max_y];
xc=[xc,max_x];
end
end
end
figure;
imshow(originalmap);
hold on;
plot(xc,yc,'R*');

axis on;
title('图像中显示特征点');
xlabel('图像的列数');
ylabel('图像的行数');

图四Forstner算子程序相关截图
（2） 原始灰度影像

图五Forstner原始灰度影像
（3） 程序运行结果

图六 Forstner程序运行结果
五、实验体会
柳少冬：此次实验让我对Forstner算子以及Moravec算子实现过程有了进一步了解，对这两种点特征提取算法有了更深刻的认识。Forstner和Moravec程序的大致流程看了看，但是还有些许问题。这两个程序需要修改的不多，只改了改两种算子所用图片的路径以及部分内容，总体操作比较顺利，结果也容易得到。这次实验通过程序运行，得到经两种点特征提取的图片，从照片中可以清晰对比不同。总之，此次实验比较顺利。