基于分水嶺與圖割的自動(dòng)分割方法

劉 榮

(北京航空航天大學(xué) 機(jī)械工程及自動(dòng)化學(xué)院，北京100191)

彭艷敏

(天津醫(yī)科大學(xué) 醫(yī)學(xué)影像學(xué)院，天津300203)

唐 粲 程 勝

(昆山市工業(yè)技術(shù)研究院，昆山215300)

Liu Rong

(School of Mechanical Engineering and Automation，Beijing University of Aeronautics and Astronautics，Beijing 100191，China)

Peng Yanmin

(School of Medical Imaging，Tianjin Medical University，Tianjin 300203，China)

Tang Can Cheng Sheng

(Kunshan Industrial Technology Research Institute，Kunshan 215300，China)

隨著計(jì)算機(jī)可視化信息的發(fā)展，機(jī)器人輔助視覺定位手術(shù)越來越多的應(yīng)用到臨床上.圖像分割不僅是醫(yī)生診斷的基礎(chǔ)，還是視覺定位的關(guān)鍵.醫(yī)學(xué)圖像由于人體復(fù)雜的特殊性，很多通用的圖像分割算法對(duì)醫(yī)學(xué)圖像并不適合，都有一定的局限性.區(qū)域增長(zhǎng)［1］、分裂合并［2］等方法都是基于區(qū)域的信息，雖然簡(jiǎn)單，但提取出物體的邊界模糊.基于活動(dòng)輪廓［3］的算法如snakes、變形和最短路徑法是基于邊界的，很難擴(kuò)展到三維空間上.文獻(xiàn)［4］提出基于邊界和區(qū)域的算法，即圖割，能夠很好的提取出物體的邊界，在計(jì)算機(jī)視覺方面發(fā)揮越來越重要的作用.但該文獻(xiàn)在尋找最小能量的過程中是基于像素的，存在海量計(jì)算問題.文獻(xiàn)［5］提出的分水嶺算法一般都存在過分割現(xiàn)象，提取出的物體沒有太多實(shí)際意義，通常還需要區(qū)域合并.但分水嶺算法能夠很好的定位對(duì)象的邊緣，并且能夠保持每個(gè)小區(qū)域中的微小差異［6］.

本文對(duì)內(nèi)外輪廓之間的圖像進(jìn)行預(yù)分割后，對(duì)分割后的小區(qū)域利用圖割進(jìn)行分割，由對(duì)每個(gè)像素操作變成對(duì)小區(qū)域操作;通過實(shí)驗(yàn)驗(yàn)證本文算法在運(yùn)行時(shí)間和運(yùn)行效果上都優(yōu)于文獻(xiàn)［7］.

1 算法基礎(chǔ)

1.1 圖割的相關(guān)理論

圖像分割可以轉(zhuǎn)換成求解能量函數(shù)最小值問題.圖割是利用圖論中的最小割算法求能量函數(shù)全局最優(yōu)解的算法.利用圖割來解決問題需要做到以下3點(diǎn):

1)構(gòu)造圖結(jié)構(gòu).把一幅圖像轉(zhuǎn)變成有向圖G或者無向圖，圖的頂點(diǎn)為圖像中的每個(gè)像素點(diǎn)，兩個(gè)相鄰像素在圖中表示為一條邊，像素之間的亮度、梯度等信息可以轉(zhuǎn)變成邊的容量.

2)構(gòu)造能量函數(shù).圖割綜合了邊界和區(qū)域的信息，圖像分割問題可以轉(zhuǎn)換為求取使能量函數(shù)取到極小值時(shí)的標(biāo)號(hào)場(chǎng).該能量函數(shù)一般形式可以表示為

其中，Vp，q(fp，fq)表示平滑項(xiàng)能量，是對(duì)不連續(xù)性的一種懲罰，為相鄰像素之間能量;Dp(fp)是數(shù)據(jù)項(xiàng)能量，是對(duì)當(dāng)前標(biāo)號(hào)與觀測(cè)數(shù)據(jù)間不符合程度的一種懲罰.

3)利用最大流/最小割求解能量函數(shù)全局最小值.得到的能量函數(shù)的最小值就是對(duì)圖像的精確分割結(jié)果.在1.3節(jié)詳細(xì)介紹本文采用的最大流/最小割算法.

