吳琳娜

(上海宇航系統(tǒng)工程研究所,上海 201108)

0 引言

圖像分割是將圖像劃分為有意義區(qū)域的處理技術(shù),它可顯著減少隨后的圖像識別、分析等階段的數(shù)據(jù)處理量。因分割中出現(xiàn)的誤差會傳播到更高層次處理階段,故其精準(zhǔn)度至關(guān)重要。圖像分割技術(shù)在醫(yī)學(xué)診斷及運(yùn)動分析、結(jié)構(gòu)分析、衛(wèi)星遙感圖像分析等領(lǐng)域有極為廣泛的應(yīng)用。工程中常對圖像分割算法有一定的速度要求,并希望減少人工干預(yù),因此圖像分割算法的效率和自動化程度同樣不可忽略。

目前,常用的三維圖像分割算法各有優(yōu)缺點(diǎn)。如分類器算法是一種監(jiān)督性算法,具有迭代特點(diǎn),計算量較小,可用于多通道圖像,但需要大量人工干預(yù),算法自動化難度較大;三維邊緣檢測方法對不同區(qū)域?qū)Ρ榷葟?qiáng)烈的數(shù)據(jù)分割效果很好,但即使檢測出邊緣,卻很難檢測出邊緣與感興趣區(qū)域邊界的關(guān)系;區(qū)域增長算法是一種根據(jù)預(yù)先定義的連接標(biāo)準(zhǔn)從三維圖像數(shù)據(jù)中提取連續(xù)區(qū)域的技術(shù),算法的實(shí)現(xiàn)非常簡單,但需要人工交互的選取種子點(diǎn),對噪聲和局部體效應(yīng)敏感,分割出的區(qū)域在空間上離散[1]。針對在噪聲較高、待分割的區(qū)域的邊界和背景區(qū)域的對比度非常低條件下分割圖像,以快速準(zhǔn)確提取感興趣區(qū)域(ROI)的邊界、形狀和體積等信息,本文對不同分割算法的組合進(jìn)行了研究。

1 數(shù)學(xué)模型

文獻(xiàn)[2]首次提出了運(yùn)用活動輪廓模型進(jìn)行2D圖像分割的方法。Snake模型以曲線演化方式實(shí)現(xiàn)圖像分割。已知三維圖像中所有點(diǎn)的灰度值信息,u0(x,y,z)為圖像每點(diǎn)(x,y,z)的灰度值。令S為ROI的邊界;inside S,outside S分別表示ROI內(nèi)、外部的點(diǎn)。MUMFORD等提出了能量方程：三維圖像在以S為輪廓邊界時能量函數(shù)可表示為

式中：area S為ROI輪廓的表面積;s+,s-分別為ROI內(nèi)、外部所有體素灰度值的均值,s-即為圖像中視作ROI背景所有體素灰度值的均值;υ,λ1,λ2為權(quán)重系數(shù)[3]。式(1)中,等號右端第二部分是對ROI內(nèi)部所有點(diǎn)的積分,是ROI內(nèi)部體素灰度值的方差;第三部分是對ROI外部所有點(diǎn)的積分,是ROI外部體素灰度值的方差。對所分割的圖像有如下的先驗(yàn)知識,先收縮ROI的表面積,并盡可能擬合連通的區(qū)域;其次,最小化ROI內(nèi)部和外部體素灰度值的方差之和,使內(nèi)部和外部的體素灰度值特征盡可能均衡。由此,圖像分割轉(zhuǎn)化成能量方程的最小化

基于以上概念,文獻(xiàn)[4]提出了無邊緣的活動輪廓法,但在解決本文問題時該經(jīng)典算法的效果并不理想。實(shí)驗(yàn)中發(fā)現(xiàn)該算法在一定程度上消除了離散的點(diǎn)和區(qū)域,但同時也將部分ROI外的體素劃進(jìn)ROI內(nèi)部。

為進(jìn)一步利用圖像中除灰度值以外的信息,考慮通過變化提取更多的圖像信息(目前能量方程僅包含單通道圖像的信息)。GABOR在1946年就提出了Gabor方程,之后DAUGMAN將其推廣到了二維[5-6]。本文研究的圖像因采集方式的關(guān)系,圖像ROI區(qū)域具備紋理特征。Gabor方程常用于有一定紋理特征的圖像[7]。Gabor變換可分割空間頻域、元素密度、方向、相位、能量有區(qū)別的區(qū)域[8]。二維Gabor方程或稱Gabor濾波器可表示為

