遙感提取特征點

1、遙感影像特征點提取一、 基于Moravec算子的特征點提取1. Moravec算子的原理及算法公式該算子是通過逐像元量測與其鄰元的灰度差，搜索相鄰像元之間具有高反差的點，具體方法有以下幾種。（1） 計算各像元的有利值，如圖所示，在5×5的窗口內(nèi)沿著圖示四個方向分別計算相鄰像元間灰度差之平方和V1，V2，V3，及V4，取其中最小值作為該像元的有利值：其中：式中， 代表像元的灰度值，為以像元計的窗口大小，如圖所示，為像元在整塊影像中位置序號。（2） 給定一個閾值，確定待定點的有利點。如果有利值大于給定的閾值，則將該像元作為候選點。閾值一般為經(jīng)驗值。（3） 抑制局部非最大。在一定大小窗口內(nèi)

2、（例如5×5，7×7，,9×9像元等），將上一步所選的候選點與其周圍的候選點比較，若該像元的有利非窗口中最大值，則去掉；否則，該像元被確定為特征點，這一步的目的在于避免紋理豐富的區(qū)域產(chǎn)生束點，用于抑制局部非最大的窗口大小取決于所需的有利點密度。綜上所述，Moravec算子是在四個主要方向上選擇具有最大最小灰度方差的點作為特征點。2. 基于MATLAB的算法編程clear all;close all;clcimg=double(imread('1001.jpg');h w=size(img);imshow(img,)imgn=zeros(h,w);n

3、=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 %垂直，水平，對角，反對角四個方向領(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)2; end pix=min(V); %四個方向中選最小值 imgn(y

4、,x)=pix; endendT=mean(imgn(:); %設(shè)閾值，小于均值置零ind=find(imgn<T);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 endendfigure;imshow(img,);3. 運行結(jié)果 1001特征點 1002特征點二、Harris角點檢測算子1、算法公式（1） Harris算子用高斯函數(shù)代替二值窗口

5、函數(shù)，對離中心點越近的像素賦予越大的權(quán)重，以減少噪聲影響。（2） Moravec算子只考慮了每隔45度方向，Harris算子用Taylor展開去近似任意方向。將圖像窗口平移u,v產(chǎn)生灰度變化E(u,v)。由得于是對于局部微小的移動量u,v，可以近似得到下面的表達：其中M是2×2的矩陣，可由圖像的導數(shù)求得：式中， 為x方向的差分， 為y方向的差分， 為高斯函數(shù)。如果和都很小，圖像窗口在所有方向上移動都無明顯灰度變化（3） Harris采用了一種新的角點判斷方法。通過M的兩個特征值的大小對圖像點進行分類。但是解特征向量需要比較多的計算量，且兩個特征值的和等于矩陣M的跡，兩個特征值的積等于

6、矩陣M的行列式。所以用下式來判定角點質(zhì)量。（K常取0.040.06）（4）Harris算法總結(jié)1：對每一像素點計算相關(guān)矩陣M2：計算每像素點的Harris角點響應(yīng)。3：在w×w范圍內(nèi)尋找極大值點，若Harris角點響應(yīng)大于閥值，則視為角點。Harris算子對灰度的平移是不變的，因為只有差分，對旋轉(zhuǎn)也有不變性，但是對尺度很敏感，在一個尺度下是角點，在另一個尺度下可能就不是了。二 MATLAB代碼clear;Image = imread('1001.jpg'); % 讀取圖像Image = im2uint8(rgb2gray(Image);dx = -1 0 1;-1 0

7、 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

8、 = size(Image,2); Corner = zeros(nrow,ncol); %矩陣Corner用來保存候選角點位置,初值全零，值為1的點是角點%真正的角點在137和138行由(row_ave,column_ave)得到%參數(shù)t:點(i,j)八鄰域的“相似度”參數(shù)，只有中心點與鄰域其他八個點的像素值之差在%（-t,+t）之間，才確認它們?yōu)橄嗨泣c，相似點不在候選角點之列t=20;%并沒有全部檢測圖像每個點，而是除去了邊界上boundary個像素，%因為我們感興趣的角點并不出現(xiàn)在邊界上boundary=8;for i=boundary:nrow-boundary+1 for j=bou

9、ndary:ncol-boundary+1 nlike=0; %相似點個數(shù) if Image(i-1,j-1)>Image(i,j)-t && Image(i-1,j-1)<Image(i,j)+t nlike=nlike+1; end if Image(i-1,j)>Image(i,j)-t && Image(i-1,j)<Image(i,j)+t nlike=nlike+1; end if Image(i-1,j+1)>Image(i,j)-t && Image(i-1,j+1)<Image(i,j)+t

10、nlike=nlike+1; end if Image(i,j-1)>Image(i,j)-t && Image(i,j-1)<Image(i,j)+t nlike=nlike+1; end if Image(i,j+1)>Image(i,j)-t && Image(i,j+1)<Image(i,j)+t nlike=nlike+1; end if Image(i+1,j-1)>Image(i,j)-t && Image(i+1,j-1)<Image(i,j)+t nlike=nlike+1; end if I

11、mage(i+1,j)>Image(i,j)-t && Image(i+1,j)<Image(i,j)+t nlike=nlike+1; end if Image(i+1,j+1)>Image(i,j)-t && Image(i+1,j+1)<Image(i,j)+t nlike=nlike+1; end if nlike>=2 && nlike<=6 Corner(i,j)=1;%如果周圍有0，1，7，8個相似與中心的（i,j） %那(i,j)就不是角點，所以，直接忽略 end; end;end;CRF =

12、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)驗值for i = boundary:nrow-boundary+1 for j = boundary:ncol-boundary+1 if Corner(i,j)=1 %只關(guān)注候選點 M = A(i,j) C(i,j); C(i,j) B(i,

13、j); CRF(i,j) = det(M)-t*(trace(M)2; if CRF(i,j) > CRFmax CRFmax = CRF(i,j); end; endend; end; %CRFmaxcount = 0; % 用來記錄角點的個數(shù)t=0.01; % 下面通過一個3*3的窗口來判斷當前位置是否為角點for i = boundary:nrow-boundary+1 for j = boundary:ncol-boundary+1 if Corner(i,j)=1 %只關(guān)注候選點的八鄰域 if CRF(i,j) > t*CRFmax && CRF(i,j)

14、 >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;%這個

15、是角點，count加1 else % 如果當前位置（i,j）不是角點，則在Corner(i,j)中刪除對該候選角點的記錄 Corner(i,j) = 0; end; end; end; end; % disp('角點個數(shù)');% disp(count)figure,imshow(Image); % display Intensity Imagehold 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

16、(i,j)=1 for x=i-3:i+3 %7*7鄰域 for y=j-3:j+3 if Corner(x,y)=1% 用算數(shù)平均數(shù)作為角點坐標，如果改用幾何平均數(shù)求點的平均坐標，對角點的提取意義不大 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.'); endend; end; 三 運行結(jié)果 1001特征點 1002特征點 （二）（三）Susan算子一 算法公式（1）借助圖3-1來解

17、釋Susan檢測的原理，其中圖片是白色背景，有一個顏色比較暗淡的矩形。用一個圓形模板在圖像上移動，若模板內(nèi)的像素灰度與模板中心的像素（被稱為核Nucleus）灰度值小于一定的閾值，則認為該點與核Nucleus具有相同的灰度，滿足該條件的像素組成的區(qū)域就稱為USAN。在圖片上有5個圓形區(qū)域。圓形區(qū)域表示的是掩碼區(qū)域。把圓形區(qū)域內(nèi)的每一個位置的像素值與圓心處的像素值相比較，那么圓中的的像素可以分為兩類，一類是像素值與圓心處的像素值相近的，另一類是像素值與圓心的處的像素值相差比較大的。 圖3-1 圖3-2如果將模板中各個像素的灰度都與模板中心的核像素的灰度進行比較，那么就會發(fā)現(xiàn)總有一部分模板區(qū)域和灰

18、度與核像素的灰度相同或相似，這部分區(qū)域可以稱為USAN。USAN區(qū)域包含很多與圖像結(jié)構(gòu)有關(guān)的信息。利用這種區(qū)域的尺寸、重心、二階矩的分析，可以得到圖像中的角點，邊緣等信息。從上圖所示，當核像素處在圖像中的灰度一致區(qū)域時，USAN的面積會達到最大。第e個模板就是屬于這種情況。（2）Susan進行角點檢測時，遵循了常規(guī)的思路：使用一個窗口在圖像上逐點滑動，在每一個位置上計算出一個角點量，再進行局部極大值抑制，得到最終的角點。我們這里使用的窗口是圓形窗口，最小的窗口是3×3的，此次使用的是37個像素的圓形窗口，如圖3-2。（3）在角點檢測中，有兩種類型的閾值，一種用來約束角點的數(shù)量，另一種

19、用來約束角點的質(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，一般取值為1/2(Max(n)），

20、也即取值為USAN最大值的一半，進行閾值化，得到角點響應(yīng)。3：使用非極大值抑制來尋找角點。通過上面的方式得到的角點，存在很大偽角點。為了去除偽角點，Susan算子可以由以下方法實現(xiàn)：1:計算USAN區(qū)域的重心，然后計算重心和模板中心的距離，如果距離較小則不是正確的角點；2：判斷USAN區(qū)域的重心和模板中心的連線所經(jīng)過的像素都是否屬于USAN區(qū)域的像素，如果屬于那么這個模板中心的點就是角點。二 MATLAB代碼clear all;close all;clc;img=imread('1001.jpg');img=rgb2gray(img);imshow(img);m n=size(

21、img);img=double(img);t=45; %模板中心像素灰度和周圍灰度差別的閾值，自己設(shè)置usan=; %當前像素和周圍在像素差別在t以下的個數(shù)%這里用了37個像素的模板for i=4:m-3 %沒有在外圍擴展圖像，最終圖像會縮小 for j=4:n-3 tmp=img(i-3:i+3,j-3:j+3); %先構(gòu)造7*7的模板，49個像素 c=0; for p=1:7 for q=1:7 if (p-4)2+(q-4)2<=12 %在其中篩選，最終模板類似一個圓形 % usan(k)=usan(k)+exp(-(img(i,j)-tmp(p,q)/t)6); if abs(i

22、mg(i,j)-tmp(p,q)<t %判斷灰度是否相近，t是自己設(shè)置的 c=c+1; end end end end usan=usan c; endendg=2*max(usan)/3; %確定角點提取的數(shù)量，值比較高時會提取出邊緣，自己設(shè)置for i=1:length(usan) if usan(i)<g usan(i)=g-usan(i); else usan(i)=0; endendimgn=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) imgn(i,j-1) imgn(i,j+1) max(imgn(i+1,j-1:j+1); re(i,j)=1; else re(i,j)=0; end endendfigure;imshow(re=1);三 運行結(jié)果 1001特征點 1002特征點（四）三種方法的比較：（1）Morave