先驗(yàn)形狀約束的超聲前列腺圖像分割方法

劉維平，楊新

0 引言

前列腺癌是我國男性中發(fā)病率呈增長趨勢(shì)的惡性腫瘤，早期診斷對(duì)于及時(shí)治療有重要的意義。醫(yī)院確診前列腺癌的主要手段是經(jīng)直腸的超聲活組織穿刺取樣檢驗(yàn)，而結(jié)合多種成像模式的穿刺活檢能夠提高前列腺癌的檢出率。核磁共振成像對(duì)前列腺及其內(nèi)部結(jié)構(gòu)的成像效果遠(yuǎn)好于超聲圖像，便于醫(yī)生觀察以確定癌癥位置。若將核磁共振成像中疑似癌癥的位置映射到超聲圖像中作有針對(duì)性的穿刺，就需要對(duì)超聲和核磁共振成像的前列腺分別分割，然后通過配準(zhǔn)分割所得的輪廓使得超聲和核磁共振成像前列腺的切片對(duì)應(yīng)，方便與醫(yī)生觀察和比較。因此，超聲前列腺分割對(duì)前列腺癌癥的診斷有著重要的意義。除此之外，其在前列腺特異抗原計(jì)算、超聲聚焦、超聲引導(dǎo)的前列腺短程放療等應(yīng)用領(lǐng)域，同樣具有極大的潛力。

近年來，基于先驗(yàn)知識(shí)的分割方法在超聲前列腺分割中受到了廣泛的關(guān)注。該方法的特點(diǎn)在于在傳統(tǒng)的分割方法中引入了先驗(yàn)知識(shí)，對(duì)于待分割目標(biāo)在復(fù)雜背景下產(chǎn)生的分割結(jié)果的不確定性，可依靠目標(biāo)的先驗(yàn)知識(shí)得到解決，即利用先驗(yàn)知識(shí)對(duì)目標(biāo)進(jìn)行恢復(fù)和重構(gòu)的過程。

Shen[1]提出處理局部形變的形狀模型，采用旋轉(zhuǎn)不變特征作為超聲圖像屬性驅(qū)動(dòng)前列腺形狀模型的形變。為了避免一些優(yōu)化過程中出現(xiàn)的局部能量最小，Shen采用多分辨率的策略，在低分辨率中粗略分割，在高分辨率中精確分割。

Xu[2]提出利用超聲前列腺的先驗(yàn)直方圖和先驗(yàn)形狀模型實(shí)現(xiàn)分割。先驗(yàn)直方圖和形狀都是通過機(jī)器學(xué)習(xí)的方法從訓(xùn)練集中學(xué)習(xí)出來，通過一個(gè)手動(dòng)調(diào)節(jié)的參數(shù)平衡這兩項(xiàng)對(duì)曲線演化的影響。

Yan[3]采用穩(wěn)定、正確的特征提取的結(jié)果來引導(dǎo)形狀模型的形變。在特征提取的時(shí)候?qū)μ卣鬟M(jìn)行初步的篩選，那些不好的或者不穩(wěn)定的特征提取結(jié)果是不對(duì)形狀模型的形變產(chǎn)生任何影響的。

本文采用基于先驗(yàn)形狀的活動(dòng)輪廓模型對(duì)超聲前列腺圖像分割。將前列腺的形狀通過機(jī)器學(xué)習(xí)的方法求得形狀模型，并將其作為先驗(yàn)知識(shí)引入活動(dòng)輪廓模型，使得曲線滿足先驗(yàn)形狀的約束。

1 局部高斯擬合的活動(dòng)輪廓

活動(dòng)輪廓模型是一種圖像分割方法，其主要思想是在圖像域定義一個(gè)曲線或曲面，并且在由曲線或曲面自身相關(guān)的內(nèi)力，以及由圖像數(shù)據(jù)定義的外力的作用下移動(dòng)。通常情況下，外力的定義是基于邊緣特征或者區(qū)域特征的，它們分別采用了待分割圖像的局部和全局信息。為了兼顧全局和局部特征對(duì)演化曲線的影響，近些年來，一些學(xué)者提出對(duì)背景和目標(biāo)做局部統(tǒng)計(jì)的方法。最為代表性的思路是基于局部擬合方法的能量函數(shù)，用以解決灰度分布不均勻的圖像分割問題。

Wang[4]提出基于局部高斯擬合的活動(dòng)輪廓模型，對(duì)于圖像上的某點(diǎn)，定義該點(diǎn)處的能量函數(shù)為

