基于分水嶺與圖割的自動分割方法

劉 榮

(北京航空航天大學 機械工程及自動化學院，北京100191)

彭艷敏

(天津醫科大學 醫學影像學院，天津300203)

唐 粲 程 勝

(昆山市工業技術研究院，昆山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)

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

本文對內外輪廓之間的圖像進行預分割后，對分割后的小區域利用圖割進行分割，由對每個像素操作變成對小區域操作;通過實驗驗證本文算法在運行時間和運行效果上都優于文獻［7］.

1 算法基礎

1.1 圖割的相關理論

圖像分割可以轉換成求解能量函數最小值問題.圖割是利用圖論中的最小割算法求能量函數全局最優解的算法.利用圖割來解決問題需要做到以下3點:

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

2)構造能量函數.圖割綜合了邊界和區域的信息，圖像分割問題可以轉換為求取使能量函數取到極小值時的標號場.該能量函數一般形式可以表示為

其中，Vp，q(fp，fq)表示平滑項能量，是對不連續性的一種懲罰，為相鄰像素之間能量;Dp(fp)是數據項能量，是對當前標號與觀測數據間不符合程度的一種懲罰.

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

1.2 多源點多匯點簡化圖

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

圖1 多源點和多匯點圖

圖2 簡化源點和匯點后的圖

1.3 最大流/最小割

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

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

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

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

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

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

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

圖3 從源點s到匯點t的一條路徑

1.4 分水嶺算法

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

算法分為排序和泛洪兩部分.文獻［5］中有相關算法的偽代碼.

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

利用分水嶺和圖割對圖像進行分割，算法步驟如下:

1)在待分割的目標物體內用框圖選擇一個區域，作為內輪廓I;

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

3)從內輪廓到外輪廓，利用分水嶺對圖像進行預分割，從而得到一系列的小區域;

4)把含有內輪廓的小區域作為源點s，含有外輪廓的小區域作為匯點t，兩相鄰像素間的流量關系為f，轉換成單源點和單匯點，構成新的圖G;

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

圖4～圖10為算法的執行步驟.其中圖5為了方便后面的運算假設了圖中的區域標號.

2.1 圖像預處理分割建立圖結構

圖4 內輪廓和外輪廓

圖5 對圖像進行預處理

圖6 對于分割后的區域建立圖結構

圖7 單源點和單匯點圖 圖8 最大流/最小割

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

圖10 分割圖

利用分水嶺算法對內外輪廓之間的區域進行預分割后，根據這些分割后的小區域建立新的圖結構G.

首先設置一個二維數組L(l，w)=3，記錄圖像中每個點初始狀態下的標號，默認狀態下每個點的標號為3，l和w分別為圖像的長度和寬度.選擇待分割物體的內輪廓，同時設置L(i，j)=0，(i，j)為內輪廓上的任意一個像素的位置;選擇待分割物體外輪廓，同時設置 L(i，j)=1，(i，j)為外輪廓上的任意一個像素的位置.內外輪廓之間的區域L(i，j)=2為待分割區域，則 L(i，j)∈{0，1，2，3}.L(i，j)=3 的區域將不予考慮.

應用分水嶺算法對圖像L(i，j)=2的區域進行預分割后，圖像會分成n個小區域，每個小區域具有不同的標號 f［k］，k=1，2，…，n，對這些小區域建立圖結構，即尋找每個小區域鄰接區域的過程，算法步驟如下:

1)k=1;

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

3)k=k+1，如果 k＞n，算法結束，否則轉向步驟2)繼續尋找下一個標號的鄰接點.

對內外輪廓之間的區域進行預分割后，根據每個區域的鄰接區域關系建立新的圖G，建立各個節點的鄰接表.

2.2 邊的流量

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

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

2.3 三維分割

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

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

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

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

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

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

3 實驗與分析

本文算法是由C++和Matlab實現的，所有的實驗結果都是在同一臺計算機上實現的，配置為:2 GB內存，Intel(R)Core(TM)2 Duo處理器，XP操作系統.

3.1 二維分割

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

表1為2種算法的運行時間比較，可見，本文的算法要遠遠超過文獻［7］算法的效率.算法的運行時間還和初始輪廓提取有關，當內外輪廓之間的區域比較大時，算法耗時較長，當內外輪廓之間的區域比較小時，耗時較短.

圖11 本文算法分割效果圖

圖12 應用文獻［7］算法的圖像分割效果

表1 運行時間 s

3.2 三維分割

利用2.3部分的三維分割算法，可實現對CT序列圖像的三維分割，如圖13所示.先在第1張CT上交互的選擇內外輪廓，如圖13a所示.利用二維分割算法可以提取目標物體的輪廓，然后根據三維分割的算法，可以自動提取出目標物體在其它CT中的輪廓，如圖13b所示.

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

圖13 三維分割

4 結論

本文把分水嶺和圖割結合起來，首先手動選擇了內外輪廓，然后對內外輪廓部分進行預處理，最后應用最大流算法求最小能量.該算法能夠達到比較準確的分割效果，同時還提高了運算效率.由于本文仍然是采用區域和邊界相結合的方法，所以在分割效果上能夠比較精確的得到目標物體的輪廓.運算效率通過以下幾個方面得以提高:內外輪廓的選取縮小了查找范圍，把預分割后的小區域作為一個像素點;把多源點和多匯點結構圖轉換成單源點和單匯點.然而，由于對細小的物體選擇內外輪廓不太方便，因此本文算法對面積較大的物體提取效果較好，而對細小物體提取輪廓并不太好.

References)

［1］曹蕾，路利軍，楊蕊夢，等.基于區域增長的肺結節自適應形態分割［J］.南方醫科大學學報，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］趙鋒，趙榮椿.分裂-合并方法在圖象分割、目標提取中的應用［J］.西北工業大學學報，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］孔丁科，汪國昭.基于區域相似性的活動輪廓SAR圖像分割［J］.計算機輔助設計與圖形學學報，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