謝勤嵐，潘先攀

(中南民族大學(xué) 生物醫(yī)學(xué)工程學(xué)院，武漢 430074)

人體組織器官，如肝臟、腎臟和脾的位置、體積以及形態(tài)的變化與多種疾病息息相關(guān)，分析這些器官的位置、體積與形態(tài)可為多種疾病的臨床研究提供相關(guān)支持[1].對(duì)這些組織器官的精準(zhǔn)分割是進(jìn)行相關(guān)定量分析的前提，因此，尋找一種有效的分割算法來實(shí)現(xiàn)腹部CT圖像的分割顯得極其重要.

針對(duì)醫(yī)學(xué)圖像分割問題，研究人員提出了許多分割方法.趙于前[2]在提取肝臟初始輪廓的基礎(chǔ)上，采用5個(gè)不同方向的結(jié)構(gòu)元素檢測(cè)不同方向的邊緣，并結(jié)合不同尺度的結(jié)構(gòu)元素對(duì)梯度圖像進(jìn)行加權(quán)平均，實(shí)現(xiàn)肝臟邊緣的準(zhǔn)確檢測(cè);何蘭[3]等人通過3D卷積神經(jīng)網(wǎng)絡(luò)對(duì)數(shù)據(jù)集訓(xùn)練提取肝臟，以此為形狀先驗(yàn)結(jié)合Graph cuts算法分割出肝臟，但是該算法對(duì)訓(xùn)練數(shù)據(jù)的數(shù)量需求比較大，常常會(huì)因?yàn)槿鄙僮銐虻挠?xùn)練數(shù)據(jù)而造成難優(yōu)化的問題；潘曉佩[4]提出一種基于Snake balloon模型和Canny算子的融合算法，將Canny算子提取的邊界結(jié)果作為Snake balloon模型在超聲腎臟圖像分割中的收斂信息，使Snake balloon模型的迭代輪廓停在腎臟圖像的弱邊界上，有效分割出腎臟；Jayaram[5]等人通過定向主動(dòng)外觀模型提取目標(biāo)器官的大致形狀，將該形狀作為形狀先驗(yàn)融入圖割算法，實(shí)現(xiàn)肝臟、腎臟和脾臟的分割，但主動(dòng)外觀模型對(duì)目標(biāo)特征點(diǎn)標(biāo)定的準(zhǔn)確度要求比較高.

本文根據(jù)腹部CT圖像的成像特點(diǎn)，提出了一種自適應(yīng)形狀約束的Graph cuts分割算法，實(shí)現(xiàn)從腹部CT圖像中提取肝臟、腎臟以及脾等組織器官的目標(biāo)，其方法流程如圖1.首先使用基于多圖譜配準(zhǔn)的分割方法處理待分割圖像，得到目標(biāo)區(qū)域的初始分割結(jié)果，將它作為待分割區(qū)域的形狀先驗(yàn)；然后由形狀先驗(yàn)的輪廓計(jì)算符號(hào)距離圖，將其通過形狀比較函數(shù)加入到Graph cuts能量函數(shù)中；接著根據(jù)配準(zhǔn)分割過程中產(chǎn)生的目標(biāo)概率圖選擇能量函數(shù)中每條邊的形狀約束項(xiàng)系數(shù)，實(shí)現(xiàn)自適應(yīng)形狀約束；最后使用最大流最小割算法[6]優(yōu)化Graph cuts能量函數(shù)，得到分割結(jié)果.實(shí)驗(yàn)表明，該算法可以從對(duì)比度低，灰度不均，組織邊界模糊，背景復(fù)雜的CT圖像中較為準(zhǔn)確高效地分割出目標(biāo)組織器官.

圖1 方法流程圖Fig.1 Method flow chart

1 Graph cuts算法

Graph cuts是Boykov[7]等人提出的一種交互式圖像分割算法，該算法將圖像的分割問題轉(zhuǎn)化為最小化能量泛函的問題.Graph cuts能量函數(shù)如方程(1)：

E(A)=λR(A)+B(A),

(1)

其中，R(A)表示區(qū)域項(xiàng)能量，B(A)表示邊界項(xiàng)能量.A=(A1,A2,…,AP)表示一個(gè)二值向量，AP表示給像素點(diǎn)P分配的標(biāo)記，當(dāng)P標(biāo)記為目標(biāo)，AP值為1，當(dāng)P標(biāo)記為背景，則AP值為0.λ是一個(gè)非負(fù)常數(shù)，用于平衡邊界項(xiàng)B與區(qū)域項(xiàng)R之間的相對(duì)重要性關(guān)系.方程(1)可以進(jìn)一步寫成：