式中：gσ(x,y)=為二維高斯方程;F,θ分別為頻率和相位參數(shù),決定了該變換作用的方向。此處：σ2為方差,其值越大,變化影響的鄰近區(qū)域就越大。u0→(u0)σ,F,θ是空間點(diǎn)u0(x,y,z)在Gabor變換之下的值。該變換的特點(diǎn)是與θ向垂直的邊界會被突出,而對所作用的圖像來說,變換結(jié)果會變得平滑,起到降低噪聲和增強(qiáng)局部紋理特征的作用。每選取一個θ對圖像進(jìn)行變化,就能獲得圖像在一個新通道上的信息。通過N次Gabor變化所得N通道的信息可用能量方程表示為

能，其是滿足顧客入住需求的關(guān)鍵。這個過程中酒店人員利用酒店前臺管理系統(tǒng)對顧客信息進(jìn)行錄入，客人正式入住后邊可通過相應(yīng)系統(tǒng)為客人做其他需求服務(wù)辦理。通常酒店前臺管理系統(tǒng)本身靈活性也相對較高，比如其在權(quán)限操作上，不同管理人員所具有的操作權(quán)限也有一定差異性。結(jié)合實(shí)際來看當(dāng)前酒店前臺管理系統(tǒng)服務(wù)需求主要體現(xiàn)在房間預(yù)訂、登記入住、換房、加房、修改房價等特性，注重前臺管理系統(tǒng)完善性和功能要點(diǎn)的靈活設(shè)置，是整個酒店管理系統(tǒng)得以流暢運(yùn)行的關(guān)鍵[1]。

分別為圖像在第i個通道ROI內(nèi)、外部變換值的均值;s0+,s0-分別為ROI內(nèi)、外部在每個通道上變換值均值的向量,且,S更新時兩組向量值將重新計算;λi,γi為權(quán)重系數(shù),取值源于實(shí)驗(yàn)或經(jīng)驗(yàn)數(shù)據(jù)。當(dāng)υ的取值固定后,對比實(shí)驗(yàn)結(jié)果選取合適的λi,γi值。

為表示輪廓S,本文引入水平集方程。水平集無顯式的定義邊界,僅規(guī)定邊界上的值φ(x,y,z)=0,ROI內(nèi)部φ(x,y,z)＞0,而外部取負(fù)值。最簡單的度量水平集的函數(shù)是Dirac方程,inside S取值為1,outside S取值為0。BAYER研究了多種近似Dirac方程的方法[9]。本文采用其中正則化的Heaviside函數(shù)：H2,ε(z)=,δ2,ε(z)=H′2,ε(z)=擬合Dirac方程的δ函數(shù)表達(dá)。此處：ε為當(dāng)z在0附近時H2,ε(z)的斜度。

為優(yōu)化能量模型并簡化計算,本文引入統(tǒng)計模型。以2D圖像為例,作如下假設(shè)：

a)輪廓(分割邊界)的長度與其成為邊界的概率滿足關(guān)系p(P(Ω))∝e-υ|C|,υ＞0。此處：|C|為輪廓的長度,即輪廓Ω越短,成為圖像邊界的概率越大。

b)無關(guān)分布,最大后驗(yàn)概率(MAP)可表示為

式中：I為待分割圖像的已知信息。

c)ROI內(nèi)、外部的像素灰度值均滿足正態(tài)分布,如圖1所示。兩部分灰度值均值接近,且灰度值的取值范圍有重疊,但各自滿足正態(tài)分布。

對式(4)求負(fù)對數(shù)。令E(Ω)=-lg(p(P(Ω)|I)),Ω1為inside S的區(qū)域,Ω2為outside S,則有

圖1 ROI內(nèi)外灰度取值分布Fig.1 Distribution of gray valueinside/outside ROI

基于這些概率假設(shè),最大化p(P(Ω)|I)和最小化E(Ω)等價?？砂l(fā)現(xiàn)此處E(Ω)也可視作2D圖像能量方程的表示。這些假設(shè)和推導(dǎo)對三維圖像同樣有效。對三維圖像的能量方程,綜合概率統(tǒng)計的模型和水平集的表達(dá)式為

式中：

