方差能量最小化活動輪廓模型超聲圖像分割

李海燕，鄒天寧，李支堯，張榆鋒，陳建華，施心陵

（1.云南大學(xué) 信息學(xué)院，云南 昆明650091；2.昆明醫(yī)科大學(xué) 第三附屬醫(yī)院乳腺外科，云南 昆明650106；3.昆明醫(yī)科大學(xué) 第三附屬醫(yī)院超聲科，云南 昆明650106）

0 引 言

超聲成像因?qū)崟r(shí)、無創(chuàng)、低價(jià)及對人體危害小等特點(diǎn)被廣泛用于臨床醫(yī)學(xué)診斷。然而，超聲成像由于其固有的干涉現(xiàn)象導(dǎo)致圖像含有大量的斑點(diǎn)噪聲而且信噪比和對比度較低，給超聲圖像分割帶來極大困難進(jìn)而導(dǎo)致臨床診斷的不精確或者錯(cuò)誤。因此有效的自動分割技術(shù)對于超聲圖像的后期處理及診斷具有重要意義。

在眾多的圖像分割技術(shù)中，活動輪廓模型因其良好的分割性能而被廣泛應(yīng)用于醫(yī)學(xué)圖像分割。跟據(jù)圖像特征和目標(biāo)先驗(yàn)知識，現(xiàn)有的活動輪廓模型可分為基于邊界的活動輪廓模型［1］和基于區(qū)域的活動模型［2，3］?；谶吔绲哪Ｐ屠脠D像梯度信息控制輪廓曲線的演化停止在目標(biāo)邊界，這類算法對圖像的噪聲和弱邊緣敏感，會導(dǎo)致過多的虛假邊緣和邊緣泄露，容易分割失敗?；趨^(qū)域的模型依靠全局信息引導(dǎo)輪廓的演化，一定程度上可以處理圖像的灰度不均勻問題，但由于需要把水平集函數(shù)初始化或重新初始化為一個(gè)符號距離函數(shù)，使得計(jì)算成本相當(dāng)高［4］。李［4，5］提出了局部像素信息最小化的活動輪廓模型，能有效地處理圖像灰度不均勻問題，但易受斑點(diǎn)噪聲的影響。文獻(xiàn)［6］提出用局部強(qiáng)度的均值和方差描述局部的像素分布，可以有效地克服圖像像素不均勻的問題和易受噪聲影響的問題，但局部強(qiáng)度的均值和方差依賴于模型的先驗(yàn)知識。

針對上述活動輪廓模型存在的不足，提出了基于方差能量最小化的活動輪廓模型。該模型利用不同的高斯分布的均值和方差描述局部的強(qiáng)度。首先，利用內(nèi)核函數(shù)定義待分割目標(biāo)附近能量信息的高斯分布；然后，基于空間變量函數(shù)的局部均值和方差構(gòu)建局部高斯分布能量并在全局圖像域整合局部能量，得到局部高斯分布能量。最后，將殘差能量最小化能量函數(shù)并入水平集的標(biāo)準(zhǔn)變量方程，在能量泛函極小化驅(qū)動水平集曲線演化到目標(biāo)邊界的過程中，利用局部像素信息求取高斯分布的均值和方差識別不同區(qū)域亮度的差異。提出算法的創(chuàng)新點(diǎn)在于：①提出算法定義了雙積分，首先用水平集函數(shù)的內(nèi)核函數(shù)定義該點(diǎn)附近的圖像信息局部高斯分布。然后，在水平集的演化公式中，將局部能量再次積分為二重積分。②提出算法根據(jù)變微分原理推導(dǎo)目標(biāo)區(qū)域亮度高斯能量方程的均值和方差，而不依賴先驗(yàn)知識。

1 高斯方差能量最小化的活動輪廓模型

設(shè)I：Ω→R 表示給定的圖像區(qū)域，C 是封閉的零水平集函數(shù)（Φ＝0）的曲線。通過近似的Heaviside函數(shù)的平滑函數(shù)（Hi＝H（Φ），Ho＝1－H（Φ））把圖像區(qū)域Ω分為Ωin和Ωout兩部分。能量函數(shù)極小化時(shí)，目標(biāo)邊界曲線C 近似等于零水平集函數(shù)Φ。圖像的強(qiáng)度P（I）用概率密度函數(shù)（pdf）的最大似然函數(shù)來表示。它的能量泛函數(shù)為

式 （1）中的第一個(gè)積分表示圖像數(shù)據(jù)，第二個(gè)積分用以控制曲線平滑的演化。λ是權(quán)重系數(shù)，為正常數(shù)。以u 和σ為參數(shù)的概率密度P（I（x））定義為

