遥感提取特征点.docx
《遥感提取特征点.docx》由会员分享，可在线阅读，更多相关《遥感提取特征点.docx（14页珍藏版）》请在冰豆网上搜索。
遥感提取特征点
遥感影像特征点提取
基于Moravec算子的特征点提取
1.Moravec算子的原理及算法公式
该算子是通过逐像元量测与其邻元的灰度差，搜索相邻像元之间具有高反差的点，具体
方法有以下几种。
（1）计算各像元的有利值，如图所示，在5X5的窗口内沿着图示四个方向分别计算相邻像元间灰度差之平方和Vi,V2,V3,及V4,取其中最小值作为该像元的有利值：
IVmin二mi
Vi='（Gi,j-Gii,j）
i
V2八（Gi,j—Gi,ji）2
i
2
V3八（Gi,j-Gi.i,j1）
i
V4=為（Gi,j-Gi.i,j-i）2
i
i=m-k,,mk—i;
j=n-k,,nk-i;k=W/2。
Gi,j代表像元R,j的灰度值，W为以像元计的窗口大小，如图所示，W=5,m,n为像
元在整块影像中位置序号。
（2）给定一个阈值，确定待定点的有利点。
如果有利值大于给定的阈值，则将该像元作为候选点。
阈值一般为经验值。
（3）抑制局部非最大。
在一定大小窗口内（例如5X5,7X7,,9X9像元等），将上
一步所选的候选点与其周围的候选点比较，若该像元的有利非窗口中最大值，则去掉；否则,
该像元被确定为特征点，这一步的目的在于避免纹理丰富的区域产生束点，用于抑制局部非
最大的窗口大小取决于所需的有利点密度。
综上所述，Moravec算子是在四个主要方向上选择具有最大一最小灰度方差的点作为特征点。
2.基于MATLAB的算法编程
clearall;closeall;clc
img=double（imread（'1001.jpg'））;
[hw]=size（img）;imshow（img,[]）imgn=zeros（h,w）;
n=4;
fory=1+n:
h-n
forx=1+n:
w-n
sq=img（y_n:
y+n,x_n:
x+n）;
V=zeros（1,4）;
fori=2:
2*n+1%垂直，水平，对角，反对角四个万向领域灰度差的平万和
V
（1）=V
（1）+（sq（i,n+1）-sq（i-1,门+1））人2;
V
（2）=V
（2）+（sq（n+1,i）-sq（n+1,i-1））A2;
W3）=V（3）+（sq（i,i）-sq（i-1,i-1））A2;
V（4）=V（4）+（sq（i,（2*n+1）-（i-1））-sq（i-1,（2*n+1）-（i-2）））A2;
end
pix=min（V）;%四个方向中选最小值
imgn（y,x）=pix;
end
end
T=mean（imgn（:
））;%设阈值，小于均值置零
ind=find（imgnimgn（ind）=0;
fory=1+n:
h-n%选局部最大且非零值作为特征点
forx=1+n:
w-n
sq=imgn（y-n:
y+n,x-n:
x+n）;
ifmax（sq（:
））==imgn（y,x）&&imgn（y,x）~=0
img（y,x）=255;
end
end
end
figure;
imshow（img,[]）;
3・运行结果
1002特征点
1001特征点
Harris角点检测算子
1算法公式
（1）Harris算子用高斯函数代替二值窗口函数，对离中心点越近的像素赋予越大的权重，以减少噪声影响。
E（u,v）口[u,v]M
其中M是2X2的矩阵，可由图像的导数求得：
但是解特征向量需要比较多的计算量，且两个特征值的和等于矩阵M的迹，两个
特征值的积等于矩阵M的行列式。
所以用下式来判定角点质量。
（K常取0.04—0.06）
R=detM-k（traceM）
（4）Harris算法总结
1对每一像素点计算相关矩阵M
A二w（x,y）：
I
B=w（x,y）_Iy
C=D=w（x,y）i（I
M〔CD
2:
计算每像素点的Harris角点响应。
R（（AB-CD）-k（AB）
3:
在wxw范围内寻找极大值点，若Harris角点响应大于阀值，则视为角点。
Harris算子对灰度的平移是不变的，因为只有差分，对旋转也有不变性，但是对尺度很敏感，在一个尺度下是角点，在另一个尺度下可能就不是了。
二MATLAB代码
clear;
%读取图像
Image=imread（'1001.jpg'）;
Image=im2uint8（rgb2gray（lmage））;
dx=[-101;-101;-101];%dx:
横向Prewitt差分模版
Ix2=filter2（dx,Image）.A2;
Iy2=filter2（dx',lmage）.A2;
Ixy=filter2（dx,Image）.*filter2（dx',Image）;
%生成9*9高斯窗口。
窗口越大，探测到的角点越少。
h=fspecial（'gaussian',9,2）;
A=filter2（h,Ix2）;%用高斯窗口差分Ix2得到A
B=filter2（h,Iy2）;
C=filter2（h,Ixy）;
nrow=size（Image,1）;
ncol=size（Image,2）;
Corner=zeros（nrow,ncol）;
%矩阵Corner用来保存候选角点位置，初值全零，值为1的点是角点
%真正的角点在137和138行由（row_ave,column_ave）得到
%参数t:
点（i,j）八邻域的相似度”参数，只有中心点与邻域其他八个点的像素值之差在
%（-t,+t）之间，才确认它们为相似点，相似点不在候选角点之列
t=20;
%并没有全部检测图像每个点，而是除去了边界上boundary个像素，
%因为我们感兴趣的角点并不出现在边界上
boundary=8;
fori=boundary:
nrow-boundary+1
forj=boundary:
ncol-boundary+1
nlike=0;%相似点个数
iflmage（i-1,j-1）>lmage（i,j）-t&&lmage（i_1,j_1）nlike=nlike+1;
end
ifImage（i-1,j）>Image（i,j）-t&&Image（i_1,j）nlike=nlike+1;
end
ifImage（i-1,j+1）>Image（i,j）-t&&Image（i-1,j+1）nlike=nlike+1;
end
ifImage（i,j-1）>Image（i,j）-t&&Image（i,j_1）nlike=nlike+1;
end
ifImage（i,j+1）>Image（i,j）-t&&Image（i,j+1）nlike=nlike+1;
end
ifImage（i+1,j-1）>Image（i,j）-t&&Image（i+1,j-1）nlike=nlike+1;
end
ifImage（i+1,j）>Image（i,j）-t&&Image（i+1,j）nlike=nlike+1;
end
ifImage（i+1,j+1）>Image（i,j）-t&&Image（i+1,j+1）nlike=nlike+1;
end
ifnlike>=2&&nlike<=6
Corner（i,j）=1;%如果周围有0,1,7,8个相似与中心的（i,j）
%那（i,j）就不是角点，所以，直接忽略
end;
end;
end;
CRF=zeros（nrow,ncol）;%CRF用来保存角点响应函数值，初值全零
CRFmax=0;%图像中角点响应函数的最大值，作阈值之用
t=0.05;
%计算CRF
%工程上常用CRF（i,j）=det（M）/trace（M）计算CRF，那么此时应该将下面第105行的
%比例系数t设置大一些，t=0.1对采集的这几幅图像来说是一个比较合理的经验值
fori=boundary:
nrow-boundary+1
forj=boundary:
ncol-boundary+1
ifCorner（i,j）==1%只关注候选点
M=[A（i,j）C（i,j）;
C（i,j）B（i,j）];
CRF（i,j）=det（M）-t*（trace（M））A2;
ifCRF（i,j）>CRFmax
CRFmax=CRF（i,j）;
end;
end
end;
end;
%CRFmax
count=0;%用来记录角点的个数
t=0.01;
%下面通过一个3*3的窗口来判断当前位置是否为角点
fori=boundary:
nrow-boundary+1
forj=boundary:
ncol-boundary+1
ifCorner（i,j）==1%只关注候选点的八邻域
ifCRF（i,j）>t*CRFmax&&CRF（i,j）>CRF（i-1,j-1）……
&&CRF（i,j）>CRF（i-1,j）&&CRF（i,j）>CRF（i-1,j+1）……
&&CRF（i,j）>CRF（i,j-1）&&CRF（i,j）>CRF（i,j+1）……
&&CRF（i,j）>CRF（i+1,j-1）&&CRF（i,j）>CRF（i+1,j）……
&&CRF（i,j）>CRF（i+1,j+1）
count=count+1;%这个是角点，count力口1
else%如果当前位置（i,j）不是角点，则在Corner（i,j）中删除对该候选角点的记录
Corner（i,j）=0;
end;
end;
end;
end;
%disp（'角点个数'）;
%disp（count）
figure,imshow（lmage）;%displayIntensityImage
holdon;
%toc（t1）
fori=boundary:
nrow-boundary+1
forj=boundary:
ncol-boundary+1
column_ave=0;
row_ave=0;
k=0;
ifCorner（i,j）==1
forx=i-3:
i+3%7*7邻域
fory=j-3:
j+3
ifCorner（x,y）==1
%用算数平均数作为角点坐标，如果改用几何平均数求点的平均坐标，对角点的提取
意义不大
row_ave=row_ave+x;
column_ave=column_ave+y;k=k+1;
end
end
end
end
ifk>0%周围不止一个角点
plot（column_ave/k,row_ave/k,'g.'）;end
end;
end;
三运行结果
1001特征点
1002特征点
升级会员 