此處：φ為水平集函數(shù),定義了點(diǎn)(x,y,z)所處的區(qū)域;H(φ(x,y,z))為水平集函數(shù)的正則化表達(dá),取值為(0,1),以免數(shù)值范圍過大;δ(φ(x,y,z))為H(φ(x,y,z))的導(dǎo)數(shù)。

2 數(shù)值計算

2.1 算法性能

三維圖像處理中,須考慮數(shù)值計算量。圖像大小1 024×1 024×1 024時,點(diǎn)的存儲空間可達(dá)GB級別的存儲量,在確保不丟失ROI信息的前提下,應(yīng)采用預(yù)提取或重新采樣法減小圖像尺寸。水平集已廣泛用于描述輪廓和分割物體中,但由于其計算的復(fù)雜度與圖像的體素成正比,而普通的水平集法是用尺寸固定的網(wǎng)格點(diǎn)計算未知點(diǎn),用偏微分方程推進(jìn)水平集的移動,導(dǎo)致迭代周期過長。

ADALSTEINSSON提出了窄帶技術(shù)以加快水平集方程的移動速度,大幅降低額外的計算開銷,進(jìn)而提高算法的效率[10]。方法要點(diǎn)是考慮將輪廓在一定距離范圍內(nèi)的體素點(diǎn)看作一個帶狀,如圖2所示。在每步迭代中,窄帶會由最新的輪廓重新計算而得。SETHIAN提出了快速步進(jìn)法計算距離[11]。該法的要點(diǎn)是在零水平集中建立可適應(yīng)的網(wǎng)格。

圖2僅考慮零水平集的在管道內(nèi)部鄰近點(diǎn),計算復(fù)雜度從O(n3)降至O(n2)。此處：n為圖像在一個維度上的點(diǎn)數(shù)。同樣,對內(nèi)存的需求也由三階降為二階。窄帶方法的缺點(diǎn)是,管道永遠(yuǎn)不能覆蓋整個圖像區(qū)域。對管道,初始化方法有多種,最簡單的是將圖像中心作為球心,以指定長度為半徑畫出球形作為初始管道,如半值過大會導(dǎo)致計算量增加,半徑過小又會導(dǎo)致管道與ROI無交集,從而混淆ROI的內(nèi)部和外部。另外引進(jìn)重新計算距離技術(shù),可避免水平集方程演進(jìn)過程中因過于平坦而粘滯在局部區(qū)域。

圖2 不使用和使用窄帶的水平集Fig.2 Level Set without/with narrow band

2.2 收斂條件

通常,偏微分方程求解采用標(biāo)準(zhǔn)的半隱式方法和AOS方法。AOS特點(diǎn)是：適于并行計算,可直接推導(dǎo)到更高階的問題無需修改算法的結(jié)構(gòu)。WEICKERT等證明了AOS法效率較顯式方法高10倍[12]。對實(shí)時性能要求,必須為算法收斂性和收斂速度考慮合理的收斂條件。CHEN-VESE給出的收斂條件是在管道內(nèi),點(diǎn)的水平集的變化絕對值均值小于某收斂值,即輪廓在迭代中“幾乎”不再移動時算法收斂,此時能量方程也趨于穩(wěn)定。同時還會設(shè)定一個固定的最大迭代次數(shù)值以控制迭代次數(shù)。

2.3 計算步驟

求解水平集的數(shù)值實(shí)現(xiàn)是基于對偏微分方程的有限差分逼近。通過改變水平集方程來驅(qū)動輪廓的運(yùn)動。其步驟如下：

步驟a)由于窄帶方法假設(shè)管道邊界上的網(wǎng)格點(diǎn)和管道鄰近點(diǎn)有相同的符號,每次迭代前邊界上的值均需重新初始化。

步驟b)每次迭代中,重新計算每個通道的水平集統(tǒng)計信息。

步驟c)用新的水平集統(tǒng)計信息值更新水平集的力場值,對應(yīng)偏微分方程用向前歐拉步進(jìn)法更新。平均曲率的運(yùn)動用加性分裂算子(AOS)逼近。

步驟d)重新計算水平集距離,更新窄帶。

步驟e)在下次迭代開始前,根據(jù)收斂條件判斷循環(huán)是否停止。

3 實(shí)驗(yàn)結(jié)果及分析