假設(shè)估計(jì)值來自于圖像域Ω的每個(gè)像素所圍繞的中心局部區(qū)域，則相應(yīng)于能量式（1）的內(nèi)部區(qū)域方程為

其中

式中：K（·）是定義在位置x 空間局部性的內(nèi)核函數(shù)，在本分割算法中，用標(biāo)準(zhǔn)偏差σK的高斯核函數(shù)表示，即

采用Gteaux導(dǎo)數(shù)計(jì)算水平集函數(shù)Φ：Ω →R 的梯度下降。對于任意的像素x，有

式中：ψ（x）可以是任意的測試函數(shù)。

由式 （7）可得下列方程

可以使用遞歸的高斯濾波實(shí)現(xiàn)式 （14），然后根據(jù)式（14）對能量Ω0即方程 （1）的第二項(xiàng)的歐拉－拉格朗日方程進(jìn)行推導(dǎo)，為了使式 （1）能量最小，對水平集函數(shù)Φ使用標(biāo)準(zhǔn)梯度下降法，即

式 （15）的推導(dǎo)過程可參見文獻(xiàn) ［5］，它的穩(wěn)態(tài)解的零水平集，即為所求的目標(biāo)曲線C。其中δε為Heaviside函數(shù)的平滑狄拉克函數(shù)［7］，e1和e2的方程為

2 算法實(shí)現(xiàn)

（1）通過零水平集Φ0初始化水平集函數(shù)Φ；

（2）用式（4）和式（5）計(jì)算目標(biāo)區(qū)域的均值ui（）x 和方差（）x ；

（3）用式（15）的偏微分方程，求解水平集函數(shù)Φn到；

（4）檢查收斂函數(shù)Φ，確定結(jié)果是否達(dá)到穩(wěn)定狀態(tài)，否則返回（2）繼續(xù)循環(huán)計(jì)算。

其迭代結(jié)果的零水平集曲線就是所求的目標(biāo)邊界曲線C。

（5）根據(jù)結(jié)果圖與原圖像的相對位置關(guān)系，得到的目標(biāo)邊界，疊加到原圖像上。

其中，Φ0為初始條件。為了讓提出的方法可以更好的應(yīng)用于超聲圖像且對初始位置不敏感，我們將初始條件設(shè)置為

式中：C0——事先定義好的常數(shù)，而Ω0為初始邊界曲線，將其設(shè)置為

式中：w、h——圖像的寬、高。

3 實(shí)驗(yàn)結(jié)果與性能評估

用仿真超聲圖像和臨床淋巴癌超聲圖像為測試圖像，對提出分割算法的有效性和普適性進(jìn)行驗(yàn)證。實(shí)驗(yàn)中主要通過主觀 “分割效果”和客觀的區(qū)域面積誤差參數(shù)和邊界誤差參數(shù)［9］評價(jià)算法的優(yōu)劣。此外，將提出算法與經(jīng)典的分割算法：Tsai［3］等人的算法，Li.Wang［8］等人的算法和Liu［9］等人的算法以及LBF［5］算法進(jìn)行比較。本文算法與另外4種算法都是在實(shí)驗(yàn)平臺為Matlab7.1，計(jì)算機(jī)配置為：Pentium （R）dual－Core E5400 2.70GHz CPU，1GB內(nèi)存，Windows XP操作系統(tǒng)下進(jìn)行仿真的。

3.1 仿真超聲圖像分割

圖1是仿真超聲圖像的5種算法的分割結(jié)果。圖1（a）是在原圖上疊加了斑點(diǎn)噪聲的仿真超聲圖像，斑點(diǎn)產(chǎn)生模型詳見文獻(xiàn)［10］。圖1（b）～圖1（f）分別是LBF 算法、Li.Wang等人算法、Liu等人的算法、Tsai等人算法和本文算法的分割結(jié)果，分割輪廓用紅色線標(biāo)注。由圖1可以看出本文算法得到較精確的分割結(jié)果。圖1（b）和圖1（c）兩種算法受噪聲影響無法得到目標(biāo)輪廓。圖1（d）和圖1（e）的分割結(jié)果因受噪聲影響，在迭代次數(shù)達(dá)到1000次時(shí)曲線仍無法完全演化至目標(biāo)邊界導(dǎo)致得到的輪廓有較大誤差。

圖1 仿真超聲圖像分割

3.2 臨床超聲圖像分割

臨床超聲圖像采集自Philips Envisor 2540A 超聲設(shè)備，由昆明醫(yī)科大學(xué)第三附屬醫(yī)院提供，在得到病人許可后使用。