(2)

其中，N為像素P的領(lǐng)域點(diǎn)的集合.區(qū)域項(xiàng)能量Rp(Ap)表示給像素P分配“前景”或“背景” 標(biāo)記所產(chǎn)生懲罰.

1.1 區(qū)域項(xiàng)能量

方程(2)中，定義區(qū)域項(xiàng)Rp(Ap)可以使用負(fù)對(duì)數(shù)似然性模型[7]，如式(3)和(4)：

Rp(o)=-lnp(Ap|o),

(3)

Rp(b)=-lnp(Ap|b),

(4)

式中o代表目標(biāo)標(biāo)記，b代表背景標(biāo)記.若P標(biāo)記為目標(biāo)，則

(5)

其中，σ表示人工交互過程中標(biāo)記的所有目標(biāo)像素的標(biāo)準(zhǔn)差，μ表示標(biāo)記的所有目標(biāo)像素的均值，給P標(biāo)記為背景則同理.

Derive the transfer function of the integrator in the z-domain as:

1.2 邊界項(xiàng)能量

方程(2)中，邊界項(xiàng)能量定義為：

(6)

dist(p,q)表示像素點(diǎn)p和q之間的歐氏距離，σ是一個(gè)非負(fù)常數(shù)，Ip、Iq分別代表p、q點(diǎn)的像素值.該邊界項(xiàng)能量表示給像素p和q分配標(biāo)簽所產(chǎn)生的懲罰.

2 自適應(yīng)形狀約束的Graph cuts算法

2.1 形狀先驗(yàn)的獲取

本文采用基于多圖譜配準(zhǔn)的分割方法[8]獲取目標(biāo)器官的形狀先驗(yàn).基于多圖譜配準(zhǔn)的分割方法，是指將多個(gè)圖譜分別與待分割圖像進(jìn)行配準(zhǔn)，獲取多組形變參數(shù)，分別對(duì)各圖譜對(duì)應(yīng)的標(biāo)記圖像進(jìn)行形變，再采用合適的圖像融合技術(shù)，對(duì)所有形變后的標(biāo)記圖像進(jìn)行融合，得到分割的結(jié)果.

圖2 基于多圖譜配準(zhǔn)的分割流程圖Fig.2 Segmentation flow chart based on multi-atlas registration

2.2 符號(hào)距離圖

為在Graph cuts算法中加入目標(biāo)器官的形狀先驗(yàn)，本文將配準(zhǔn)分割結(jié)果的輪廓作為零水平集，計(jì)算符號(hào)距離圖，將其通過形狀比較函數(shù)加入到能量方程(2)中.

符號(hào)距離函數(shù)(Signed Distance Function，SDF)定義了空間中任意一個(gè)點(diǎn)到特定區(qū)域邊界的距離，該邊界對(duì)應(yīng)于水平集函數(shù)中的零水平集.符號(hào)距離函數(shù)同時(shí)對(duì)距離的符號(hào)進(jìn)行定義：點(diǎn)在區(qū)域邊界內(nèi)部為正，外部為負(fù)，位于邊界上時(shí)為0.假設(shè)d是空間X的一種度量，那么SDF用數(shù)學(xué)公式表達(dá)為：

(7)

式中(x)為符號(hào)距離函數(shù)，Ω為空間中的某個(gè)區(qū)域，?Ω是區(qū)域Ω的邊界，d(x, ?Ω)為點(diǎn)x到區(qū)域邊界?Ω的最短距離.

本文符號(hào)距離函數(shù)的計(jì)算使用了C.R.Maurer提出的線性時(shí)間算法[9].在算法中我們使用無符號(hào)距離函數(shù)，故在所得的符號(hào)距離圖中，對(duì)每個(gè)點(diǎn)的符號(hào)距離取絕對(duì)值.如圖3所示，(a)為多圖譜配準(zhǔn)分割結(jié)果，(b)為由(a)計(jì)算的符號(hào)距離圖像，(c)為無符號(hào)距離圖像.

圖3 符號(hào)距離圖的計(jì)算Fig.3 Calculation of signed distance map

在Graph cuts算法中加入形狀約束后，能量方程(2)改可改寫為：

(8)

其中Sp,q表示形狀約束項(xiàng)能量，ω是衡量形狀約束項(xiàng)能量的在能量函數(shù)中比重的常數(shù).

在本文方法中，我們引用Freedman和Zhang[10]的方法，將形狀先驗(yàn)以水平集的形式結(jié)合到Graph cuts算法中.該方法中，

(9)

2.3 自適應(yīng)形狀約束

方程(8)中形狀約束項(xiàng)的權(quán)重系數(shù)ω可寫為：

