遥感提取特征点

遥感影像特征点提取基于Moravec算子的特征点提取1.Moravec算子的原理及算法公式该算子是通过逐像元量测与其邻元的灰度差，搜索相邻像元之间具有高反差的点，具体方法有以下几种。计算各像元的有利值，如图所示，在5X5的窗口内沿着图示四个方向分别计算相邻像元间灰度差之平方和V],v2,v3，及v4，取其中最小值作为该像元的有利值：IV二min{V,V,V,V}min1234其中：二乙(G-G)2i,ji+1,j二工(G-G)2i,ji,j+1二工(G-G)2i,ji+1,j+1二工(G-G)2i,ji+1,j-1式中，i=m一k,…,m+k一1;j=n一k，…,n+k一1;

k=W/2。\1Lkk//■TJ7p\\2G代表像元P的灰度值，W为以像元计的窗口大小，如图所示，W=5,m,n为像i,ji,j元在整块影像中位置序号。给定一个阈值，确定待定点的有利点。如果有利值大于给定的阈值，则将该像元作为候选点。阈值一般为经验值。抑制局部非最大。在一定大小窗口内(例如5X5，7X7，,9X9像元等)，将上一步所选的候选点与其周围的候选点比较，若该像元的有利非窗口中最大值，则去掉；否则，该像元被确定为特征点，这一步的目的在于避免纹理丰富的区域产生束点，用于抑制局部非最大的窗口大小取决于所需的有利点密度。综上所述，Moravec算子是在四个主要方向上选择具有最大一最小灰度方差的点作为特征点。2.基于MATLAB的算法编程clearall;closeall;clcimg=double(imread('1001.jpg'));[hw]=size(img);imshow(img,[])imgn=zeros(h,w);n=4;fory=1+n:h-nforx=1+n:w-nsq=img(y-n:y+n,x-n:x+n);V=zeros(1,4);fori=2:2*n+lfori=2:2*n+l%垂直，水平，对角，反对角四个方向领域灰度差的平方和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)))A2;endpix=min(V);imgn(y,x)=pix;endendT=mean(imgn(:));ind=find(imgnvT);imgn(ind)=0;fory=1+n:h-nforx=1+n:w-n%四个方向中选最小值%四个方向中选最小值%设阈值，小于均值置零%选局部最大且非零值作为特征点ifmax(sq(:))==imgn(y,x)&&imgn(y,x)〜=0img(y,x)=255;endendendfigure;imshow(img,[]);<f■*ij.<T3・运行结果<f■*ij.<T二、Harris角点检测算子1、算法公式(1)Harris算子用高斯函数代替二值窗口函数，对离中心点越近的像素赋予越大的权重，以减少噪声影响。w(x,y)二(2)Moravec算子只考虑了每隔45度方向，Harris算子用Taylor展开去近似任意方向。将图像窗口平移[u,v]产生灰度变<E(u,v)oE(u,v)二乙w(x,y)[I(x+u,y+v)-1(x,y)]2x,yI(x+u,y+v)=I(x,y)+1u+1v+O(u2,v2)xyE(u,v)=工w(x,y)[Iu+1v+O(u2,v2)]2I2xI2x11xyu11I2vLxyy」x,yu+1v」=[u,v]xy于是对于局部微小的移动量[u,v],可以近似得到下面的表达:E(u,v)=[u,v]M其中M是2X2的矩阵，可由图像的导数求得：M=2w(x,yM=2w(x,y)x,yxxyxyI2式中，I为x方向的差分，/为y方向的差分，w(x，y)为高斯函数。xy(3)Harris采用了一种新的角点判断方法。通过M的两个特征值“»的大小对图像点进行分类。12/•L<Comer,:丸1和鯨都较大且数值相当久]〜-團像窗口在所有方向上移动都严生明显灰度殳化cxEdge?,7如果九和九都很小，图像窗口在所124cTlat;\十―Mg严有方向上移动都无明显灰度变化region\1但是解特征向量需要比较多的计算量，且两个特征值的和等于矩阵M的迹，两个特征值的积等于矩阵M的行列式。所以用下式来判定角点质量。(K常取0.04—0.06)R二detM—k(traceM)2Harris算法总结1：对每一像素点计算相关矩阵MA=w(x,y)012xB=w(x,y)012yC=D=w(x,y)0(II)L-IxyABM=CD2：计算每像素点的Harris角点响应。R=(AB—CD)2—k(A+B)23在wXw范围内寻找极大值点，若Harris角点响应大于阀值，则视为角点。Harris算子对灰度的平移是不变的，因为只有差分，对旋转也有不变性，但是对尺度很敏感，在一个尺度下是角点，在另一个尺度下可能就不是了。二MATLAB代码clear;Image=imread('lOOl.jpg');%读取图像Image=im2uint8(rgb2gray(Image));dx=[-101;-101;-101];%dx:横向Prewitt差分模版Ix2=filter2(dx,Image).人2;Iy2=filter2(dx',Image).A2;Ixy=filter2(dx,Image).*filter2(dx',Image);%生成9*9高斯窗口。窗口越大，探测到的角点越少。h=fspecial('gaussian',9,2);A=filter2(h,Ix2);%用高斯窗口差分Ix2得到AB=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)得到%参数匸点(i,j)八邻域的“相似度”参数，只有中心点与邻域其他八个点的像素值之差在%(-t,+t)之间，才确认它们为相似点，相似点不在候选角点之列t=20;%并没有全部检测图像每个点，而是除去了边界上boundary个像素，%因为我们感兴趣的角点并不出现在边界上boundary=8;fori=boundary:nrow-boundary+lforj=boundary:ncol-boundary+1nlike=0;%相似点个数ifImage(i-1,j-1)>Image(i,j)-t&&Image(i-1,j-1)vImage(i,j)+tnlike=nlike+1;endifImage(i-1,j)>Image(i,j)-t&&Image(i-1,j)vImage(i,j)+tnlike=nlike+1;endifImage(i-1,j+1)>Image(i,j)-t&&Image(i-1,j+1)vImage(i,j)+tnlike=nlike+1;endifImage(i,j-1)>Image(i,j)-t&&Image(i,j-1)vImage(i,j)+tnlike=nlike+1;endifImage(i,j+1)>Image(i,j)-t&&Image(i,j+1)vImage(i,j)+tnlike=nlike+1;endifImage(i+1,j-1)>Image(i,j)-t&&Image(i+1,j-1)vImage(i,j)+tnlike=nlike+1;endifImage(i+1,j)>Image(i,j)-t&&Image(i+1,j)vImage(i,j)+tnlike=nlike+1;endifImage(i+1,j+1)>Image(i,j)-t&&Image(i+1,j+1)vImage(i,j)+tnlike=nlike+1;endifnlike>=2&&nlike<=6Corner(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+1forj=boundary:ncol-boundary+1ifCorner(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;ifCRF(i,j)>CRFmaxCRFmax=CRF(i,j);end;endend;end;%CRFmaxcount=0;%用来记录角点的个数t=0.01;%下面通过一个3*3的窗口来判断当前位置是否为角点fori=boundary:nrow-boundary+1forj=boundary:ncol-boundary+1ifCorner(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加1else%如果当前位置(i,j)不是角点，则在Corner(i,j)中删除对该候选角点的记录Corner(i,j)=0;end;end;end;end;%disp('角点个数')；%disp(count)figure,imshow(Image);%displayIntensityImageholdon;%toc(t1)fori=boundary:nrow-boundary+1forj=boundary:ncol-boundary+1column_ave=O;row_ave=0;k=0;ifCorner(i,j)==1forx=i-3:i+3%7*7邻域fory=j-3:j+3ifCorner(x,y)==1%用算数平均数作为角点坐标，如果改用几何平均数求点的平均坐标，对角点的提取意义不大row_ave=row_ave+x;