其中，Kσ(x)為截?cái)嗟母咚购撕瘮?shù)，σ為方差，pi,x(I(y ))，i= 1,2，為像素的灰度在區(qū)域Ω1，Ω2中的概率密度分布。點(diǎn)x處能量函數(shù)被限制在以該點(diǎn)為中心的某個(gè)窗口中，截?cái)喔咚购说姆讲顩Q定了窗口ρ的大小，并且，圖像中的點(diǎn)y對(duì)能量函數(shù)的貢獻(xiàn)隨著y與x距離的增大而減小，該截?cái)喔咚购撕瘮?shù)如下所示：

其中，c是歸一化的常量，通常情況下，令ρ=2σ。

在Wang提出的方法中，假設(shè)概率密度分布 pi,x(I(y ))滿足高斯分布，即

其中， ui( x)和 σi(x)2分別為點(diǎn)x附近第i個(gè)局部區(qū)域灰度的均值和方差。

整個(gè)目標(biāo)的邊界通過最優(yōu)化圖像域Ω中所有局部高斯擬合能量函數(shù)得到，即

其中，φ為水平集函數(shù)，H(·)為 Heaviside函數(shù)，公式(5)是公式(4)的水平集表示。對(duì)該能量函數(shù)求極值并結(jié)合梯度下降方法可得到演化方程和某點(diǎn)處演化曲線內(nèi)部和外部灰度值的局部擬合方程。該方法能夠有效地克服圖像亮度不一致的問題，并且在紋理圖像分割中有著較好的效果。

2 先驗(yàn)形狀模型約束的活動(dòng)輪廓模型

影響超聲前列腺圖像分割的主要原因有如下幾個(gè)方面：1）斑點(diǎn)噪聲是超聲圖像中常見的干擾，導(dǎo)致分割困難。2）在超聲前列腺區(qū)域中，灰度或者紋理分布是不均勻的，有些部位灰度級(jí)較低，而另外一些部位灰度值較高；在超聲前列腺外部，有些地方的灰度分布和前列腺區(qū)域內(nèi)部的分布一樣，故難以采用全局的紋理或者灰度特征區(qū)別前列腺的內(nèi)外部分。3）超聲偽影也會(huì)對(duì)前列腺分割產(chǎn)生嚴(yán)重的影響。因此，超聲圖像的分割迄今仍然是一個(gè)難題。

在超聲前列腺圖像分割中，單純采用基于邊緣或者區(qū)域的傳統(tǒng)活動(dòng)輪廓模型的方法很難得到較好的分割結(jié)果，因此引入前列腺的先驗(yàn)性形狀作為傳統(tǒng)活動(dòng)輪廓模型的約束。對(duì)于先驗(yàn)形狀約束的圖像分割方法，首先要將待分割目標(biāo)形狀的訓(xùn)練集配準(zhǔn)，然后將配準(zhǔn)后的形狀采用機(jī)器學(xué)習(xí)的方法提取形狀特征并且對(duì)形狀建模，最后將形狀模型嵌入到傳統(tǒng)活動(dòng)輪廓模型當(dāng)中。

主成分分析方法(Principle Component Analysis, PCA)常常用于提取形狀模型的形狀特征，而形狀形變能力主要由形狀模型形變參數(shù)的概率密度估計(jì)表示。當(dāng)目標(biāo)的形狀變化不是很大的情況下，可以假設(shè)形變參數(shù)的概率密度函數(shù)滿足高斯分布，而當(dāng)形變很大的情況下，就需要采用非參數(shù)的方法估計(jì)形變參數(shù)[5]。但是超聲前列腺形狀的形變不是嚴(yán)重，因此針對(duì)這個(gè)具體問題，采用高斯分布的形狀模型。

2.1 形狀配準(zhǔn)

形狀配準(zhǔn)是機(jī)器視覺中的一個(gè)非常復(fù)雜的問題，其基本思想可以描述如下：給定目標(biāo)形狀、當(dāng)前要配準(zhǔn)的形狀和相似性度量，求解某變換使得變換后兩個(gè)形狀的相似性度量最小。Paragios[6]采用符號(hào)距離函數(shù)描述形狀，通過最優(yōu)化兩形狀的平方誤差和從而實(shí)現(xiàn)形狀的剛性配準(zhǔn)。本章采用這種方法將形狀訓(xùn)練集中的圖像配準(zhǔn)。

假設(shè)目標(biāo)形狀為ψT(X , Y )，當(dāng)前形狀為ψC(U , V )，目的是找到兩者間的剛性變換

