遥感提取特征点

遥感影像特征点提取

一、 基于Moravec 算子的特征点提取

1. Moravec 算子的原理及算法公式

该算子是通过逐像元量测与其邻元的灰度差，搜索相邻像元之间具有高反差的点，具体方法有以下几种。

（1）计算各像元的有利值，如图所示，在5×5的窗口内沿着图示四个方向分别计算相邻像元间灰度差之平方和V 1，V 2，V 3，及V 4，取其中最小值作为该像元的有利值：

其中：

式中，

j i G ,代表像元j i P ,的灰度值，W 为以像元计的窗口大小，如图所示，n m W ,5，=为像元在整块影像中位置序号。

（2）给定一个阈值，确定待定点的有利点。如果有利值大于给定的阈值，则将该像元作为候选点。阈值一般为经验值。

（3）抑制局部非最大。在一定大小窗口内（例如5×5，7×7，,9×9像元等），将上一步所选的候选点与其周围的候选点比较，若该像元的有利非窗口中最大值，则去掉；否则，该像元被确定为特征点，这一步的目的在于避免纹理丰富的区域产生束点，用于抑制局部非最大的窗口大小取决于所需的有利点密度。

综上所述，Moravec 算子是在四个主要方向上选择具有最大—最小灰度方差的点作为特征点。

2. 基于MATLAB 的算法编程

clear all;close all;clc

img=double(imread('1001.jpg'));

[h w]=size(img);

imshow(img,[])

imgn=zeros(h,w);

n=4;

for y=1+n:h-n

for x=1+n:w-n

sq=img(y-n:y+n,x-n:x+n);

V=zeros(1,4);

}V ,V ,V ,min{V 4321min =IV ∑+-=i j i j i G G

V 21-,1,4)(∑++-=i j i j i G G V 21,1,3)(∑+-=i j i j i G G V 21,,2)(∑+-=i j i j i G G V 2,1,1)(;1,,-+-=k m k m i ;1,,-+-=k n k n j 。2/W k

=

for i=2:2*n+1 %垂直，水平，对角，反对角四个方向领域灰度差的平方和V(1)=V(1)+(sq(i,n+1)-sq(i-1,n+1))^2;

V(2)=V(2)+(sq(n+1,i)-sq(n+1,i-1))^2;

V(3)=V(3)+(sq(i,i)-sq(i-1,i-1))^2;

V(4)=V(4)+(sq(i,(2*n+1)-(i-1))-sq(i-1,(2*n+1)-(i-2)))^2;

end

pix=min(V); %四个方向中选最小值

imgn(y,x)=pix;

end

end

T=mean(imgn(:)); %设阈值，小于均值置零

ind=find(imgn

imgn(ind)=0;

for y=1+n:h-n %选局部最大且非零值作为特征点

for x=1+n:w-n

sq=imgn(y-n:y+n,x-n:x+n);

if max(sq(:))==imgn(y,x) && imgn(y,x)~=0

img(y,x)=255;

end

end

end

figure;

imshow(img,[]);

3.运行结果

1001特征点1002特征点

二、Harris 角点检测算子

1、算法公式

（1）Harris 算子用高斯函数代替二值窗口函数，对离中心点越近的像素赋予越大的权重，以减少噪声影响。

（2）Moravec 算子只考虑了每隔45度方向，Harris 算子用Taylor 展开去近似任意方向。将图像窗口平移[u,v]产生灰度变化E(u,v)。

由

得 于是对于局部微小的移动量[u,v]，可以近似得到下面的表达：

其中M 是2×2的矩阵，可由图像的导数求得： 式中， 为x 方向的差分， 为y 方向的差分，

为高斯函数。 （3）Harris 采用了一种新的角点判断方法。通过M 的两个特征值21λλ，的大小对图像点进行分类。

∑

⎥⎥⎦⎤⎢⎢⎣⎡=y x y y x y x x I I I I I I y x w M ,22),(∑

-++=y x y x I v y u x I y x w v u E ,2

)],(),()[,(),()

,(),(,22v u O v I u I y x I v y u x I y x +++=++）（∑++=y x y x v u O v I u I y x w v u E ,2

22)],()[,(),(⎥⎦

⎤⎢⎣⎡=v u M v u v u E ],[,）（[]⎥⎦

⎤⎢⎣⎡⎥⎥⎦⎤⎢⎢⎣⎡=+v u I I I I I I v u v I u I y y x y x x y x 222],[x I y I ),(y x w 如果1λ和2λ都很小，图像窗口在所有方向上移动都无明显灰度变化 2222)(221),(σπσy x e y x w +-=

但是解特征向量需要比较多的计算量，且两个特征值的和等于矩阵M 的迹，两个特征值的积等于矩阵M 的行列式。所以用下式来判定角点质量。（K 常取0.04—0.06）

（4）Harris 算法总结

1：对每一像素点计算相关矩阵M

2：计算每像素点的Harris 角点响应。

3：在w ×w 范围内寻找极大值点，若Harris 角点响应大于阀值，则视为角点。

Harris 算子对灰度的平移是不变的，因为只有差分，对旋转也有不变性，但是对尺度很敏感，在一个尺度下是角点，在另一个尺度下可能就不是了。

二 MATLAB 代码

clear;

Image = imread('1001.jpg'); % 读取图像

Image = im2uint8(rgb2gray(Image));

dx = [-1 0 1;-1 0 1;-1 0 1]; %dx ：横向Prewitt 差分模版

Ix2 = filter2(dx,Image).^2;

Iy2 = filter2(dx',Image).^2;

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;

2

2)(B A k CD AB R +--=）（2),(x I y x w A ⊗=2),(y I y x w B ⊗=）（y x I I y x w D C ⊗==),(⎥⎦⎤⎢⎣⎡=D C B A M 2)(det traceM k M R -=