column_ave=column_ave+y;k=k+l;endendendendifk>0%周围不止一个角点plot(column_ave/k,row_ave/k,'g.');endend;end;三运行结果PlFigure2|□||BFileEditViewInsertToolsDesktopWindowHelp*□已日釦氐|驚決黔®毘尿”|忌|!3匡1|・啦「1002特征点1001特征点1002特征点（二）（三）Susan算子」算法公式（1）借助图3-1来解释Susan检测的原理，其中图片是白色背景，有一个颜色比较暗淡的矩形。用一个圆形模板在图像上移动，若模板内的像素灰度与模板中心的像素（被称为核Nucleus）灰度值小于一定的阈值，则认为该点与核Nucleus具有相同的灰度，满足该条件的像素组成的区域就称为USAN。在图片上有5个圆形区域。圆形区域表示的是掩码区域。把圆形区域内的每一个位置的像素值与圆心处的像素值相比较，那么圆中的的像素可以分为两类，一类是像素值与圆心处的像素值相近的，另一类是像素值与圆心的处的像素值相差比较大的。~~―——I――r图3-1图3-2如果将模板中各个像素的灰度都与模板中心的核像素的灰度进行比较，那么就会发现总有一部分模板区域和灰度与核像素的灰度相同或相似，这部分区域可以称为USAN。USAN区域包含很多与图像结构有关的信息。利用这种区域的尺寸、重心、二阶矩的分析，可以得到图像中的角点，边缘等信息。从上图所示，当核像素处在图像中的灰度一致区域时，USAN的面积会达到最大。第e个模板就是属于这种情况。（2）Susan进行角点检测时，遵循了常规的思路：使用一个窗口在图像上逐点滑动，在每一个位置上计算出一个角点量，再进行局部极大值抑制，得到最终的角点。我们这里使用的窗口是圆形窗口，最小的窗口是3X3的，此次使用的是37个像素的圆形窗口，如图3-2。（3）在角点检测中，有两种类型的阈值，一种用来约束角点的数量，另一种用来约束角点的质量。当然，一个阈值不能完全做到只影响质量或数量，只是会有一个侧重点。阈值g是角点质量的，尽管也会影响数量，但是相对来说更侧重于影响质量（角点的形状）。例如，g值减小，那么Susan会更加侧重于检测到更加“尖锐”的角点，所以，可以更加自己的实际需求来确定阈值g；而阈值t是角点的数量，当t减小时，会检测到更多的角点，所以,阈值t可以在不影响角点质量的情况下，控制检测到的角点的数量，如果图像的对比度比较低，可以修改t值以适应变化。下面简单叙述下利用Susan算子检测角点的步骤：1：利用圆形模板遍历图像，计算每点处的USAN值。2：设置一阈值g,—般取值为l/2(Max(n)),也即取值为USAN最大值的一半，进行阈值化，得到角点响应。3：使用非极大值抑制来寻找角点。通过上面的方式得到的角点，存在很大伪角点。为了去除伪角点，Susan算子可以由以下方法实现：1:计算USAN区域的重心，然后计算重心和模板中心的距离，如果距离较小则不是正确的角点；2：判断USAN区域的重心和模板中心的连线所经过的像素都是否属于USAN区域的像素，如果属于那么这个模板中心的点就是角点。二MATLAB代码clearall;closeall;clc;img=imread('1001.jpg');img=rgb2gray(img);imshow(img);[mn]=size(img);img=double(img);t=45;%模板中心像素灰度和周围灰度差别的阈值，自己设置usan=[];%当前像素和周围在像素差别在t以下的个数%这里用了37个像素的模板fori=4:m-3%没有在外围扩展图像，最终图像会缩小forj=4:n-3tmp=img(i-3:i+3,j-3:j+3);%先构造7*7的模板，49个像素c=0;forp=1:7forq=1:7if(p-4)A2+(q-4)A2<=12%在其中筛选，最终模板类似一个圆形%usan(k)=usan(k)+exp(-(((img(i,j)-tmp(p,q))/t)人6));ifabs(img(i,j)-tmp(p,q))vt%判断灰度是否相近，t是自己设置的c=c+1;endendendendusan=[usanc];endendg=2*max(usan)/3;%确定角点提取的数量，值比较高时会提取出边缘，自己设置fori=1:length(usan)ifusan(i)<gusan(i)=g-usan(i);elseusan(i)=0;end

endimgn=reshape(usan,[n-6,m-6])';figure;imshow(imgn)%非极大抑制[mn]=size(imgn);re=zeros(m,n);fori=2:m-lforj=2:n-lifimgn(i,j)>max([max(imgn(i-l,j-l:j+l))imgn(i,j-l)imgn(i,j+l)max(imgn(i+l,j-l:j+l))]);re(i,j)=l;elsere(i,j)=0;endendendfigure;imshow(re=l);三运行结果®Figure31002特征点1001特征点（四）三种方法的比较:（1）Moravec