ωp,q=e-(αp-αq)2,

(10)

其中α為目標(biāo)概率圖，αp、αq分別表示p、q點(diǎn)屬于目標(biāo)的概率，即α中對(duì)應(yīng)p、q點(diǎn)的值.

將式(9)、(10)帶入方程(8)可得新的能量函數(shù)方程：

(11)

通過使用最大流最小割算法優(yōu)化能量方程(11)，便可得到最終分割結(jié)果.

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

為驗(yàn)證本文提出的算法的可靠性及有效性，我們將分割方法在3Dircadb數(shù)據(jù)集上進(jìn)行實(shí)驗(yàn)，依次分割出肝臟、腎臟和脾臟，并將分割結(jié)果與專家手動(dòng)分割的金標(biāo)準(zhǔn)比較.3Dircadb數(shù)據(jù)集包含了20組病人的腹部CT圖像，每個(gè)病人的影像切片大小為512×512像素，層數(shù)為91～260層.本文算法代碼的實(shí)現(xiàn)基于ITK( Insight Segmentation and Registration Toolkit)醫(yī)學(xué)圖像處理軟件包.

本文使用國(guó)際通用的Dice系數(shù)(Dice coefficient)、體積相似度(Volume similarity)和假陽(yáng)性誤差(False Positive Error)[12]三種評(píng)價(jià)指標(biāo)對(duì)分割結(jié)果進(jìn)行評(píng)估，其計(jì)算公式分別如下：

(12)

(13)

(14)

其中，S表示分割結(jié)果，T表示專家手動(dòng)分割的金標(biāo)準(zhǔn)，ST為分割結(jié)果減去金標(biāo)準(zhǔn)與分割結(jié)果重疊的部分，│·│表示體素.三種評(píng)價(jià)指標(biāo)的中，Dice系數(shù)越大，表示分割結(jié)果越接近金標(biāo)準(zhǔn)，即分割結(jié)果越準(zhǔn)確，相反，體積相似度和假陽(yáng)性誤差越小，說明分割結(jié)果越準(zhǔn)確.為方便比較不同方法分割結(jié)果的VS均值，本文對(duì)體積相似度VS取絕對(duì)值.

圖4展示了部分實(shí)驗(yàn)結(jié)果.從上至下，第一行為3Dircadb數(shù)據(jù)集中第3個(gè)病人的肝臟分割，第二行為第3個(gè)病人的腎臟分割，第三行為第11個(gè)病人的脾臟分割.從左到右，(a)列為待分割的原始圖像，紅色區(qū)域?yàn)槟繕?biāo)區(qū)域，(b)列為傳統(tǒng)圖割算法的分割結(jié)果，(c)列為多圖譜配準(zhǔn)分割的結(jié)果， (d)列為本文分割方法的分割結(jié)果.

從圖4 (b)列的分割結(jié)果來看，當(dāng)目標(biāo)輪廓不規(guī)則或者背景比較復(fù)雜時(shí)，傳統(tǒng)的Graph cuts算法容易出現(xiàn)欠分割或者過分割現(xiàn)象，導(dǎo)致這種結(jié)果的原因可以從Graph cuts算法的能量方程中分析：當(dāng)邊界項(xiàng)在能量方程中占比比較高，而兩個(gè)不連通的區(qū)域同屬于一個(gè)目標(biāo)時(shí)，分割結(jié)果可能不同時(shí)包含兩個(gè)區(qū)域，這會(huì)導(dǎo)致欠分割，如(b)列中的肝臟分割；當(dāng)區(qū)域項(xiàng)在能量方程中占比比較高，而圖像中其他非目標(biāo)區(qū)域的灰度與目標(biāo)區(qū)域相似時(shí)，分割結(jié)果會(huì)包含這一非目標(biāo)區(qū)域，這會(huì)造成過分割，如(b)列中的腎臟分割和脾臟分割.在(c)列分割結(jié)果中，基于多圖譜配準(zhǔn)的分割方法雖然能大致分割出目標(biāo)區(qū)域，但是仍然不能避免過分割現(xiàn)象.從(d)列的分割結(jié)果來看，本文方法通過加入形狀先驗(yàn)約束分割結(jié)果后，欠分割和過分割現(xiàn)象得到了很大的改善，且相比配準(zhǔn)分割結(jié)果，本文方法的分割結(jié)果在輪廓上保留了更多細(xì)節(jié).