其中，F(xiàn)(·)為前向映射；s，θ，Tx，Ty分別為從坐標(biāo)系(U, V)前向變換到坐標(biāo)系(X, Y)的尺度參數(shù)，旋轉(zhuǎn)參數(shù)，x軸方向的平移參數(shù)，y軸方向的平移參數(shù)，使得兩個(gè)形狀的平方誤差和最小，即最優(yōu)化如下能量函數(shù)

其中，Φ為坐標(biāo)系(X, Y)空間；B(·)為后向映射，即為 F(·)的逆映射

采用梯度下降的方法對(duì)公式(12)求最優(yōu)，得到配準(zhǔn)參數(shù)的迭代方程如下

通過迭代公式(13)～(16)直至穩(wěn)定，可將兩個(gè)以符號(hào)距離函數(shù)表示的形狀配準(zhǔn)。在訓(xùn)練形狀模型的時(shí)候，給定形狀的訓(xùn)練集，選取其中任意一個(gè)形狀的作為目標(biāo)形狀，將其余的形狀與該形狀采用上述方法配準(zhǔn)即可得到配準(zhǔn)后的訓(xùn)練集。

2.2 訓(xùn)練超聲前列腺的形狀模型

其中，U的列向量描述了形狀模型變化的正交系，Λ是相應(yīng)奇異值的對(duì)角矩陣。新形狀可以用k個(gè)主分量和k維列向量系數(shù)γ計(jì)算如下：

其中，Uk是d×k維矩陣，包含了k個(gè)最大特征值對(duì)應(yīng)的U的列向量。

一個(gè)確定k的方法是選擇滿足公式

最小的k。假設(shè)參數(shù)向量γj滿足高斯分布，并且滿足-3≤γj≤3，γj是 γ的第j個(gè)分量，λj是(1 n) MMT的第 j個(gè)特征值。通過調(diào)整γ使形狀模型產(chǎn)生形變。

2.3 基于先驗(yàn)形狀模型的前列腺分割

定義基于先驗(yàn)形狀模型的活動(dòng)輪廓的能量函數(shù)如下：

其中，ETotal為總的能量函數(shù)，ELGDF為局部高斯擬合的活動(dòng)輪廓模型的能量函數(shù)，EShape為形狀模型的能量函數(shù)。分割的目的是通過最優(yōu)化以上的總的能量函數(shù)來驅(qū)動(dòng)水平集曲線演化。λ為權(quán)值系數(shù)，用以平衡兩個(gè)能量項(xiàng)對(duì)曲線演化的影響。

在基于先驗(yàn)形狀的活動(dòng)輪廓模型中，需要最優(yōu)化形狀模型與演化曲線的平方誤差和，使得形狀模型與演化曲線最匹配。因此，整個(gè)系統(tǒng)中有兩套參數(shù)需要被優(yōu)化：配準(zhǔn)參數(shù)和形變參數(shù)。構(gòu)造的能量函數(shù)如下所示：

其中，B(·)為后向映射，包括尺度s，旋轉(zhuǎn)θ，x和y方向上的平移參數(shù)分別為Tx和Ty。配準(zhǔn)參數(shù)的迭代方程如下

其中F(·)為前向映射。通過梯度下降方法得到形變參數(shù)的迭代方程：

其中，jγ為形狀模型的第j個(gè)形變參數(shù)，j=1,…,k。

對(duì)能量函數(shù)求最優(yōu)解并且利用梯度先將方法可以得到如下公式

其中，δ(φ)為Dirac函數(shù)。公式 (23) 為曲線演化方程，公式 (24)～ (27) 為計(jì)算局部高斯擬合方法的參數(shù)計(jì)算公式。

其中，α=0.1為設(shè)定的參數(shù)。

給出先驗(yàn)形狀約束的活動(dòng)輪廓模型的具體計(jì)算步驟：

1.給定演化方程的終止條件

2.先驗(yàn)形狀約束的前列腺分割算法步驟：

3.初始化活動(dòng)輪廓和形狀模型；

4.將形狀模型與演化曲線配準(zhǔn)，公式(17)～(20)；

5.優(yōu)化形變參數(shù)，公式(22)；

6.迭代曲線演化方程，公式(23)～(27)；

7.判斷曲線是否滿足終止條件(28)，如果不滿足則轉(zhuǎn)到第4步；如果滿足則結(jié)束算法。

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

超聲數(shù)據(jù)庫有12個(gè)病例，每個(gè)病例有6個(gè)2維超聲前列腺圖像的切片，每個(gè)切片的大小為 700×400。這些超聲圖像數(shù)據(jù)來自于新華醫(yī)院超聲科，并且由醫(yī)生手動(dòng)分割完畢。采用留一法，即每個(gè)病歷留下一個(gè)作為分割結(jié)果的測(cè)試集，其余的圖像的手動(dòng)分割結(jié)果作為形狀模型的訓(xùn)練集。因此，在這個(gè)模型當(dāng)中，采用60幅前列腺的形狀作為形狀模型的訓(xùn)練集，12幅圖像作為算法的測(cè)試。