實(shí)驗(yàn)取4個方向作為Gabor變換的方向,均位于與Z軸垂直的平面中,θ的取值為0,,即N=4。設(shè)4個通道Gabor變換的權(quán)重參數(shù)為λi==200,i=0,1,2,3。通過對一組數(shù)據(jù)進(jìn)行20組實(shí)驗(yàn),選取不同的λ值,對比實(shí)驗(yàn)結(jié)果的誤差,發(fā)現(xiàn)當(dāng)λ=800時,分割誤差值最小。在15組三維圖像的數(shù)據(jù)的分割運(yùn)算中,平均在20次迭代之內(nèi),算法就能達(dá)到收斂,平均運(yùn)算時間小于5 s,應(yīng)用中能實(shí)現(xiàn)快速分割。

用戶用工具手動畫出的輪廓如圖3(a)所示,本文算法提取的圖像如圖3(b)所示。用HUTTENLOCHER的Hausdorff距離判別兩條輪廓曲線的誤差。其中,15組數(shù)據(jù)的分割圖像結(jié)果的Hausdorff距離平均值為0.966 5,表明算法提取的輪廓與驗(yàn)證結(jié)果的像素偏離很小,準(zhǔn)確度達(dá)到期望結(jié)果。

圖3 三維圖像分割結(jié)果的一個二維截面Fig.3 2D section from 3D segmentation results

4 結(jié)束語

本文根據(jù)常用三維圖像分割算法,提出了綜合圖像紋理和統(tǒng)計學(xué)信息的活動輪廓算法,并進(jìn)行了實(shí)驗(yàn)驗(yàn)證。算法的重點(diǎn)是紋理提取變換在多通道中的合成,可提供較原始圖像更完整的分割信息;水平集方法克服了參數(shù)化表示輪廓的缺陷,半隱式的數(shù)值算法可顯著加大每次迭代的步長,加快收斂速度,收斂準(zhǔn)則在迭代趨于穩(wěn)定時輪廓移動較小,可減少總迭代次數(shù)。在三維活動輪廓模型采用了Gabor變換,算法在效率和準(zhǔn)確度均達(dá)到實(shí)用要求。算法的缺陷是不同采集方式獲得的圖像的λ取值并不確定,故還需研究不同類型圖像的不同收斂性和精確度的λ值確定方法。后續(xù)研究還可包括知識庫建立,以及典型圖像特征求解及系列參數(shù)值存貯。

[1]LI S Z,Markov random field modeling in computer vision[M].Tokyo：Springer-Verlag,1995.

[2]MUMFORD D,SHAH J,Optimal approximations by piecewise smooth functions and associated variational problems[R].Center for Intelligent Control Systems,1988.

[3]CHAN T F,VESE L A.Active contours without edges[J].Image Processing,IEEE Transactions on,2001,10(2)：266-277.

[4]DAUGMAN JG.Uncertainty relation for resolution in space,spatial f requency,and orientation optimized by two-dimensional visual cortical filters[J].Optical Society of America,Journal,A：Optics and Image Science,1985(2)：1160-1169.

[5]SANDBERG B,CHAN T,VESE L.A level-set and Gabor-based active contour algorithm or segmenting textured images[R].UCLA Department of Mathematics CAM Report,2002：02-39.

[6]CLARK M,BOVIK A C,GEISLER W S.Texture segmentation using Gabor modulation/demodulation[J].Pattern Recognition Letters,1987,6(4)：261-267.

[7]BEYER R P,LEVEQUE R J.Analysis of a one-dimentional model for the immersed boundary method[J].SIAM J Number Anal,1992,29(2)：332-364.

[8]ADALSTEINSSON D.A fast level set method for propagating interfaces[D].Oakland：University of California,1994.

[9]SETHIAN J A.Level set methods and fast marching methods[J].Computer Vision and Materials Science,Cambridge University Press,1999.

[10]WEICKERT J,ROMENY B,VIERGEVER M A.Efficient and reliable schemes for nonlinear diffusion filtering[J].IEEE Transactions on Image Processing,1998,7(3)：398-410.

[11]HUTTENLOCHER D P,KLANDERMAN G A,RUCK LIDGE W J.Comparing images using the Hausdorff distance[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,1993,15(9)：850-863.

[12]CHEN Y,TAGARE H D,THIRUVENKADAM S,et al.Using prior shapes in geometric active contours in a variational framework[J].International Journal of Computer Vision,2002,50(3)：315-328.