為更直觀地展示本文方法的分割效果，在原始圖像上疊加了金標(biāo)準(zhǔn)的輪廓和本文方法分割結(jié)果的輪廓，如圖5.從上至下，第1到第3行分別為肝臟、腎臟和脾臟的分割結(jié)果，從左到右分別為相應(yīng)分割結(jié)果的橫截面、矢狀面和冠狀面.綠色線表示金標(biāo)準(zhǔn)的輪廓，紅色線表示本文方法分割結(jié)果的輪廓，紅色線越接近綠色線，分割的結(jié)果越準(zhǔn)確.結(jié)果表明，本文方法分割結(jié)果的輪廓與金標(biāo)準(zhǔn)的輪廓重疊程度很高，說明本文提出的分割方法是有效的.

圖4 實(shí)驗(yàn)結(jié)果展示1Fig.4 Experimental results 1

圖5 實(shí)驗(yàn)結(jié)果2Fig.5 Experimental results 2

為分析形狀約束在Graph cuts算法中的抗噪聲干擾能力，實(shí)驗(yàn)分別使用傳統(tǒng)Graph cuts算法和本文方法對(duì)具有低對(duì)比度和高噪聲特征的合成圖像分割，并對(duì)實(shí)驗(yàn)結(jié)果對(duì)比分析.實(shí)驗(yàn)結(jié)果如圖6所示，第一行為3Dircadb數(shù)據(jù)集中第9個(gè)病人的肝臟分割，第二行為第6個(gè)病人的腎臟分割，第三行為第11個(gè)病人的脾臟分割. (a)列為低對(duì)比度的待分割圖像， (b)列為(a)中加入高斯噪聲后的圖像，(c)列為傳統(tǒng)Graph cuts算法對(duì)(b)的分割結(jié)果， (d)列為本文分割方法對(duì)(b)的分割結(jié)果.圖中綠色線表示目標(biāo)區(qū)域輪廓，紅色線表示算法分割結(jié)果的輪廓.從實(shí)驗(yàn)結(jié)果可以看出，由于圖像對(duì)比度低，傳統(tǒng)Graph cuts算法分割結(jié)果的錯(cuò)分現(xiàn)象比較嚴(yán)重，尤其是脾臟的分割，灰度相似的背景區(qū)域會(huì)被歸入目標(biāo)區(qū)域；圖像中的噪聲導(dǎo)致傳統(tǒng)Graph cuts算法分割結(jié)果產(chǎn)生許多空洞，而且背景像素去除的不是很干凈.相比之下，本文方法能較好地解決以上問題，可以得到相對(duì)理想的分割效果.

圖6 實(shí)驗(yàn)結(jié)果3Fig.6 Experimental results 3

實(shí)驗(yàn)對(duì)比了5種不同分割方法的分割效果，每種方法都對(duì)3Dircadb數(shù)據(jù)集中的20幅CT圖像(pat1-pat20)分割，并統(tǒng)計(jì)每種分割方法的分割精度均值，統(tǒng)計(jì)結(jié)果如表1.Graph cuts為傳統(tǒng)Graph cuts算法分割，Level Set[13]為水平集方法分割，MARS (Multi Atlas Registration Segmentation)為基于多圖譜配準(zhǔn)的分割，LS-GC為使用水平集分割結(jié)果作為形狀先驗(yàn)的Graph cuts分割(注：形狀約束項(xiàng)系數(shù)由模糊處理后的待分割圖像自適應(yīng)選擇)，MARS -GC為本文分割方法，即使用多圖譜配準(zhǔn)分割結(jié)果作為形狀先驗(yàn)的Graph cuts分割.從統(tǒng)計(jì)結(jié)果來看，不論是肝臟分割、腎臟分割，還是脾臟分割，相比其他幾種分割方法，本文分割方法的分割精度更高，錯(cuò)誤率更低，分割效果更好.

表1 不同分割方法的分割精度均值統(tǒng)計(jì)Tab.1 Segmentation accuracy mean statistics of different segmentation methods

4 結(jié)語

本文根據(jù)腹部CT圖像具有對(duì)比度低，灰度不均，低信噪比，不同組織之間或軟組織與病灶之間邊界模糊等特點(diǎn)，提出了一種自適應(yīng)形狀約束的Graph cuts算法，采用基于多圖譜配準(zhǔn)的分割方法分割原圖像得到初始分割結(jié)果，將初始分割結(jié)果作為形狀先驗(yàn)加入到圖割算法中，約束分割結(jié)果的形狀.實(shí)驗(yàn)結(jié)果表明，該方法能夠很好地解決傳統(tǒng)Graph cuts算法分割時(shí)出現(xiàn)的欠分割和過分割問題，且相比水平集分割和基于多圖譜配準(zhǔn)分割，本文分割方法能夠得到更好的分割結(jié)果和更高的分割精度.