將配準(zhǔn)后的前列腺形狀作為訓(xùn)練集來訓(xùn)練前列腺的形狀模型。通過分別調(diào)整公式(13) 中的jγ，j=1,…,5，得到前列腺形狀模型的形變結(jié)果，如圖1所示：

圖1 前列腺形狀模型的形變

圖1可以看出，這種形狀模型能夠較好地模擬前列腺形變。當(dāng)然，有學(xué)者提出更為復(fù)雜的形狀模型，用以提取形變較大的形狀，例如基于核密度估計(jì)的形狀模型。但是在二維前列腺分割中，基于PCA和高斯分布的形狀模型足以應(yīng)付前列腺的形變，因此，針對(duì)特定問題采用了該種形狀模型。

在超聲前列腺的分割中，通過人為在前列腺圖像中選取上下左右四點(diǎn)擬合成一個(gè)閉合區(qū)域作為活動(dòng)輪廓模型演化曲線的初始化，如圖2(a)所示。先將λ賦予較小的值，再將該值在每次迭代曲線演化方程的時(shí)候乘以一個(gè)大于 1的參數(shù)。這樣使得形狀約束逐漸作用于曲線演化，避免了能量函數(shù)的局部最小。圖2(b)～(d)為曲線演化和形狀模型配準(zhǔn)、形變的過程，可以看到曲線最終收斂到形狀模型附近，得到了較好的分割效果，如圖2(a) ～(d)所示：

圖2 前列腺分割過程

從另外一些超聲前列腺的分割結(jié)果中可以看出，盡管超聲圖像中的噪聲、偽影、灰度分布不均勻等問題很嚴(yán)重，提出的方法仍然可以較好地提取前列腺的邊緣，如圖3所示：

圖3 其它分割結(jié)果

為了對(duì)分割算法的準(zhǔn)確程度有一個(gè)客觀的評(píng)測(cè)，醫(yī)學(xué)圖像中，普遍采用的是和臨床醫(yī)生的手工分割的標(biāo)準(zhǔn)輪廓進(jìn)行比較。采用 Dice相似性系數(shù) (Dice Similarity Coefficient,DSC) 作為衡量標(biāo)準(zhǔn)。Dice相似性系數(shù)定義如下

其中，As和Ag分別是算法分割結(jié)果和醫(yī)生手動(dòng)分割結(jié)果的面積。算法的分割結(jié)果與醫(yī)生手動(dòng)分割結(jié)果的Dice相似性系數(shù)，如表1所示：

表1 算法分割結(jié)果和醫(yī)生手動(dòng)分割結(jié)果的Dice相似性系數(shù)

4 結(jié)論

提出一種半監(jiān)督的、基于先驗(yàn)形狀約束和高斯擬合活動(dòng)輪廓的分割模型，主要應(yīng)用于超聲前列腺圖像分割。通過將分割結(jié)果與醫(yī)生手動(dòng)分割結(jié)果的對(duì)比，驗(yàn)證了所提出方法的有效性。

[1]Dinggang Shen, Yiqiang Zhan, Christos Davatzikos.Segmentation of prostate boundaries from ultrsound images using statistical shape model [J].IEEE Transactions on Medical Imaging, 2003, 22 (4): 539-551.

[2]Robert S.Xu, Oleg V.Michailovich, Igor Solovey, etc.A probability tracking approach to segmentation of ultrasound prostate images using weak shape priors [C].In Proceeding of SPIE Medical Imaging, 2010.

[3]Pingkun Yan, Sheng Xu, Baris Turbey, etc.Discrete deformable model guided by partial active shape model for TRUS image segmentation [J].IEEE Transactions on Biomedical Engineering, 2010, 57 (5): 1158-1166.

[4]Li Wang, Lei He, Arabinda Mishra, etc.Active contours driven by local Gaussian distribution fitting energy [J].Signal Processing, 2009, 89 (12): 2435-2447.

[5]Weiping Liu, Yanfeng Shang, Xin Yang, etc.A shape prior constraint for implicit active contours, Pattern Recognition Letters [J].2011, 32(15): 1937-1947

[6]Nikos Paragios, Mikael Rousson, Visvanathan Ramesh.Non-rigid registration using distance functions [J].Computer Vision and Image Understanding, 2003, 89(2-3): 142-165.