摘 要:在Meyer算法的基礎上,提出一種改進的標記分水嶺遙感影像分割方法,該方法針對高空間分辨率遙感影像的特點,依據梯度影像的分布特征自動提取合適的標記影像。浸沒過程中,非標記像素按照梯度值由小到大進行處理,并使用正反兩個隊列記錄當前處理的像素。實驗證明,將該算法用于高分辨率遙感影像分割,不僅獲得高質量的分割結果,而且具有極高的運行效率與空間利用效率。
關鍵詞:標記分水嶺分割; 遙感影像分割; Butterworth低通濾波
中圖法分類號:TP391
文獻標志碼:A
文章編號:1001-3695(2010)02-0760-04
doi:10.3969/j.issn.1001-3695.2010.02.100
Improved marked-based watershed algorithm for segmentation of high spatial resolution remote sensing imagery
ZHANG Gui-feng, WU Zhao-cong, YI Li-na
(School of Remote Sensing Information Engineering, Wuhan University, Wuhan 430079, China)
Abstract:This paper proposed an improved marked-based watershed image segmentation method based on Meyer’s algorithm. The method extracted the markers automatically according to the gradient distribution. In the flooding process,processed the unmarked pixels sequentially with an ascending order according to the gradient value, and used the forward order priority queue and reverse order priority queue to record the current processed pixel. Experiment results show the proposed algorithm performs well with high efficiency and low computer memory cost.
Key words: marked-based watershed; remote sensing imagery segmentation; Butterworth low-pass filtering
0 引言
面向對象的高分辨率遙感影像處理方法,不僅利用了對象的光譜信息,而且還可以利用目標地物的上下文信息、結構信息和形狀信息[1]。對象的獲取是面向對象的高分辨率遙感影像處理方法的前提條件,一般是通過影像分割來實現。目前主要有三類影像分割方法,即基于全局閾值的分割方法、基于邊緣的分割方法以及基于區域的影像分割方法。
分水嶺分割方法作為一種基于區域的影像分割方法,具有簡單、快速,能檢測出弱邊緣對象及分割后對象具有完整邊界等諸多優點,被廣泛地應用于各種領域。但因為分水嶺分割是在梯度影像上進行的,受暗噪聲等因素的影響存在偽局部最小區域,會出現過分割現象。標記分水嶺影像分割是解決過分割現象的一個好方法,它能抑制偽局部極小值區域。標記分水嶺依據標記影像修改梯度影像,消除了梯度影像上偽局部極小值區域[2~5]。修改梯度影像之后再進行普通的分水嶺分割效率較低,為解決這一問題,Meyer提出了一種使用等級隊列的基于浸沒的標記分水嶺分割方法[6],這一方法是在浸沒過程中抑制偽局部極小值區域。標記分水嶺影像分割中,標記的提取是關鍵的一步,標記影像是否合適直接決定了分割質量。文獻[6~8]中針對不同的影像提出了多種穩健的標記影像自動提取方法。
本文針對高分辨率地物種類多、紋理豐富及具有大量噪聲等特點,提出了一種依據影像上梯度分布特征的標記影像自動提取方法。另外,本文還改進了傳統的Meyer標記分水嶺影像分割方法,浸沒過程中使用了正反兩個隊列記錄當前處理像素,解決了Meyer算法中隊列過多的問題[7],有效地利用內存,提高了影像分割速度。
1 標記分水嶺影像分割
分水嶺影像分割是一種形態學的影像分割工具,它可以分為兩類,即基于浸沒模型的分水嶺分割和基于地形距離的分水嶺分割。因為執行效率高于基于地形距離的分水嶺分割方法,基于浸沒的分水嶺分割方法的應用更廣泛。本文所提的標記分水嶺分割方法也是一種基于浸沒模型的分水嶺分割方法。
1.1 基本概念
在介紹標記分水嶺分割之前,先介紹幾個本文使用到的基本概念。
a)極小值區域[9]。它是影像上的一個連通區域,連通區域中的像素值相等,且區域中每一像素值均不大于其鄰域像素值。
b)梯度影像。它由原始影像與梯度算子卷積計算得到,其尺寸與原始影像一致。本文中為計算方便,每個像素的梯度均由整數表示,且取值控制在0~255。
c)標記[9]。影像上連通區域中像素的集合,這些像素是對象的內部像素,每一個標記對應著一個對象。
d)標記影像[9]。它是一個二值影像,尺寸與原始影像一致,影像中前景像素值為1,是屬于標記的像素,背景像素值為0。影像上每一個連通的前景區域對應著一個標記。
1.2 基于浸沒模型的分水嶺分割
基于浸沒模型的分水嶺分割算法將梯度影像當做一個地形表面,分割時在地形表面的每個局部(局部極小值區域)最低點處刺一個小孔,然后將地形表面垂直按入湖水中,水從洞中進入并淹沒表面。隨著水位的上升,從不同局部最低點進入的水將會匯合,為防止匯合,在這些可能匯合的地方筑建大壩,這些大壩就是分水嶺[10]。
1.3 標記分水嶺分割
因暗噪聲等的影響,梯度影像上存在大量的局部極小值區域,分水嶺分割后存在過分割現象。目前解決過分割現象的方法主要有三種,即影像預處理、分割后處理以及標記分水嶺分割。影像預處理是通過低通濾波等方法,去除偽局部極小值區域;分割后處理主要通過合并相鄰的過分割斑塊來達到目的;而標記分水嶺分割方法則是利用標記影像抑制偽局部極小值區域。本文正是利用標記分水嶺分割方法來減少高分辨率遙感影像分水嶺分割時存在的過分割現象。
標記分水嶺方法主要包括標記影像的自動提取以及像素的標記兩個主要步驟。標記影像的自動提取指的是自動從影像上提取出每個對象內部像素形成標記影像。像素的標記有兩種方法:a)將提取的標記強制作為原梯度影像的局部極小值,同時消除其他局部極小值,然后在修改后的梯度影像上進行普通的分水嶺分割算法。b)直接在原始的梯度影像上進行分割,在分割過程中抑制偽局部極小值區域,Meyer標記分水嶺分割就是這種方法。
1.4 Meyer標記分水嶺影像分割算法
Meyer算法采用具有不同優先級的隊列和先進先出的數據結構來實現,每個像素的優先級由像素取值決定,相同優先級的像素按進入隊列的順序由先到后處理,不同優先級的像素按優先級由高到低處理。其算法流程[7]可概述如下:
首先初始化,初始化包括給標記影像上每個連通區域賦予惟一的標記號以及初始化每個連通區域對應的優先級隊列。每一像素的優先級是由像素梯度值定義,梯度值越大優先級越低,反之則越高。初始化某一連通區域的優先級隊列就是將該區域塊相鄰的非標記像素按該像素的取值壓入到相應的優先隊列中,例如某像素梯度為128,則將其壓入對應梯度為128的優先隊列中。
非標記像素的標記是一個區域增長的過程,主要包括以下幾個步驟:
a)對于每個連通區域(對應著一個對象)取出最高優先級隊列的最后一個元素p,將其標記為該對象的序號。
b)對于p的鄰域像素q,如果其未標記并未壓入優先隊列中,則將其壓入該對象的對應優先級隊列中。如果q的優先級高于當前處理的優先級,則將其壓入當前處理的優先級隊列中。利用此操作可以抑制偽局部極小值區域。
由上述Meyer算法的簡介可知,該算法所需隊列個數與對象個數、優先級個數以及影像大小有關,不適合數據量大、地物復雜的遙感影像分割。
2 改進的標記分水嶺遙感影像分割方法
2.1 影像濾波
考慮到高空間分辨率遙感影像噪聲明顯,采用傳統的低通濾波方法容易破壞對象的輪廓,因此采用同組濾波技術對原始影像進行濾波,確保在有效去除噪聲的同時較好地保留對象的輪廓信息[1]。
同組濾波算法的主要思想是找出濾波窗口中與中心像素特征距離相近的鄰域像素作為同組成員參加濾波。算法首先根據濾波窗口中鄰域像素與中心像素特征距離值進行升序排列,再通過Filsher判別準則找出像素的同組成員,即與中心像素值相近的點,然后用濾波窗口中屬于同組成員的像素加權平均值作為中心像素的特征值[1,11]。
2.2 梯度信息提取
梯度影像更能反映變化趨勢,因此分水嶺分割一般是在梯度影像上進行的。原始影像經同組濾波去噪后,即可用形態梯度運算來求取梯度影像。形態梯度影像等于原始影像膨脹變換后的影像減去腐蝕變換后的結果,其公式[2,3]為
MG=fB-fΘB
其中:MG為梯度影像;、Θ分別表示數學形態學的膨脹和腐蝕運算;而B則為結構元素。結構元素的尺寸對梯度圖像結果有一定影響,本文采用3×3的平板結構元素來求取形態學梯度。彩色影像梯度計算時,首先計算每個通道的梯度,然后將三個通道的梯度平均值作為對應像素的梯度[3]。
2.3 標記的自動提取及標記
考慮到高分辨率遙感影像上同類地物內部像素具有同質性,對象內部梯度值較小,因此使用閾值法分割梯度影像來提取標記影像。但是高分辨率遙感影像上對象內部像素灰度值可能存在變化,同時一些對象具有豐富的紋理,這就導致某些對象內部的梯度值有可能大于其他對象邊界處的梯度??梢娫谠继荻扔跋裆鲜褂脝我婚撝颠M行分割顯然不能獲得合適的標記影像。
鑒于此,本文使用Butterworth低通濾波方法對梯度影像進行低通濾波。Butterworth濾波既能實現有效的低通濾波,又不會出現明顯的振鈴效應,濾波結果能反映每個區域的梯度分布[3,12,13]。用原始梯度影像減去濾波后的影像,所得的影像即可用全局閾值法獲取比較合理的標記影像。因此,定義原始梯度影像為MG,則本文標記影像的提取由以下幾個步驟組成:
a)使用Butterworth濾波器對MG進行濾波,得到濾波后的影像BG。
b)計算MG與BG的差值影像DG=MG-BG。
c)計算DG的中位數T,以中位數為閾值將影像二值化,形成標記影像BW。其中BW=(DG d)影像BW上,每個連通區域對應一個對象。針對某些無意義的連通區域只包含比較少的像素,可預先設定一區域面積閾值(像素個數),刪除包含像素個數小于閾值的區域,即將該區域內的像素值全賦為0。所得的BW即為標記影像,該影像上前景像素為標記像素,其他為未非標記像素。 2.4 像素的標記 像素的標記指的是給像素賦予一個特定的序號,主要包括兩個部分,即標記像素的標記和非標記像素的標記。 本文所提算法與Meyer一樣,也是一個基于浸沒的分水嶺分割算法。與Meyer不同的是,本文算法在標記影像標記后,并不是直接由每個區域的邊緣像素開始向外增長,而是依像素的優先級順序處理每一個像素,然后按順序查找每個非標記像素鄰域中已經被標記的像素,如果存在,則將該鄰域像素對應的對象序號賦給該非標記像素。像素的優先級由梯度進行定義,梯度越小,優先級越高,反之則越低。 2.4.1 非標記像素排序 非標記像素被標記前,需按像素原始梯度MG進行排序,本文算法計算的梯度值在0~255,因此可以利用基于直方圖的方法進行排序。經過排序,非標記像素按梯度值由小到大排列,其中同一梯度級的像素連續排列。因此可以利用直方圖統計出順序排列中每個梯度級第一個像素的位置。若梯度影像包含n個梯度級,則設hist為梯度直方圖,psum為存儲每個梯度級已掃描像素個數的數組,pOrder為存儲排序后像素的坐標數組(像素的行值與列數的乘積再與列值之和),則基于直方圖排序的主要步驟如下: a)統計梯度影像直方圖hist(標記像素不參加統計),hist第i個元素對存儲的是第i級梯度像素的個數。 b)將hist進行累加計算,其計算偽代碼為for( int i = 1; i < n; i++ ) hist[i] += hist[i-1]。 c)從最低梯度級開始,其起始像素位置分別為0,hist[0],hist[1],…hist,[n-2]。 d)初始化每一灰度級已處理好的像素個數psum(i)=0;i=0,1,…,n-1。 e)掃描影像中非種子影像,若其梯度為第i個灰度級,則其順序為hist(i)+psum(i),修改pOrder第hist(i)+psum(i)個元素,存儲該像素的坐標,并將psum第i個元素加1。直至所有非標記像素都掃描到即停止。 2.4.2 標記像素標記 對于標記影像上各區域,采取同一區域中各像素標記號相同、不同區域中的標記號不同的原則對標記影像上各標記像素進行標記。本文使用種子跟蹤算法跟蹤出每一個連通(4-或8-連通),然后給每個連通區域賦予一個惟一的序號,該連通區域內包含的像素也標記為該序號。若定義初始的標記號regionNo=1,所有像素都設置為未處理,然后掃描標記影像,若掃描到的為未處理的標記像素p,則標記其所屬連通區域包括兩步:a)利用種子生長法跟蹤屬于p所在區域的所有像素;b)將跟蹤出來的像素序號全部標記為regionNo,且將這些像素設置為已處理,并將regionNo值加1。依照這兩個步驟處理每一個未處理的標記像素。 2.4.3 標記非標記像素 與Meyer方法相異的是,本文所提算法依據像素梯度的值由小到大進行標記(標記某像素指的是給該像素賦予惟一的對象號)。標記非標記像素時,是通過查找鄰域是否有已經被標記的像素,如果有則將其鄰域像素的標記號賦給該像素。實際標記過程中,因處理有先后次序,在標記某一梯度水平的像素時,有些本應被標記的像素可能未被標記。為解決這一問題,本文使用正序與逆序兩個隊列進行處理。正序隊列記錄所有當前待處理的像素,而逆序隊列則記錄正序隊列中目前不能被標記的像素。 在2.4.1節中已對未標記影像進行了排序,并將每一個非標記像素的影像坐標存儲在數組pOrder中,若以queue1代表正序隊列,queue2代表逆序隊列,則非標記像素的標記主要過程如下: a)設置pOrder中當前處理元素序號pCurNo為0。 b)從pOrder區域中取出第pCurNo個元素并獲得該元素對應像素的梯度,將pOrder中對應梯度為該梯度的元素全壓入隊列queue1中,若共有k個元素是該梯度水平,則將pCurNo值修改為pCurNo+k,此時queue1中像素均為當前待處理像素。 c)彈出queue1隊尾像素p,若p未標記,查找p鄰域是否有已標記的像素,如果有且標記號為code,將p標記為code,同時以p為種子跟蹤出與p在一個區域的所有當前待標記的像素,并將這些像素均標記為code。如果沒有查找到,則將p壓入隊列queue2。重復這一步直到隊列queue1為空。 d)彈出queue2隊尾像素q,若q未標記,查找q鄰域是否有已標記的像素,如果有且標記號為code,將q標記為code,同時以q為種子跟蹤出與q在一個區域的所有當前待標記的像素,并將這些像素均標記為code。如果沒有查找到,則將q壓入隊列queue1。重復這一步直到隊列queue2為空。此時queue1中像素將與梯度級更高的像素一起處理。 (e)重復(b)~(d),直到pOrder中元素全部處理完。 (f)經過(a)~(e)的處理,queue2中還有元素為未處理,因此還需重復運行(c)、(d),直到兩個隊列中元素全為0。 3 實驗分析 利用本文所提的標記影像提取方法可以獲得合適的標記影像,為驗證該標記影像提取方法的有效性,選擇武漢郊區秋季拍攝的大小為400×402的分辨率為0.6m的Quickbird全色影像(圖1)進行標記提取與分割。圖2是圖1的分割結果。該分割結果是使用文獻[3]所提的標記自動提取方法得到的結果。從圖中可以看出,樹林等紋理豐富的地方明顯存在過分割現象,造成這種現象的原因是紋理豐富地區梯度較大,如圖3所示。如果閾值過小,則紋理豐富的對象中可提取出多個由標記像素組成的塊;如果提高閾值,則提取出來的影像中,相鄰的對比度小的對象有可能被合為一個對象,分割結果中就會出現欠分割現象。圖4是利用本文所述標記自動提取方法得到的分割結果,圖中面積閾值與圖2中一樣,都是20,圖4的分割結果明顯優于圖2,而且其具有較好的分割結果。 與Meyer標記分水嶺方法相比,本文所提標記分水嶺分割算法只使用了兩個隊列,而且隊列最大長度已知,有效地提高了內存利用效率并節約了內存的開辟時間。圖4、6、8均是利用本文所提的算法分割不同影像的結果,使用相同標記影像時,Meyer算法分割的結果與本文算法所得結果相同。從圖中可以看出,基于本文所提的標記影像提取算法與標記分水嶺分割算法能夠消除過分割現象,分割獲得的斑塊與實際情況基本相符。 表1是在CPU為P4 2.7 GHz,內存為1 GB上的分割結果對比。從表中可以看出,當影像尺寸較小,包含對象也比較少時,Meyer的執行效率優于本文算法;當影像尺寸較大時,本文算法效率明顯高于Meyer算法。當影像較小時,包含的對象少,程序運行過程中所需優先級隊列個數及內存空間都比較少,算法在內存處理方面花費時間不多,因此Meyer 算法效率比較高;但是當影像較大,包含的對象個數增多時,所需優先級隊列個數呈線性增加,這時不僅需要較大的內存來維持程序的運行,而且此時每個隊列的長度也比較難以確定,運行過程中需要花費大量的時間來調整每個隊列的長度以使隊列不至于溢出[8]。而本文算法在程序一開始就可分配確定長度的隊列,與Meyer算法相比,其在內存管理方面的花銷減少,因此本文算法在處理大圖像時效率優于Meyer算法。當影像尺寸特別大時,Meyer的優先級隊列所需內存巨大,程序很有可能不能獲得充分的內存空間,這樣就會導致不能獲得最終的分割結果,而本文所提算法所需內存空間明顯少于Meyer算法。 4 結束語 本文所提算法可以提取合適的標記影像,分割效果與Meyer算法效果類似。對于中小尺寸影像的分割,Meyer算法執行效率優于本文所提的算法;但分割大尺寸影像時,本文算法耗時小于Meyer算法。實驗證明,本文算法適合大尺寸、內容豐富的遙感影像分割,并獲得較好的分割結果。 參考文獻: [1]巫兆聰,覃茂運. 一種顧及幾何特征的云模型遙感影像分割方法[J]. 武漢大學學報, 2008,33(9):940-942. [2]劉洲峰. 改進的分水嶺分割算法在圖像分割中的應用[J]. 中原工學院學報, 2008,19(1):13-15. [3]陳波,張友靜,陳亮. 標記分水嶺算法及區域合并的遙感影像分割[J].國土資源遙感, 2007, 18(2):35-38. [4]張鯤,王士同.分水嶺算法和基于MRF的層次聚類相結合的混合無監督圖像分割算法[J].計算機應用, 2007,27(3):673-676. [5]孫穎,何國金.基于標記分水嶺算法的高分辨率遙感影像分割方法[J].科學技術與工程, 2008,18(11):2776-2781. [6]ROERDINK J B T M,MEIJSTER A. The watershed transform: definitions, algorithm and parallelization strategies[J]. Fundamenta Informaticae,2001,41(1-2): 187-228. [7]MEYER F. Color image segmentation[C]//Proc of International Conference on Image Processing and Its Applications,1992:303-306. [8]GAO Hai,LIN Wei-si, XUE Ping. Marked-based image segmentation relying on disjoint set union [J]. Signal Processing: Image Communication,2006,21(3):100-102. [9]SERRA J,VINCENT L. An overview of morphological filtering[J]. Circuits Systems Signal Process, 1992, 11(1):555-557. [10]VINCENT L,SOILLE P. Watersheds in digital spaces:an efficient algorithm based on immersion simulations[J].IEEE Trans on Pattern Anal,1991,13(6): 583-589. [11]陳忠.高分辨率遙感影像分類技術研究[D].北京:中國科學院,2006. [12]王啟志.神經網絡辨識的自適應逆控制[J].華僑大學學報, 2005,26(4):397-400. [13]戴金偉,陳增祿,毛惠豐,等. Butterworth帶通濾波器設計[J].西安工程科技學院學報, 2007,21(3):367-370. [14]陳忠,趙忠明,宮鵬.一種快速高分辨率遙感影像分割算法[J].計算機應用研究, 2006,23(10):154-155.