1.2 多源點(diǎn)多匯點(diǎn)簡(jiǎn)化圖

一幅圖像可看作有向圖 G=〈V，ε〉，其中V是頂點(diǎn)集，代表圖像中的每個(gè)像素，ε是相鄰頂點(diǎn)之間邊的集合，即鄰接邊的集合.對(duì)多源點(diǎn)到多匯點(diǎn)的有向圖，可以采用文獻(xiàn)［7］中把多個(gè)節(jié)點(diǎn)簡(jiǎn)化成一個(gè)節(jié)點(diǎn)，匯合邊的容量，刪除多余邊簡(jiǎn)化圖結(jié)構(gòu)的方式，即把一系列的源點(diǎn) Vs={s1，s2，s3，…，sn}簡(jiǎn)化成一個(gè)新的節(jié)點(diǎn)V1，把一系列的匯點(diǎn)Vt={t1，t2，…，tn}簡(jiǎn)化成一個(gè)新的節(jié)點(diǎn) V2，對(duì)于每一個(gè)與源點(diǎn)相連的節(jié)點(diǎn)，其邊的容量為所有與源點(diǎn)相連邊的容量之和，同理與匯點(diǎn)相連的節(jié)點(diǎn)的容量，為這個(gè)節(jié)點(diǎn)與所有匯點(diǎn)相連邊的容量之和.這樣有向圖G簡(jiǎn)化成一個(gè)新的有向圖G′=〈V′，ε′〉，其中 V′=V －Vs－ Vt+V1+V2.對(duì)有向圖G=〈V，ε〉和它的簡(jiǎn)化圖 G′=〈V′，ε′〉進(jìn)行最大流/最小割求得最小能量是相同的，文獻(xiàn)［7］已給出相關(guān)證明.圖1為多源點(diǎn)和多匯點(diǎn)圖，根據(jù)上述規(guī)則可以簡(jiǎn)化成單源點(diǎn)和單匯點(diǎn)圖，如圖2所示.

圖1 多源點(diǎn)和多匯點(diǎn)圖

圖2 簡(jiǎn)化源點(diǎn)和匯點(diǎn)后的圖

1.3 最大流/最小割

最大流/最小割算法是用來求流量網(wǎng)絡(luò)的最小能量或近似最小能量的算法.給定一個(gè)流量網(wǎng)絡(luò)G=(V，ε)，V是頂點(diǎn)集，ε是相鄰頂點(diǎn)邊的集合，源點(diǎn)s，匯點(diǎn)t，f是邊的流量函數(shù).

最大流/最小割的基本思想:分別以源點(diǎn)s和匯點(diǎn)t來構(gòu)造搜索樹S和T，通過生長(zhǎng)樹S和T來尋找增廣路徑［8］.它可以分成4步來做:

1)尋找增廣路徑(growth stage);

2)處理增廣路徑(augmentation stage);

3)組合深林成樹(adoption stage);

4)若步驟1)中尋找不到增廣路徑，則找到最小割，算法退出;否則，繼續(xù)步驟1)～步驟3).

由上述算法思想可以找到一條從源點(diǎn)s到匯點(diǎn)t的增廣路徑，如圖3所示.與源點(diǎn)和匯點(diǎn)均不相連的點(diǎn)被歸為背景區(qū)域.

圖3 從源點(diǎn)s到匯點(diǎn)t的一條路徑

1.4 分水嶺算法

分水嶺分割思想來源于地形學(xué)，在眾多現(xiàn)有的順序分水嶺算法中，Vincent與Soille在1991年提出的基于沉浸模擬的算法是最有效的算法，也是目前應(yīng)用最廣的算法.把圖像看作地形學(xué)上被水覆蓋的自然地貌，將圖像中各點(diǎn)的梯度值視為該點(diǎn)的高度，在圖像的每個(gè)區(qū)域的最小值的位置打一個(gè)洞，然后向圖像形成的地表面中緩慢注水，水面將逐漸浸沒地面，從而形成一個(gè)個(gè)集水盆地.在來自兩個(gè)不同集水盆地的水將要發(fā)生匯合，則在匯合處建一水壩.在浸沒過程的最后，每個(gè)集水盆地最終都會(huì)被水壩包圍.集水盆地的邊界，也就是水壩，則為圖像的分水嶺.

