遙感提取特征點

遙感影像特征點提取基于Moravec算子的特征點提取1.Moravec算子的原理及算法公式該算子是通過逐像元量測與其鄰元的灰度差，搜索相鄰像元之間具有高反差的點，具體方法有以下幾種。計算各像元的有利值，如圖所示，在5X5的窗口內(nèi)沿著圖示四個方向分別計算相鄰像元間灰度差之平方和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元在整塊影像中位置序號。給定一個閾值，確定待定點的有利點。如果有利值大于給定的閾值，則將該像元作為候選點。閾值一般為經(jīng)驗值。抑制局部非最大。在一定大小窗口內(nèi)(例如5X5，7X7，,9X9像元等)，將上一步所選的候選點與其周圍的候選點比較，若該像元的有利非窗口中最大值，則去掉；否則，該像元被確定為特征點，這一步的目的在于避免紋理豐富的區(qū)域產(chǎn)生束點，用于抑制局部非最大的窗口大小取決于所需的有利點密度。綜上所述，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%垂直，水平，對角，反對角四個方向領(lǐng)域灰度差的平方和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%四個方向中選最小值%四個方向中選最小值%設(shè)閾值，小于均值置零%選局部最大且非零值作為特征點ifmax(sq(:))==imgn(y,x)&&imgn(y,x)?=0img(y,x)=255;endendendfigure;imshow(img,[]);<f■*ij.<T3?運行結(jié)果<f■*ij.<T二、Harris角點檢測算子1、算法公式(1)Harris算子用高斯函數(shù)代替二值窗口函數(shù)，對離中心點越近的像素賦予越大的權(quán)重，以減少噪聲影響。w(x,y)二(2)Moravec算子只考慮了每隔45度方向，Harris算子用Taylor展開去近似任意方向。將圖像窗口平移[u,v]產(chǎn)生灰度變<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的矩陣，可由圖像的導(dǎo)數(shù)求得：M=2w(x,yM=2w(x,y)x,yxxyxyI2式中，I為x方向的差分，/為y方向的差分，w(x，y)為高斯函數(shù)。xy(3)Harris采用了一種新的角點判斷方法。通過M的兩個特征值“?的大小對圖像點進行分類。12/?L<Comer,:丸1和鯨都較大且數(shù)值相當久]?-團像窗口在所有方向上移動都嚴生明顯灰度殳化cxEdge?,7如果九和九都很小，圖像窗口在所124cTlat;\十―Mg嚴有方向上移動都無明顯灰度變化region\1但是解特征向量需要比較多的計算量，且兩個特征值的和等于矩陣M的跡，兩個特征值的積等于矩陣M的行列式。所以用下式來判定角點質(zhì)量。(K常取0.04—0.06)R二detM—k(traceM)2Harris算法總結(jié)1：對每一像素點計算相關(guān)矩陣MA=w(x,y)012xB=w(x,y)012yC=D=w(x,y)0(II)L-IxyABM=CD2：計算每像素點的Harris角點響應(yīng)。R=(AB—CD)2—k(A+B)23在wXw范圍內(nèi)尋找極大值點，若Harris角點響應(yīng)大于閥值，則視為角點。Harris算子對灰度的平移是不變的，因為只有差分，對旋轉(zhuǎn)也有不變性，但是對尺度很敏感，在一個尺度下是角點，在另一個尺度下可能就不是了。二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)得到%參數(shù)匸點(i,j)八鄰域的“相似度”參數(shù)，只有中心點與鄰域其他八個點的像素值之差在%(-t,+t)之間，才確認它們?yōu)橄嗨泣c，相似點不在候選角點之列t=20;%并沒有全部檢測圖像每個點，而是除去了邊界上boundary個像素，%因為我們感興趣的角點并不出現(xiàn)在邊界上boundary=8;fori=boundary:nrow-boundary+lforj=boundary:ncol-boundary+1nlike=0;%相似點個數(shù)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用來保存角點響應(yīng)函數(shù)值，初值全零CRFmax=0;%圖像中角點響應(yīng)函數(shù)的最大值，作閾值之用t=0.05;%計算CRF%工程上常用CRF(i,j)=det(M)/trace(M)計算CRF，那么此時應(yīng)該將下面第105行的%比例系數(shù)t設(shè)置大一些，t=0.1對采集的這幾幅圖像來說是一個比較合理的經(jīng)驗值fori=boundary:nrow-boundary+1forj=boundary:ncol-boundary+1ifCorner(i,j)==1%只關(guān)注候選點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;%用來記錄角點的個數(shù)t=0.01;%下面通過一個3*3的窗口來判斷當前位置是否為角點fori=boundary:nrow-boundary+1forj=boundary:ncol-boundary+1ifCorner(i,j)==1%只關(guān)注候選點的八鄰域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('角點個數(shù)')；%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%用算數(shù)平均數(shù)作為角點坐標，如果改用幾何平均數(shù)求點的平均坐標，對角點的提取意義不大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;三運行結(jié)果PlFigure2|□||BFileEditViewInsertToolsDesktopWindowHelp*□已日釦氐|驚決黔?毘尿”|忌|!3匡1|?啦「1002特征點1001特征點1002特征點（二）（三）Susan算子」算法公式（1）借助圖3-1來解釋Susan檢測的原理，其中圖片是白色背景，有一個顏色比較暗淡的矩形。用一個圓形模板在圖像上移動，若模板內(nèi)的像素灰度與模板中心的像素（被稱為核Nucleus）灰度值小于一定的閾值，則認為該點與核Nucleus具有相同的灰度，滿足該條件的像素組成的區(qū)域就稱為USAN。在圖片上有5個圓形區(qū)域。圓形區(qū)域表示的是掩碼區(qū)域。把圓形區(qū)域內(nèi)的每一個位置的像素值與圓心處的像素值相比較，那么圓中的的像素可以分為兩類，一類是像素值與圓心處的像素值相近的，另一類是像素值與圓心的處的像素值相差比較大的。~~―——I――r圖3-1圖3-2如果將模板中各個像素的灰度都與模板中心的核像素的灰度進行比較，那么就會發(fā)現(xiàn)總有一部分模板區(qū)域和灰度與核像素的灰度相同或相似，這部分區(qū)域可以稱為USAN。USAN區(qū)域包含很多與圖像結(jié)構(gòu)有關(guān)的信息。利用這種區(qū)域的尺寸、重心、二階矩的分析，可以得到圖像中的角點，邊緣等信息。從上圖所示，當核像素處在圖像中的灰度一致區(qū)域時，USAN的面積會達到最大。第e個模板就是屬于這種情況。（2）Susan進行角點檢測時，遵循了常規(guī)的思路：使用一個窗口在圖像上逐點滑動，在每一個位置上計算出一個角點量，再進行局部極大值抑制，得到最終的角點。我們這里使用的窗口是圓形窗口，最小的窗口是3X3的，此次使用的是37個像素的圓形窗口，如圖3-2。（3）在角點檢測中，有兩種類型的閾值，一種用來約束角點的數(shù)量，另一種用來約束角點的質(zhì)量。當然，一個閾值不能完全做到只影響質(zhì)量或數(shù)量，只是會有一個側(cè)重點。閾值g是角點質(zhì)量的，盡管也會影響數(shù)量，但是相對來說更側(cè)重于影響質(zhì)量（角點的形狀）。例如，g值減小，那么Susan會更加側(cè)重于檢測到更加“尖銳”的角點，所以，可以更加自己的實際需求來確定閾值g；而閾值t是角點的數(shù)量，當t減小時，會檢測到更多的角點，所以,閾值t可以在不影響角點質(zhì)量的情況下，控制檢測到的角點的數(shù)量，如果圖像的對比度比較低，可以修改t值以適應(yīng)變化。下面簡單敘述下利用Susan算子檢測角點的步驟：1：利用圓形模板遍歷圖像，計算每點處的USAN值。2：設(shè)置一閾值g,—般取值為l/2(Max(n)),也即取值為USAN最大值的一半，進行閾值化，得到角點響應(yīng)。3：使用非極大值抑制來尋找角點。通過上面的方式得到的角點，存在很大偽角點。為了去除偽角點，Susan算子可以由以下方法實現(xiàn)：1:計算USAN區(qū)域的重心，然后計算重心和模板中心的距離，如果距離較小則不是正確的角點；2：判斷USAN區(qū)域的重心和模板中心的連線所經(jīng)過的像素都是否屬于USAN區(qū)域的像素，如果屬于那么這個模板中心的點就是角點。二MATLAB代碼clearall;closeall;clc;img=imread('1001.jpg');img=rgb2gray(img);imshow(img);[mn]=size(img);img=double(img);t=45;%模板中心像素灰度和周圍灰度差別的閾值，自己設(shè)置usan=[];%當前像素和周圍在像素差別在t以下的個數(shù)%這里用了37個像素的模板fori=4:m-3%沒有在外圍擴展圖像，最終圖像會縮小forj=4:n-3tmp=img(i-3:i+3,j-3:j+3);%先構(gòu)造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是自己設(shè)置的c=c+1;endendendendusan=[usanc];endendg=2*max(usan)/3;%確定角點提取的數(shù)量，值比較高時會提取出邊緣，自己設(shè)置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);三運行結(jié)果?Figure31002特征點1001特征點（四）三種方法的比較:（1）Moravec