遥感影像特征点提取
一、 基于Moravec算子的特征点提取
1. Moravec算子的原理及算法公式
该算子是通过逐像元量测与其邻元的灰度差，搜索相邻像元之间具有高反差的点，具体方法有以下几种。
（1）计算各像元的有利值，如图所示，在5×5的窗口内沿着图示四个方向分别计算相邻像元间灰度差之平方和V1，V2，V3，及V4，取其中最小值作为该像元的有利值： IVmin=min{V1,V2,V3,V4}
其中： V1=
V2=
V3=
V4= ∑(Giii,j-Gi+1,j)2-Gi,j+1)2-Gi+1,j+1)2∑(Gii,j∑(G∑(G
ii,j2-G)i,ji+1,j-1式中， i=m-k, ,m+k-1; j=n-k, ,n+k-1;
k
=W/2。
Gi,j代表像元Pi,j的灰度值，W为以像元计的窗口大小，如图所示，W=5，m,n为像元在整块影像中位置序号。
（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);
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算子用高斯函数代替二值窗口函数，对离中心点越近的像素赋予越大的权重，以减少噪声影响。
1 -(x2+y2)2σ2w(x,y)=e2πσ2
（2）Moravec算子只考虑了每隔45度方向，Harris算子用Taylor展开去近似任意方向。将图像窗口平移[u,v]产生灰度变化E(u,v)。
E(u,v)=w(x,y)[I(x+u,y+v)-I(x,y)]2
x,y ∑
由
I（x+u,y+v）=I(x,y)+Ixu+Iyv+O(u2,v2)
得 E(u,v)=w(x,y)[Ixu+Iyv+O(u2,v2)]2
x,y 2⎡IIxIy⎤⎡u⎤2x Ixu+Iyv=[u,v]⎢2⎥⎢⎥III ⎢y⎥⎣xy⎦⎣v⎦
于是对于局部微小的移动量[u,v]，可以近似得到下面的表达：
⎡u⎤ E（u,v）=[u,v]M⎢⎥v其中M是2×2的矩阵，可由图像的导数求得： ⎣⎦
⎡Ix2IxIy⎤ M=w(x,y)⎢2⎥IIIx,y⎢xy⎥y⎦ ⎣∑[]∑
wx , y ) 为高斯函数。 式中，I x 为x方向的差分，I y 为y方向的差分， (
（3）Harris采用了一种新的角点判断方法。通过M的两个特征值λ1，λ2的大小对图像点进行分类。
但是解特征向量需要比较多的计算量，且两个特征值的和等于矩阵M的迹，两个特征值的积等于矩阵M的行列式。所以用下式来判定角点质量。（K常取0.04—0.06） 2R=detM-k(traceM)
（4）Harris算法总结
1：对每一像素点计算相关矩阵M
A=w(x,y)⊗Ix2B=w(x,y)⊗Iy2C=D=w(x,y)⊗（IxIy）⎡AB⎤M=⎢⎥⎣CD⎦2：计算每像素点的Harris角点响应。 2R=（AB-CD）-k(A+B)2
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;
for i=boundary:nrow-boundary+1
for j=boundary:ncol-boundary+1
nlike=0; %相似点个数
if Image(i-1,j-1)>Image(i,j)-t && Image(i-1,j-1)
nlike=nlike+1;
end
if Image(i-1,j)>Image(i,j)-t && Image(i-1,j)
nlike=nlike+1;
end
if Image(i-1,j+1)>Image(i,j)-t && Image(i-1,j+1)
nlike=nlike+1;
end
if Image(i,j-1)>Image(i,j)-t && Image(i,j-1)
nlike=nlike+1;
end
if Image(i,j+1)>Image(i,j)-t && Image(i,j+1)
nlike=nlike+1;
end
if Image(i+1,j-1)>Image(i,j)-t && Image(i+1,j-1)
nlike=nlike+1;
end
if Image(i+1,j)>Image(i,j)-t && Image(i+1,j)
nlike=nlike+1;
end
if Image(i+1,j+1)>Image(i,j)-t && Image(i+1,j+1)
nlike=nlike+1;
end
if nlike>=2 && nlike
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对采集的这几幅图像来说是一个比较合理的经验值 for i = boundary:nrow-boundary+1
for j = boundary:ncol-boundary+1
if Corner(i,j)==1 %只关注候选点
M = [A(i,j) C(i,j);
C(i,j) B(i,j)];
CRF(i,j) = det(M)-t*(trace(M))^2;
if CRF(i,j) > CRFmax
CRFmax = CRF(i,j);
end;
end
end;
end;
%CRFmax
count = 0; % 用来记录角点的个数
t=0.01;
% 下面通过一个3*3的窗口来判断当前位置是否为角点
for i = boundary:nrow-boundary+1
for j = boundary:ncol-boundary+1
if Corner(i,j)==1 %只关注候选点的八邻域
if CRF(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(Image); % display Intensity Image
hold on;
% toc(t1)
for i=boundary:nrow-boundary+1
for j=boundary:ncol-boundary+1
column_ave=0;
row_ave=0;
k=0;
if Corner(i,j)==1
for x=i-3:i+3 %7*7邻域
for y=j-3:j+3
if Corner(x,y)==1
% 用算数平均数作为角点坐标，如果改用几何平均数求点的平均坐标，对角点的提取意义不大
row_ave=row_ave+x;
column_ave=column_ave+y;
k=k+1;
end
end
end
end
if k>0 %周围不止一个角点
plot( column_ave/k,row_ave/k ,'g.');
end
end;
end;
三 运行结果
1001特征点 1002特征点
二、
（二）
（三）Susan算子
一 算法公式 （1）借助图3-1来解释Susan检测的原理，其中图片是白色背景，有一个颜色比较暗淡的矩形。用一个圆形模板在图像上移动，若模板内的像素灰度与模板中心的像素（被称为核Nucleus）灰度值小于一定的阈值，则认为该点与核Nucleus具有相同的灰度，满足该条件的像素组成的区域就称为USAN。在图片上有5个圆形区域。圆形区域表示的是掩码区域。把圆形区域内的每一个位置的像素值与圆心处的像素值相比较，那么圆中的的像素可以分为两类，一类是像素值与圆心处的像素值相近的，
另一类是像素值与圆心的处的像素值相差比较大的。
图3-1 图3-2
如果将模板中各个像素的灰度都与模板中心的核像素的灰度进行比较，那么就会发现总有一部分模板区域和灰度与核像素的灰度相同或相似，这部分区域可以称为USAN。USAN区域包含很多与图像结构有关的信息。利用这种区域的尺寸、重心、二阶矩的分析，可以得到图像中的角点，边缘等信息。从上图所示，当核像素处在图像中的灰度一致区域时，USAN的面积会达到最大。第e个模板就是属于这种情况。
（2）Susan进行角点检测时，遵循了常规的思路：使用一个窗口在图像上逐点滑动，在每一个位置上计算出一个角点量，再进行局部极大值抑制，得到最终的角点。我们这里使用的窗口是圆形窗口，最小的窗口是3×3的，此次使用的是37个像素的圆形窗口，如图3-2。
（3）在角点检测中，有两种类型的阈值，一种用来约束角点的数量，另一种用来约束角点的质量。当然，一个阈值不能完全做到只影响质量或数量，只是会有一个侧重点。阈值g是角点质量的，尽管也会影响数量，但是相对来说更侧重于影响质量（角点的形状）。例如，g值减小，那么Susan会更加侧重于检测到更加“尖锐”的角点，所以，可以更加自己的实际需求来确定阈值g；而阈值t，是角点的数量，当t减小时，会检测到更多的角点，所以，阈值t可以在不影响角点质量的情况下，控制检测到的角点的数量，如果图像的对比度比较低，可以修改t值以适应变化。
下面简单叙述下利用Susan算子检测角点的步骤：
1:利用圆形模板遍历图像，计算每点处的USAN值。
2：设置一阈值g，一般取值为1/2(Max(n)）， 也即取值为USAN最大值的一半，进行阈值化，得到角点响应。
3：使用非极大值抑制来寻找角点。
通过上面的方式得到的角点，存在很大伪角点。为了去除伪角点，Susan算子可以由以下方法实现：
1:计算USAN区域的重心，然后计算重心和模板中心的距离，如果距离较小则不是正确的角点；
2：判断USAN区域的重心和模板中心的连线所经过的像素都是否属于USAN区域的像素，如果属于那么这个模板中心的点就是角点。
二 MATLAB代码
clear all;close all;clc;
img=imread('1001.jpg');
img=rgb2gray(img);
imshow(img); [m n]=size(img);
img=double(img);
t=45; %模板中心像素灰度和周围灰度差别的阈值，自己设置 usan=[]; %当前像素和周围在像素差别在t以下的个数
%这里用了37个像素的模板
for i=4:m-3 %没有在外围扩展图像，最终图像会缩小
for j=4:n-3
tmp=img(i-3:i+3,j-3:j+3); %先构造7*7的模板，49个像素
c=0;
for p=1:7 for q=1:7
if (p-4)^2+(q-4)^2
if abs(img(i,j)-tmp(p,q))
end
end
end
end
usan=[usan c];
end
end
g=2*max(usan)/3; %确定角点提取的数量，值比较高时会提取出边缘，自己设置
for i=1:length(usan)
if usan(i)
usan(i)=g-usan(i);
else
usan(i)=0;
end
end imgn=reshape(usan,[n-6,m-6])';
figure; imshow(imgn)
%非极大抑制
[m n]=size(imgn);
re=zeros(m,n); for i=2:m-1
for j=2:n-1
if imgn(i,j)>max([max(imgn(i-1,j-1:j+1))
max(imgn(i+1,j-1:j+1))]);
re(i,j)=1;
else
re(i,j)=0;
end end
end
figure;
imshow(re==1); imgn(i,j-1) imgn(i,j+1)
三 运行结果
1001特征点 1002特征点
（四）三种方法的比较：
（1）Moravec算子对斜边缘的响应很强，因为只考虑了每隔45度的方向变化，而没有在全部的方向上进行考虑；同时由于窗口函数是一个二值函数,不管像素点离中心点的距离，赋于一样的权重，因此对噪声响应也一般，最终对角点的定位较准确。
（2）Harris算子是一种有效的点特征提取算子,其优点总结起来有:
1：计算简单:Harris算子中只用到灰度的一阶差分以及滤波,操作简单。
2:提取的点特征均匀而且合理:Harris算子对图像中的每个点都计算其兴趣值,然后在邻域中选择最优点。实验表明,在纹理信息丰富的区域,Harris算子可以提取出大量有用的特征点,而在纹理信息少的区域,提取的特征点则较少。
3:稳定:Harris算子的计算公式中只涉及到一阶导数，因此对图像旋转、灰度变化、噪声影响和视点变换不敏感,它也是比较稳定的一种点特征提取算子。
Harris算子的局限性有：
1:它对尺度很敏感，不具有尺度不变性。
2:提取的角点是像素级的。
（3）Susan算子的特点有：
1:在用Susan算子对边缘和角点进行检测时不需要计算微分，这使得Susan算子对噪声更加鲁棒。
2:Susan检测算子能提供不依赖于模板尺寸的边缘精度。换句话说，最小USAN区域面积的计算是个相对的概念，与模版尺寸无关，所以Susan边缘算子的性能不受模版尺寸影响。 3:控制参数的选择很简单，且任意性小，容易实现自动化选取。
（五）心得体会
通过本次实验，我对特征点提取方法的计算原理及实现过程有了深刻的认识，在本次实验中，我遇到了很多困难，但是在同学们的帮助下，我们互相商讨，这些问题都一一得到了解决，在解决困难的过程中的编程能力得到了提高，对其所涉及到的知识的印象也得到了加深。
总之，感谢老师给了我这次锻炼自己的机会，之后还要继续学习研究MATLAB，提升自己的编程能力。
遥感影像特征点提取
一、 基于Moravec算子的特征点提取
1. Moravec算子的原理及算法公式
该算子是通过逐像元量测与其邻元的灰度差，搜索相邻像元之间具有高反差的点，具体方法有以下几种。
（1）计算各像元的有利值，如图所示，在5×5的窗口内沿着图示四个方向分别计算相邻像元间灰度差之平方和V1，V2，V3，及V4，取其中最小值作为该像元的有利值： IVmin=min{V1,V2,V3,V4}
其中： V1=
V2=
V3=
V4= ∑(Giii,j-Gi+1,j)2-Gi,j+1)2-Gi+1,j+1)2∑(Gii,j∑(G∑(G
ii,j2-G)i,ji+1,j-1式中， i=m-k, ,m+k-1; j=n-k, ,n+k-1;
k
=W/2。
Gi,j代表像元Pi,j的灰度值，W为以像元计的窗口大小，如图所示，W=5，m,n为像元在整块影像中位置序号。
（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);
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算子用高斯函数代替二值窗口函数，对离中心点越近的像素赋予越大的权重，以减少噪声影响。
1 -(x2+y2)2σ2w(x,y)=e2πσ2
（2）Moravec算子只考虑了每隔45度方向，Harris算子用Taylor展开去近似任意方向。将图像窗口平移[u,v]产生灰度变化E(u,v)。
E(u,v)=w(x,y)[I(x+u,y+v)-I(x,y)]2
x,y ∑
由
I（x+u,y+v）=I(x,y)+Ixu+Iyv+O(u2,v2)
得 E(u,v)=w(x,y)[Ixu+Iyv+O(u2,v2)]2
x,y 2⎡IIxIy⎤⎡u⎤2x Ixu+Iyv=[u,v]⎢2⎥⎢⎥III ⎢y⎥⎣xy⎦⎣v⎦
于是对于局部微小的移动量[u,v]，可以近似得到下面的表达：
⎡u⎤ E（u,v）=[u,v]M⎢⎥v其中M是2×2的矩阵，可由图像的导数求得： ⎣⎦
⎡Ix2IxIy⎤ M=w(x,y)⎢2⎥IIIx,y⎢xy⎥y⎦ ⎣∑[]∑
wx , y ) 为高斯函数。 式中，I x 为x方向的差分，I y 为y方向的差分， (
（3）Harris采用了一种新的角点判断方法。通过M的两个特征值λ1，λ2的大小对图像点进行分类。
但是解特征向量需要比较多的计算量，且两个特征值的和等于矩阵M的迹，两个特征值的积等于矩阵M的行列式。所以用下式来判定角点质量。（K常取0.04—0.06） 2R=detM-k(traceM)
（4）Harris算法总结
1：对每一像素点计算相关矩阵M
A=w(x,y)⊗Ix2B=w(x,y)⊗Iy2C=D=w(x,y)⊗（IxIy）⎡AB⎤M=⎢⎥⎣CD⎦2：计算每像素点的Harris角点响应。 2R=（AB-CD）-k(A+B)2
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;
for i=boundary:nrow-boundary+1
for j=boundary:ncol-boundary+1
nlike=0; %相似点个数
if Image(i-1,j-1)>Image(i,j)-t && Image(i-1,j-1)
nlike=nlike+1;
end
if Image(i-1,j)>Image(i,j)-t && Image(i-1,j)
nlike=nlike+1;
end
if Image(i-1,j+1)>Image(i,j)-t && Image(i-1,j+1)
nlike=nlike+1;
end
if Image(i,j-1)>Image(i,j)-t && Image(i,j-1)
nlike=nlike+1;
end
if Image(i,j+1)>Image(i,j)-t && Image(i,j+1)
nlike=nlike+1;
end
if Image(i+1,j-1)>Image(i,j)-t && Image(i+1,j-1)
nlike=nlike+1;
end
if Image(i+1,j)>Image(i,j)-t && Image(i+1,j)
nlike=nlike+1;
end
if Image(i+1,j+1)>Image(i,j)-t && Image(i+1,j+1)
nlike=nlike+1;
end
if nlike>=2 && nlike
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对采集的这几幅图像来说是一个比较合理的经验值 for i = boundary:nrow-boundary+1
for j = boundary:ncol-boundary+1
if Corner(i,j)==1 %只关注候选点
M = [A(i,j) C(i,j);
C(i,j) B(i,j)];
CRF(i,j) = det(M)-t*(trace(M))^2;
if CRF(i,j) > CRFmax
CRFmax = CRF(i,j);
end;
end
end;
end;
%CRFmax
count = 0; % 用来记录角点的个数
t=0.01;
% 下面通过一个3*3的窗口来判断当前位置是否为角点
for i = boundary:nrow-boundary+1
for j = boundary:ncol-boundary+1
if Corner(i,j)==1 %只关注候选点的八邻域
if CRF(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(Image); % display Intensity Image
hold on;
% toc(t1)
for i=boundary:nrow-boundary+1
for j=boundary:ncol-boundary+1
column_ave=0;
row_ave=0;
k=0;
if Corner(i,j)==1
for x=i-3:i+3 %7*7邻域
for y=j-3:j+3
if Corner(x,y)==1
% 用算数平均数作为角点坐标，如果改用几何平均数求点的平均坐标，对角点的提取意义不大
row_ave=row_ave+x;
column_ave=column_ave+y;
k=k+1;
end
end
end
end
if k>0 %周围不止一个角点
plot( column_ave/k,row_ave/k ,'g.');
end
end;
end;
三 运行结果
1001特征点 1002特征点
二、
（二）
（三）Susan算子
一 算法公式 （1）借助图3-1来解释Susan检测的原理，其中图片是白色背景，有一个颜色比较暗淡的矩形。用一个圆形模板在图像上移动，若模板内的像素灰度与模板中心的像素（被称为核Nucleus）灰度值小于一定的阈值，则认为该点与核Nucleus具有相同的灰度，满足该条件的像素组成的区域就称为USAN。在图片上有5个圆形区域。圆形区域表示的是掩码区域。把圆形区域内的每一个位置的像素值与圆心处的像素值相比较，那么圆中的的像素可以分为两类，一类是像素值与圆心处的像素值相近的，
另一类是像素值与圆心的处的像素值相差比较大的。
图3-1 图3-2
如果将模板中各个像素的灰度都与模板中心的核像素的灰度进行比较，那么就会发现总有一部分模板区域和灰度与核像素的灰度相同或相似，这部分区域可以称为USAN。USAN区域包含很多与图像结构有关的信息。利用这种区域的尺寸、重心、二阶矩的分析，可以得到图像中的角点，边缘等信息。从上图所示，当核像素处在图像中的灰度一致区域时，USAN的面积会达到最大。第e个模板就是属于这种情况。
（2）Susan进行角点检测时，遵循了常规的思路：使用一个窗口在图像上逐点滑动，在每一个位置上计算出一个角点量，再进行局部极大值抑制，得到最终的角点。我们这里使用的窗口是圆形窗口，最小的窗口是3×3的，此次使用的是37个像素的圆形窗口，如图3-2。
（3）在角点检测中，有两种类型的阈值，一种用来约束角点的数量，另一种用来约束角点的质量。当然，一个阈值不能完全做到只影响质量或数量，只是会有一个侧重点。阈值g是角点质量的，尽管也会影响数量，但是相对来说更侧重于影响质量（角点的形状）。例如，g值减小，那么Susan会更加侧重于检测到更加“尖锐”的角点，所以，可以更加自己的实际需求来确定阈值g；而阈值t，是角点的数量，当t减小时，会检测到更多的角点，所以，阈值t可以在不影响角点质量的情况下，控制检测到的角点的数量，如果图像的对比度比较低，可以修改t值以适应变化。
下面简单叙述下利用Susan算子检测角点的步骤：
1:利用圆形模板遍历图像，计算每点处的USAN值。
2：设置一阈值g，一般取值为1/2(Max(n)）， 也即取值为USAN最大值的一半，进行阈值化，得到角点响应。
3：使用非极大值抑制来寻找角点。
通过上面的方式得到的角点，存在很大伪角点。为了去除伪角点，Susan算子可以由以下方法实现：
1:计算USAN区域的重心，然后计算重心和模板中心的距离，如果距离较小则不是正确的角点；
2：判断USAN区域的重心和模板中心的连线所经过的像素都是否属于USAN区域的像素，如果属于那么这个模板中心的点就是角点。
二 MATLAB代码
clear all;close all;clc;
img=imread('1001.jpg');
img=rgb2gray(img);
imshow(img); [m n]=size(img);
img=double(img);
t=45; %模板中心像素灰度和周围灰度差别的阈值，自己设置 usan=[]; %当前像素和周围在像素差别在t以下的个数
%这里用了37个像素的模板
for i=4:m-3 %没有在外围扩展图像，最终图像会缩小
for j=4:n-3
tmp=img(i-3:i+3,j-3:j+3); %先构造7*7的模板，49个像素
c=0;
for p=1:7 for q=1:7
if (p-4)^2+(q-4)^2
if abs(img(i,j)-tmp(p,q))
end
end
end
end
usan=[usan c];
end
end
g=2*max(usan)/3; %确定角点提取的数量，值比较高时会提取出边缘，自己设置
for i=1:length(usan)
if usan(i)
usan(i)=g-usan(i);
else
usan(i)=0;
end
end imgn=reshape(usan,[n-6,m-6])';
figure; imshow(imgn)
%非极大抑制
[m n]=size(imgn);
re=zeros(m,n); for i=2:m-1
for j=2:n-1
if imgn(i,j)>max([max(imgn(i-1,j-1:j+1))
max(imgn(i+1,j-1:j+1))]);
re(i,j)=1;
else
re(i,j)=0;
end end
end
figure;
imshow(re==1); imgn(i,j-1) imgn(i,j+1)
三 运行结果
1001特征点 1002特征点
（四）三种方法的比较：
（1）Moravec算子对斜边缘的响应很强，因为只考虑了每隔45度的方向变化，而没有在全部的方向上进行考虑；同时由于窗口函数是一个二值函数,不管像素点离中心点的距离，赋于一样的权重，因此对噪声响应也一般，最终对角点的定位较准确。
（2）Harris算子是一种有效的点特征提取算子,其优点总结起来有:
1：计算简单:Harris算子中只用到灰度的一阶差分以及滤波,操作简单。
2:提取的点特征均匀而且合理:Harris算子对图像中的每个点都计算其兴趣值,然后在邻域中选择最优点。实验表明,在纹理信息丰富的区域,Harris算子可以提取出大量有用的特征点,而在纹理信息少的区域,提取的特征点则较少。
3:稳定:Harris算子的计算公式中只涉及到一阶导数，因此对图像旋转、灰度变化、噪声影响和视点变换不敏感,它也是比较稳定的一种点特征提取算子。
Harris算子的局限性有：
1:它对尺度很敏感，不具有尺度不变性。
2:提取的角点是像素级的。
（3）Susan算子的特点有：
1:在用Susan算子对边缘和角点进行检测时不需要计算微分，这使得Susan算子对噪声更加鲁棒。
2:Susan检测算子能提供不依赖于模板尺寸的边缘精度。换句话说，最小USAN区域面积的计算是个相对的概念，与模版尺寸无关，所以Susan边缘算子的性能不受模版尺寸影响。 3:控制参数的选择很简单，且任意性小，容易实现自动化选取。
（五）心得体会
通过本次实验，我对特征点提取方法的计算原理及实现过程有了深刻的认识，在本次实验中，我遇到了很多困难，但是在同学们的帮助下，我们互相商讨，这些问题都一一得到了解决，在解决困难的过程中的编程能力得到了提高，对其所涉及到的知识的印象也得到了加深。
总之，感谢老师给了我这次锻炼自己的机会，之后还要继续学习研究MATLAB，提升自己的编程能力。