圖2是提出方法對乳腺超聲圖像的分割結(jié)果，圖2（a）是初始輪廓，圖2（b）～圖2（f）分別是提出分割算法經(jīng)過100，200，300，400和500次迭代后的分割結(jié)果，從圖中可看出，經(jīng)過500次迭代后，將目標(biāo)區(qū)域較完整地分割出來。

圖3是淋巴癌1的算法分割結(jié)果與醫(yī)生標(biāo)注區(qū)域輪廓的比較。圖3（a）為醫(yī)生標(biāo)注的淋巴癌輪廓。圖3（b）～圖3（e）分別為LBF算法、Li.Wang算法、Liu算法和Tsai算法的分割結(jié)果，分割輪廓用紅色線標(biāo)注。從圖中可看出，圖3（b）和圖3（c）容易受到噪聲的影響，檢測出虛假邊緣，無法得到準(zhǔn)確的目標(biāo)輪廓。圖3（d）和圖3（e）算法得到的曲線無法演化至目標(biāo)輪廓。相比較而言，本文算法雖然檢測出部分虛假邊緣，但是得到的目標(biāo)輪廓與醫(yī)生標(biāo)注的目標(biāo)區(qū)域非常接近。

圖2 乳腺超聲圖像分割結(jié)果

圖3 淋巴癌1分割算法產(chǎn)生的區(qū)域與醫(yī)生標(biāo)注輪廓的比較

圖4是淋巴癌2的算法分割結(jié)果與醫(yī)生標(biāo)注輪廓的比較。圖4（a）為醫(yī)生標(biāo)注的淋巴癌輪廓。圖4（b）～圖4（e）分別為LBF算法、Li.Wang算法、Liu算法和Tsai算法的分割結(jié)果。由圖3可以看出，圖（b）和圖（c）的分割結(jié)果受斑點(diǎn)噪聲的干擾，無法得到準(zhǔn)確的目標(biāo)輪廓，檢測出許多虛假邊界。圖（d）的分割結(jié)果得到的曲線明顯比醫(yī)生勾畫的小。圖（e）得到了較接近醫(yī)生標(biāo)注的輪廓，但在左下角處檢測出虛假輪廓。相比較而言，提出算法圖（f）更接近醫(yī)生的標(biāo)注區(qū)域輪廓，而且分割結(jié)果沒有包括虛假邊緣。

圖4 對淋巴癌2分割算法產(chǎn)生的區(qū)域與醫(yī)生標(biāo)注區(qū)域輪廓的比較

圖5是淋巴癌3的算法分割結(jié)果與醫(yī)生標(biāo)注區(qū)域輪廓的比較。圖5（a）為醫(yī)生標(biāo)注的淋巴癌輪廓。圖5（b）～圖5（e）分別為LBF算法、Li.Wang算法、Liu算法和Tsai算法的分割結(jié)果。由圖4可以看出，圖（b）和圖（c）的結(jié)果檢測出了許多的虛假邊界，無法得到準(zhǔn)確的分割結(jié)果。圖（d）的分割結(jié)果與醫(yī)生標(biāo)注的區(qū)域相差較大，圖（e）的輪廓較接近醫(yī)生標(biāo)注的輪廓，但仍檢測出少許虛假輪廓。相比較而言，提出算法圖（f）提取的輪廓更接近醫(yī)生的標(biāo)注區(qū)域輪廓，而且沒有包含多余的虛假邊緣。

圖5 對淋巴癌3分割算法產(chǎn)生的區(qū)域與醫(yī)生標(biāo)注區(qū)域輪廓比較

3.3 客觀分割結(jié)果比較

為客觀比較提出分割算法的合理性，采用區(qū)域面積誤差參數(shù)和邊界誤差參數(shù)［9］作為指標(biāo)，對本文算法與另外4種對照算法的分割結(jié)果與醫(yī)生標(biāo)注的病變區(qū)域輪廓進(jìn)行定量比較，判定分割結(jié)果的準(zhǔn)確性。

3.3.1 區(qū)域面積誤差參數(shù)

劉等［9］提出用區(qū)域面積誤差評估有多少病變區(qū)域被生成的病變區(qū)域正確或是錯(cuò)誤覆蓋。TP 表示真陽性的面積比，F(xiàn)P 表示假陽性的面積比，F(xiàn)N 表示假陰性面積比和SI表示相似性面積比。計(jì)算方法如下