算法分為排序和泛洪兩部分.文獻(xiàn)［5］中有相關(guān)算法的偽代碼.

2 基于圖割和分水嶺的算法

利用分水嶺和圖割對(duì)圖像進(jìn)行分割，算法步驟如下:

1)在待分割的目標(biāo)物體內(nèi)用框圖選擇一個(gè)區(qū)域，作為內(nèi)輪廓I;

2)用框圖框住待分割的物體，作為外輪廓O;

3)從內(nèi)輪廓到外輪廓，利用分水嶺對(duì)圖像進(jìn)行預(yù)分割，從而得到一系列的小區(qū)域;

4)把含有內(nèi)輪廓的小區(qū)域作為源點(diǎn)s，含有外輪廓的小區(qū)域作為匯點(diǎn)t，兩相鄰像素間的流量關(guān)系為f，轉(zhuǎn)換成單源點(diǎn)和單匯點(diǎn)，構(gòu)成新的圖G;

5)用最大流/最小割求解圖的最小能量.

圖4～圖10為算法的執(zhí)行步驟.其中圖5為了方便后面的運(yùn)算假設(shè)了圖中的區(qū)域標(biāo)號(hào).

2.1 圖像預(yù)處理分割建立圖結(jié)構(gòu)

圖4 內(nèi)輪廓和外輪廓

圖5 對(duì)圖像進(jìn)行預(yù)處理

圖6 對(duì)于分割后的區(qū)域建立圖結(jié)構(gòu)

圖7 單源點(diǎn)和單匯點(diǎn)圖 圖8 最大流/最小割

圖9 切割圖(虛線為切割線)

圖10 分割圖

利用分水嶺算法對(duì)內(nèi)外輪廓之間的區(qū)域進(jìn)行預(yù)分割后，根據(jù)這些分割后的小區(qū)域建立新的圖結(jié)構(gòu)G.

首先設(shè)置一個(gè)二維數(shù)組L(l，w)=3，記錄圖像中每個(gè)點(diǎn)初始狀態(tài)下的標(biāo)號(hào)，默認(rèn)狀態(tài)下每個(gè)點(diǎn)的標(biāo)號(hào)為3，l和w分別為圖像的長(zhǎng)度和寬度.選擇待分割物體的內(nèi)輪廓，同時(shí)設(shè)置L(i，j)=0，(i，j)為內(nèi)輪廓上的任意一個(gè)像素的位置;選擇待分割物體外輪廓，同時(shí)設(shè)置 L(i，j)=1，(i，j)為外輪廓上的任意一個(gè)像素的位置.內(nèi)外輪廓之間的區(qū)域L(i，j)=2為待分割區(qū)域，則 L(i，j)∈{0，1，2，3}.L(i，j)=3 的區(qū)域?qū)⒉挥杩紤].

應(yīng)用分水嶺算法對(duì)圖像L(i，j)=2的區(qū)域進(jìn)行預(yù)分割后，圖像會(huì)分成n個(gè)小區(qū)域，每個(gè)小區(qū)域具有不同的標(biāo)號(hào) f［k］，k=1，2，…，n，對(duì)這些小區(qū)域建立圖結(jié)構(gòu)，即尋找每個(gè)小區(qū)域鄰接區(qū)域的過程，算法步驟如下:

1)k=1;

2)對(duì)標(biāo)號(hào)為k的區(qū)域，遍歷所有標(biāo)號(hào)為k的元素，并尋找它們的鄰接點(diǎn)，如果鄰接點(diǎn)的標(biāo)號(hào)為x，那么標(biāo)號(hào)為k的區(qū)域就和標(biāo)號(hào)為x的區(qū)域鄰接;遍歷完標(biāo)號(hào)為k的區(qū)域，找到其所有的鄰接區(qū)域;如果標(biāo)號(hào)為k的區(qū)域中的像素存在L(i，j)=0的像素點(diǎn)，那么標(biāo)號(hào)為k的區(qū)域就是源點(diǎn);如果存在L(i，j)=1，那么標(biāo)號(hào)為k的區(qū)域就是匯點(diǎn);這樣可以遍歷下一個(gè)標(biāo)號(hào)的區(qū)域;

