综合训练三.docx
《综合训练三.docx》由会员分享，可在线阅读，更多相关《综合训练三.docx（17页珍藏版）》请在冰豆网上搜索。
综合训练三
辽宁工程技术大学
本科生综合训练报告书
教学单位测绘学院
专业测绘工程
名称特征提取
班级测绘14-3班
学生姓名曹辰
学号1404070302
指导教师徐辛超
一、目的要求
目的和要求：
理解遥感图像特征提取的方法及基本原理，达到培养学生独立完成遥感图像特征提取的能力；要求学生借助编程软件，实现点、线特征的提取。
通过熟悉提取过程进一步理解特征提取领会实现过程与基本原理间的对应关系。
二、仪器设备
实习时间：
6月15日
实验地点：
博文楼116
刘思琪
分配工作并完成成果汇总
曹辰
遥感图像的点特征提取，并解释语句程序的含义
高晓鹏
遥感图像的线特征提取，并解释语句程序的含义
李鹏年
通过点、线特征提取对特征提取进一步加深理解
使用仪器设备：
笔记本电脑、MATLAB7.0
三、综合训练内容
1、点特征提取
（1）Moravec点特征提取算子
clc
clear
I=imread（'C:
\Users\Administrator\Desktop\145.bmp'）;
originalmap=I;
ifisrgb（I）
I=rgb2gray（I）;%将彩色图像转换为灰度图像
end
I=double（I）;%转换为双精度图像
[row,column]=size（I）;
w1=floor（5/2）;%确定窗口大小
w2=5;
step=w1;
ir=7000;%设置阈值
num=0;
%得到特征点的坐标
xc=[];
yc=[];
tic;
fori=w1+1:
1:
row-w1-1
ifi+w2>row-w1-1
continue;
end
forj=w1+1:
1:
column-w1-1
ifj+w2>column-w1-1
continue;
end
interest=0;%兴趣值
v1=0;v2=0;v3=0;v4=0;
fork=-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]）;
ifinterest>ir%//最大兴趣值与阈值比较
num=num+1;
pointx（num）=j;
pointy（num）=i;
yc=[yc,i];
xc=[xc,j];
end
end
end
%图像中显示特征点
figure（4）;
imshow（originalmap）;
holdon;
plot（xc,yc,'R*'）;
axison;
title（'图像中显示特征点'）;
xlabel（'图像的列数'）;
ylabel（'图像的行数'）;
toc;
t=toc;
disp（['本程序的运行时间为',num2str（t）,'秒。
']）;
（2）Forstner点特征提取算子
clc
clear
I=imread（'C:
\Users\LBS\Desktop\MoravecMATLAB\girl.bmp'）;
originalmap=I;
ifisrgb（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;
fori=k+1:
row-k
forj=k+1:
column-k
%协方差矩阵
G0=0;G2=0;G1=0;
form=i-k:
i+k-1
forn=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;
iftrN==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;
fori=k:
w2:
row-w1
ifi+w2>row-w1
continue;
end
forj=k:
w2:
column-w1
ifj+w2>column-w1
continue;
end
max_w=0
max_x=0;max_y=0;
formh=1:
w2
formw=1:
w2
if（max_wmax_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）;
holdon;
plot（xc,yc,'R*'）;
axison;
title（'图像中显示特征点'）;
xlabel（'图像的列数'）;
ylabel（'图像的行数'）;
2、线特征提取
（1）Robert线特征提取算子
%%%%Robert算子
function[ans]=robert（I）;
ans=zeros（size（I））;
Isize=size（I）;
fori=1:
Isize
（1）-1
forj=1:
Isize
（2）-1
aa=I（i,j）-I（i+1,j）;
bb=I（i,j）-I（i,j+1）;
ans（i,j）=abs（aa）+abs（bb）;
ifans（i,j）>50
ans（i,j）=255;
else
ans（i,j）=0;
end
end
End
（2）Soble线特征提取算子
%%%%sobel算子
function[ans]=sobel（I）;
ans=zeros（size（I））;
Isize=size（I）;
I=double（I）;
fx=0;
fy=0;
fori=2:
Isize
（1）-1
forj=2:
Isize
（2）-1
fx=（I（i-1,j+1）+2*I（i,j+1）+I（i+1,j+1））-（I（i-1,j-1）+2*I（i,j-1）+I（i+1,j-1））;
fy=（I（i+1,j-1）+2*I（i+1,j）+I（i+1,j+1））-（I（i-1,j-1）+2*I（i-1,j）+I（i-1,j+1））;
ans（i,j）=sqrt（fx*fx+fy*fy）;%
ifans（i,j）>4%160
ans（i,j）=255;
else
ans（i,j）=0;
end
end
end
%tp=min（min（ans））;
%ans=ans+tp;
ans=uint8（ans）;
（3）prewiit线特征提取算子
%%%%prewitt算子
function[ans]=prewitt（I）;
ans=zeros（size（I））;
Isize=size（I）;
I=double（I）;
fx=0;
fy=0;
fori=2:
Isize
（1）-1
forj=2:
Isize
（2）-1
fx=（I（i-1,j+1）+I（i,j+1）+I（i+1,j+1））-（I（i-1,j-1）+I（i,j-1）+I（i+1,j-1））;
fy=（I（i+1,j-1）+I（i+1,j）+I（i+1,j+1））-（I（i-1,j-1）+I（i-1,j）+I（i-1,j+1））;
ans（i,j）=sqrt（fx*fx+fy*fy）;
ifans（i,j）>130
ans（i,j）=255;
else
ans（i,j）=0;
end
end
end
ans=uint8（ans）;
（4）拉普拉斯线特征提取算子
I=imread（'C:
\Users\Administrator\Desktop\123.bmp'）;
subplot（2,2,1）;
imshow（I）;
title（'原始图像'）;
I1=rgb2gray（I）;
subplot（2,2,2）;
imshow（I1）;
title（'灰度图像'）;
I2=edge（I1,'log'）;
subplot（2,2,3）;
imshow（I2）;
title（'log算子分割结果'）;
（5）高斯—拉普拉斯线特征提取算子
%%%%marr算子
function[ans]=marr（I）;
ans=zeros（size（I））;
Isize=size（I）;
I=double（I）;
fx=0;
fy=0;
fori=3:
Isize
（1）-2
forj=3:
Isize
（2）-2
ans（i,j）=16*I（i,j）-2*I（i-1,j）-2*I（i+1,j）-2*I（i,j-1）-2*I（i,j+1）-I（i-1,j-1）-I（i-1,j+1）-I（i+1,j-1）-I（i+1,j+1）-I（i-2,j）-I（i+2,j）-I（i,j-2）-I（i,j+2）;
ifans（i,j）>-120
ans（i,j）=0;
else
ans（i,j）=255;
end
end
end
ans=uint8（ans）;
（6）Hough直线提取算子
%%%%canny算子
function[ans]=canny（I）;
gauss_I=I;
Isize=size（I）;
ans=zeros（size（I））;
dir=zeros（size（I））;
I=double（I）;
gauss_I=double（gauss_I）;
fx=0;
fy=0;
fori=2:
Isize
（1）-1
forj=2:
Isize
（2）-1
fx=gauss_I（i,j）+gauss_I（i,j+1）-gauss_I（i+1,j）-gauss_I（i+1,j+1）;
fy=gauss_I（i,j）+gauss_I（i+1,j）-gauss_I（i,j+1）-gauss_I（i+1,j+1）;
ans（i,j）=sqrt（fx*fx+fy*fy）;
dir（i,j）=atan（fy/fx）;
end
end
fori=2:
Isize
（1）-1
forj=2:
Isize
（2）-1
ifdir（i,j）>=-pi/8&dir（i,j）ifans（i,j）<=ans（i,j-1）|ans（i,j）<=ans（i,j+1）
ans（i,j）=0;
end
end
ifdir（i,j）>=pi/8&dir（i,j）<3*pi/8
ifans（i,j）<=ans（i-1,j+1）|ans（i,j）<=ans（i+1,j-1）
ans（i,j）=0;
end
end
ifdir（i,j）>=3*pi/8|dir（i,j）<-3*pi/8
ifans（i,j）<=ans（i-1,j）|ans（i,j）<=ans（i+1,j）
ans（i,j）=0;
end
end
ifdir（i,j）<-pi/8&dir（i,j）>=3*pi/8
ifans（i,j）<=ans（i-1,j-1）|ans（i,j）<=ans（i+1,j+1）
ans（i,j）=0;
end
end
ifans（i,j）<40
ans（i,j）=0;
else
ans（i,j）=255;
end
end
end
ans=uint8（ans）;
四、综合训练成果
（1）Moravec点特征提取算子
图1
（2）Forstner点特征提取算子
图2
（3）Robert线特征提取算子
图3
（4）Soble线特征提取算子
图4
（5）Prewiit线特征提取算子
图5
（6）拉普拉斯线特征提取算子
图6
（7）高斯—拉普拉斯线特征提取算子
图7
（8）Hough直线提取算子
图8
五、综合训练体会
通过本次综合训练的学习，我对图像的特征提取有了深刻的理解。
通过进行MATLAB程序语句的研究发现算法的优缺点。
对于不同的算子有着不同的矩阵计算模板，在进行计算时会发现因为算法的优化不同会造成遥感图像提取质量产生很大差异。
我对比了网上几个论坛较好下载的MATLAB算子代码，通过比较发现对相同的算子在进行计算时对矩阵模板计算的越准确越是会产生较好的结果。
在代码编写时应注意提取线特征时确定好不同的特征线是否为同一条直线，如果不是应注意区分不同的直线。
对于不同的算法之间也存在不同的优化，虽然大部分模板都会计算梯度但是在对梯度进行应用时会有很大差别，算法之间最大的差别就是卷积模板的不同。
在这次综合训练之中我对MATLAB语句的认识有了进一步的提高。
指
导
教
师
意
见
成绩评定：
指导教师签字：
年月日
教学
单位意见
负责人签字：
（单位盖章）
年月日
备注
注：
综合训练结束时，由实习学生填写本表后，交指导教师和教学单位签署意见，最后交所在教学单位归档保管。
升级会员 