式中：Am——真實(shí)輪廓或醫(yī)生手動畫出的病變區(qū)域的像素集合。Aa——分割算法生成的病變區(qū)域像素集合。當(dāng)TP百分比較大時(shí)，更多的病變區(qū)域被生成的病變區(qū)域覆蓋。當(dāng)FP 百分比較小時(shí)，較少的非病變區(qū)域被生成的病變區(qū)域覆蓋。當(dāng)SI 百分比較大時(shí)，生成的區(qū)域更相似于醫(yī)生手工畫的區(qū)域。

3.3.2 邊界誤差參數(shù)

本文使用 Hausdorff 距離（HD）和絕對平均距離（MD）［9］兩個(gè)邊界誤差指標(biāo)分析由分割算法自動得到的輪廓與放射科醫(yī)生手畫的輪廓差異對比，HD 是兩輪廓之間最大誤差，MD 是兩輪廓的平均誤差，HD 和MD 越小說明分割的結(jié)果越好。

設(shè)定真實(shí)的或醫(yī)生勾畫的邊界和算法分割的邊界分別為Q ＝ ｛q1，q2，…，qr｝和P ＝ ｛p1，p2，…，pu｝，qr和pn分別是Q 和P 輪廓上的點(diǎn)。則輪廓P 上的點(diǎn)到輪廓Q 的最短距離為

式中：‖·‖ ——二維的歐氏距離。HD 和MD 定義為

式中：γ、μ——輪廓Q 和P 邊界像素的數(shù)量。HD 和MD 的一致性規(guī)范為

提出算法與4種對比算法的區(qū)域面積誤差參數(shù)和邊界誤差參數(shù)比較見表1。

表1 本文算法與另外4種算法的客觀指標(biāo)比較

由表1知，本文的算法得到的TP 和SI 百分率比其它4種算法得到的更高，而Hausdorff距離和邊界誤差參數(shù)比對比算法得到的低，說明本文算法分割的輪廓比其它算法得到的輪廓更接近于真實(shí)的目標(biāo)輪廓或醫(yī)生勾畫的輪廓。

4 結(jié)束語

本文提出了一種新的超聲圖像分割算法。提出算法根據(jù)不同的區(qū)域的均值和方程建立其能量函數(shù)。提出算法充分利用了全局變量驅(qū)使能量最小化，實(shí)現(xiàn)了無需人工干預(yù)的自動分割。將提出算法用于含有斑點(diǎn)噪聲的仿真超聲圖像和臨床超聲圖像進(jìn)行分割測試，并與經(jīng)典的基于活動輪廓模型的分割算法進(jìn)行主客觀評價(jià)，實(shí)驗(yàn)結(jié)果表明：提出算法有較好的分割性能，驗(yàn)證了算法的有效性和普適性。

［1］Li C，Xu C，Gui C，et al.Distance regularized level set evolution and its application to image segmentation ［J］.IEEE Transaction on Image Processing，2010 （19）：3243－3254.

［2］Mao H，Liu H，Shi P.Neighbor－constrained active contours without edges ［C］／／Anchorage，Alaska，USA：CVPR，2008：1－7.

［3］Tsai A，Yezzi A，Willsky AS.Curve evolution implementation of the mumford－Shah functional for image segmentation，denoising，interpolation and magnification ［J］.IEEE Transaction on Image Processing，2001，10 （8）：1169－1186.

［4］Li C，Kao C，Gore J，et al.Implicit active contours driven by local binary fitting energy ［C］／／Washington，DC，USA：CVPR，2007：1－7.

［5］Li C，Kao C，Gore J.Minimization of region－scalable fitting energy for image segmentation ［J］.IEEE Transaction on Image Processing，2008，17 （10）：1940－1949.

［6］Brox T.From pixels to regions：Partial differential equations in image analysis［D］.Germany：Saarland University，2005.

［7］Chan TF，Vese L.Active contours without edges［J］.IEEE Transaction on Image Processing，2001，10 （2）：266－277.

［8］Wang L，Li C，Sun Q.et al.Brain MR image segmentation using local and global intensity fitting active contours／surfaces［J］.Med Image Comput Comput Assist Interv，2008，11 （Pt 1）：384－392.

［9］LIU B，Cheng HD，Huang J，et al.Automated segmentation of ultrasonic breast lesions using statistical texture classification and active contour based on probability distance ［J］.Ultrasound in Medicine Biololgy，2009，35 （8）：1309－1324.

［10］Laporte C，Clark JJ，Arbel T.A fractal multi－dimensional ultrasound scatter distribution model［C］／／Arlington，Virginia，USA：4th IEEE International Symposium on Biomedical Imaging，2007：880－883.