3)k=k+1，如果 k＞n，算法結(jié)束，否則轉(zhuǎn)向步驟2)繼續(xù)尋找下一個(gè)標(biāo)號(hào)的鄰接點(diǎn).

對(duì)內(nèi)外輪廓之間的區(qū)域進(jìn)行預(yù)分割后，根據(jù)每個(gè)區(qū)域的鄰接區(qū)域關(guān)系建立新的圖G，建立各個(gè)節(jié)點(diǎn)的鄰接表.

2.2 邊的流量

圖像中的邊的鄰接關(guān)系一般用4鄰域或者8鄰域表示.為了避免二義性，本文邊的鄰接關(guān)系采用8鄰域表示法.相鄰頂點(diǎn)之間的邊的流量值代表相鄰頂點(diǎn)之間的相似度.本文采用文獻(xiàn)［8］中邊的容量公式:

其中 g(i，j)=exp(－gij(i)/maxk(gij(k)))gij(k)表示從i到j(luò)第k個(gè)像素的梯度幅度值，maxk(gij(k))表示最大的梯度幅度值;在預(yù)分割后的圖像中，gij(i)表示第i個(gè)區(qū)域所有像素的梯度幅度值的均值;C(i，j)表示兩個(gè)相鄰區(qū)域i和j的容量.

2.3 三維分割

在機(jī)器人輔助視覺定位系統(tǒng)中，重建出三維的病灶或人體組織，可以幫助醫(yī)生診斷病情，確定規(guī)劃路徑.三維分割和二維分割類似，算法如下:

1)第i張CT圖像上，利用本文二維分割的算法提取出輪廓，保存這個(gè)輪廓為Bi，i=1;

2)把輪廓Bi映射到第i+1張CT的圖像上，把Bi+w作為外輪廓Oi+1，Bi－w作為內(nèi)輪廓Ii+1;

3)在內(nèi)外輪廓之間利用本文算法對(duì)第i+1張CT分割，提取輪廓Bi+1;

4)i=i+1;重復(fù)步驟2)～3)的過程，直到遍歷完所有的CT圖像.

通過步驟1)～4)可完成對(duì)圖像的三維分割，然后把所提取出來物體的輪廓重建出三維物體.

3 實(shí)驗(yàn)與分析

本文算法是由C++和Matlab實(shí)現(xiàn)的，所有的實(shí)驗(yàn)結(jié)果都是在同一臺(tái)計(jì)算機(jī)上實(shí)現(xiàn)的，配置為:2 GB內(nèi)存，Intel(R)Core(TM)2 Duo處理器，XP操作系統(tǒng).

3.1 二維分割

使用本文算法對(duì)一系列圖像進(jìn)行分割，并在分割效果和運(yùn)行時(shí)間上與文獻(xiàn)［7］相比較.圖11為應(yīng)用本文算法對(duì)圖像分割的效果圖，圖12為文獻(xiàn)［7］中算法對(duì)圖像分割的效果圖.由圖11、圖12知，本文算法對(duì)圖像的分割效果優(yōu)于文獻(xiàn)［7］.同時(shí)，文獻(xiàn)［7］對(duì)初始輪廓的選擇很敏感，如果初始輪廓嚴(yán)重偏離待分割的物體，提取出的輪廓和待分割的物體相差很大;如果初始輪廓已經(jīng)很接近待分割物體，提取出物體的輪廓比較好.而本文算法，只要選擇的內(nèi)輪廓在待分割物體上，外輪廓在待分割物體上的外面，就可以提取出待分割物體的輪廓，因此本文算法受初始輪廓的影響較小.

表1為2種算法的運(yùn)行時(shí)間比較，可見，本文的算法要遠(yuǎn)遠(yuǎn)超過文獻(xiàn)［7］算法的效率.算法的運(yùn)行時(shí)間還和初始輪廓提取有關(guān)，當(dāng)內(nèi)外輪廓之間的區(qū)域比較大時(shí)，算法耗時(shí)較長(zhǎng)，當(dāng)內(nèi)外輪廓之間的區(qū)域比較小時(shí)，耗時(shí)較短.

圖11 本文算法分割效果圖

圖12 應(yīng)用文獻(xiàn)［7］算法的圖像分割效果

表1 運(yùn)行時(shí)間 s

3.2 三維分割

利用2.3部分的三維分割算法，可實(shí)現(xiàn)對(duì)CT序列圖像的三維分割，如圖13所示.先在第1張CT上交互的選擇內(nèi)外輪廓，如圖13a所示.利用二維分割算法可以提取目標(biāo)物體的輪廓，然后根據(jù)三維分割的算法，可以自動(dòng)提取出目標(biāo)物體在其它CT中的輪廓，如圖13b所示.

由于CT在掃描時(shí)間距一般在0.5～5 mm之間，同一個(gè)目標(biāo)物體的輪廓在2張CT之間相差不大.因此，在第1張CT的輪廓映射到第2張CT上的時(shí)候，已經(jīng)最大限度的把第2張CT上的輪廓提取出來.在盡可能小的范圍內(nèi)提取目標(biāo)物體輪廓，不僅輪廓效果提取的好，而且效率較高，從而實(shí)現(xiàn)整個(gè)肺部自動(dòng)分割.

圖13 三維分割

4 結(jié)論

本文把分水嶺和圖割結(jié)合起來，首先手動(dòng)選擇了內(nèi)外輪廓，然后對(duì)內(nèi)外輪廓部分進(jìn)行預(yù)處理，最后應(yīng)用最大流算法求最小能量.該算法能夠達(dá)到比較準(zhǔn)確的分割效果，同時(shí)還提高了運(yùn)算效率.由于本文仍然是采用區(qū)域和邊界相結(jié)合的方法，所以在分割效果上能夠比較精確的得到目標(biāo)物體的輪廓.運(yùn)算效率通過以下幾個(gè)方面得以提高:內(nèi)外輪廓的選取縮小了查找范圍，把預(yù)分割后的小區(qū)域作為一個(gè)像素點(diǎn);把多源點(diǎn)和多匯點(diǎn)結(jié)構(gòu)圖轉(zhuǎn)換成單源點(diǎn)和單匯點(diǎn).然而，由于對(duì)細(xì)小的物體選擇內(nèi)外輪廓不太方便，因此本文算法對(duì)面積較大的物體提取效果較好，而對(duì)細(xì)小物體提取輪廓并不太好.

References)

［1］曹蕾，路利軍，楊蕊夢(mèng)，等.基于區(qū)域增長(zhǎng)的肺結(jié)節(jié)自適應(yīng)形態(tài)分割［J］.南方醫(yī)科大學(xué)學(xué)報(bào)，2008，28(12):2109 －2125 Cao Lei，Lu Lijun，Yang Ruimeng，et al.Self-adapted segmentation of pulmonary nodule based on region growing［J］.J South Med Univ，2008，28(12):2109 －2125(in Chinese)

［2］趙鋒，趙榮椿.分裂-合并方法在圖象分割、目標(biāo)提取中的應(yīng)用［J］.西北工業(yè)大學(xué)學(xué)報(bào)，2000，18(1):116 －120 Zhao Feng，Zhao Rongchun.A new method for image segmentation［J］.Journal of North Western Polytechnical University，2000，18(1):116 －120(in Chinese)

［3］孔丁科，汪國(guó)昭.基于區(qū)域相似性的活動(dòng)輪廓SAR圖像分割［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2010，22(9):1554－1560 Kong Dingke，Wang Guozhao.Region-similarity based active contour model for SAR image segmentation［J］.Journal of Computer-Aided Design ＆ Computer Graphics，2010，22(9):1554 －1560(in Chinese)

［4］Yuri Boykov，Gareth Funkalea.Graph cuts and efficient N-D image segmentation［J］.International Journal of Computer Vision，2006，70(2):109 －131

［5］Vincent L，Soille P.Watersheds in digital spaces:an efficient algorithm based on immersion simulation［J］.IEEE Transactions on Pattern Analysis and Machine Inteligence，1991，13:583 －598

［6］Li Yin，Sun Jian，Tang Chikeung，et al.Lazy snapping［J］.ACM Transactions on Graphics，2004，23(3):303 －308

［7］Xu Ning，Ahuja Narendra，Bansal Ravi.Object segmentation using graph cuts based active contours［J］.Computer Vision and Image Understanding，2007，107:210 －224

［8］Boykov Yuri，Kolmogorov Vladimir.An experimental comparison of min-cut max-flow algorithm of energy minimization in vision［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，2004，26(9):1124 －1137