李海燕,鄒天寧,李支堯,張榆鋒,陳建華,施心陵
(1.云南大學 信息學院,云南 昆明650091;2.昆明醫(yī)科大學 第三附屬醫(yī)院乳腺外科,云南 昆明650106;3.昆明醫(yī)科大學 第三附屬醫(yī)院超聲科,云南 昆明650106)
超聲成像因實時、無創(chuàng)、低價及對人體危害小等特點被廣泛用于臨床醫(yī)學診斷。然而,超聲成像由于其固有的干涉現(xiàn)象導致圖像含有大量的斑點噪聲而且信噪比和對比度較低,給超聲圖像分割帶來極大困難進而導致臨床診斷的不精確或者錯誤。因此有效的自動分割技術對于超聲圖像的后期處理及診斷具有重要意義。
在眾多的圖像分割技術中,活動輪廓模型因其良好的分割性能而被廣泛應用于醫(yī)學圖像分割。跟據(jù)圖像特征和目標先驗知識,現(xiàn)有的活動輪廓模型可分為基于邊界的活動輪廓模型[1]和基于區(qū)域的活動模型[2,3]?;谶吔绲哪P屠脠D像梯度信息控制輪廓曲線的演化停止在目標邊界,這類算法對圖像的噪聲和弱邊緣敏感,會導致過多的虛假邊緣和邊緣泄露,容易分割失敗?;趨^(qū)域的模型依靠全局信息引導輪廓的演化,一定程度上可以處理圖像的灰度不均勻問題,但由于需要把水平集函數(shù)初始化或重新初始化為一個符號距離函數(shù),使得計算成本相當高[4]。李[4,5]提出了局部像素信息最小化的活動輪廓模型,能有效地處理圖像灰度不均勻問題,但易受斑點噪聲的影響。文獻[6]提出用局部強度的均值和方差描述局部的像素分布,可以有效地克服圖像像素不均勻的問題和易受噪聲影響的問題,但局部強度的均值和方差依賴于模型的先驗知識。
針對上述活動輪廓模型存在的不足,提出了基于方差能量最小化的活動輪廓模型。該模型利用不同的高斯分布的均值和方差描述局部的強度。首先,利用內核函數(shù)定義待分割目標附近能量信息的高斯分布;然后,基于空間變量函數(shù)的局部均值和方差構建局部高斯分布能量并在全局圖像域整合局部能量,得到局部高斯分布能量。最后,將殘差能量最小化能量函數(shù)并入水平集的標準變量方程,在能量泛函極小化驅動水平集曲線演化到目標邊界的過程中,利用局部像素信息求取高斯分布的均值和方差識別不同區(qū)域亮度的差異。提出算法的創(chuàng)新點在于:①提出算法定義了雙積分,首先用水平集函數(shù)的內核函數(shù)定義該點附近的圖像信息局部高斯分布。然后,在水平集的演化公式中,將局部能量再次積分為二重積分。②提出算法根據(jù)變微分原理推導目標區(qū)域亮度高斯能量方程的均值和方差,而不依賴先驗知識。
設I:Ω→R 表示給定的圖像區(qū)域,C 是封閉的零水平集函數(shù)(Φ=0)的曲線。通過近似的Heaviside函數(shù)的平滑函數(shù)(Hi=H(Φ),Ho=1-H(Φ))把圖像區(qū)域Ω分為Ωin和Ωout兩部分。能量函數(shù)極小化時,目標邊界曲線C 近似等于零水平集函數(shù)Φ。圖像的強度P(I)用概率密度函數(shù)(pdf)的最大似然函數(shù)來表示。它的能量泛函數(shù)為
式 (1)中的第一個積分表示圖像數(shù)據(jù),第二個積分用以控制曲線平滑的演化。λ是權重系數(shù),為正常數(shù)。以u 和σ為參數(shù)的概率密度P(I(x))定義為
假設估計值來自于圖像域Ω的每個像素所圍繞的中心局部區(qū)域,則相應于能量式(1)的內部區(qū)域方程為
其中
式中:K(·)是定義在位置x 空間局部性的內核函數(shù),在本分割算法中,用標準偏差σK的高斯核函數(shù)表示,即
采用Gteaux導數(shù)計算水平集函數(shù)Φ:Ω →R 的梯度下降。對于任意的像素x,有
式中:ψ(x)可以是任意的測試函數(shù)。
由式 (7)可得下列方程
可以使用遞歸的高斯濾波實現(xiàn)式 (14),然后根據(jù)式(14)對能量Ω0即方程 (1)的第二項的歐拉-拉格朗日方程進行推導,為了使式 (1)能量最小,對水平集函數(shù)Φ使用標準梯度下降法,即
式 (15)的推導過程可參見文獻 [5],它的穩(wěn)態(tài)解的零水平集,即為所求的目標曲線C。其中δε為Heaviside函數(shù)的平滑狄拉克函數(shù)[7],e1和e2的方程為
(1)通過零水平集Φ0初始化水平集函數(shù)Φ;
(2)用式(4)和式(5)計算目標區(qū)域的均值ui()x 和方差()x ;
(3)用式(15)的偏微分方程,求解水平集函數(shù)Φn到;
(4)檢查收斂函數(shù)Φ,確定結果是否達到穩(wěn)定狀態(tài),否則返回(2)繼續(xù)循環(huán)計算。
其迭代結果的零水平集曲線就是所求的目標邊界曲線C。
(5)根據(jù)結果圖與原圖像的相對位置關系,得到的目標邊界,疊加到原圖像上。
其中,Φ0為初始條件。為了讓提出的方法可以更好的應用于超聲圖像且對初始位置不敏感,我們將初始條件設置為
式中:C0——事先定義好的常數(shù),而Ω0為初始邊界曲線,將其設置為
式中:w、h——圖像的寬、高。
用仿真超聲圖像和臨床淋巴癌超聲圖像為測試圖像,對提出分割算法的有效性和普適性進行驗證。實驗中主要通過主觀 “分割效果”和客觀的區(qū)域面積誤差參數(shù)和邊界誤差參數(shù)[9]評價算法的優(yōu)劣。此外,將提出算法與經典的分割算法:Tsai[3]等人的算法,Li.Wang[8]等人的算法和Liu[9]等人的算法以及LBF[5]算法進行比較。本文算法與另外4種算法都是在實驗平臺為Matlab7.1,計算機配置為:Pentium (R)dual-Core E5400 2.70GHz CPU,1GB內存,Windows XP操作系統(tǒng)下進行仿真的。
圖1是仿真超聲圖像的5種算法的分割結果。圖1(a)是在原圖上疊加了斑點噪聲的仿真超聲圖像,斑點產生模型詳見文獻[10]。圖1(b)~圖1(f)分別是LBF 算法、Li.Wang等人算法、Liu等人的算法、Tsai等人算法和本文算法的分割結果,分割輪廓用紅色線標注。由圖1可以看出本文算法得到較精確的分割結果。圖1(b)和圖1(c)兩種算法受噪聲影響無法得到目標輪廓。圖1(d)和圖1(e)的分割結果因受噪聲影響,在迭代次數(shù)達到1000次時曲線仍無法完全演化至目標邊界導致得到的輪廓有較大誤差。
圖1 仿真超聲圖像分割
臨床超聲圖像采集自Philips Envisor 2540A 超聲設備,由昆明醫(yī)科大學第三附屬醫(yī)院提供,在得到病人許可后使用。
圖2是提出方法對乳腺超聲圖像的分割結果,圖2(a)是初始輪廓,圖2(b)~圖2(f)分別是提出分割算法經過100,200,300,400和500次迭代后的分割結果,從圖中可看出,經過500次迭代后,將目標區(qū)域較完整地分割出來。
圖3是淋巴癌1的算法分割結果與醫(yī)生標注區(qū)域輪廓的比較。圖3(a)為醫(yī)生標注的淋巴癌輪廓。圖3(b)~圖3(e)分別為LBF算法、Li.Wang算法、Liu算法和Tsai算法的分割結果,分割輪廓用紅色線標注。從圖中可看出,圖3(b)和圖3(c)容易受到噪聲的影響,檢測出虛假邊緣,無法得到準確的目標輪廓。圖3(d)和圖3(e)算法得到的曲線無法演化至目標輪廓。相比較而言,本文算法雖然檢測出部分虛假邊緣,但是得到的目標輪廓與醫(yī)生標注的目標區(qū)域非常接近。
圖2 乳腺超聲圖像分割結果
圖3 淋巴癌1分割算法產生的區(qū)域與醫(yī)生標注輪廓的比較
圖4是淋巴癌2的算法分割結果與醫(yī)生標注輪廓的比較。圖4(a)為醫(yī)生標注的淋巴癌輪廓。圖4(b)~圖4(e)分別為LBF算法、Li.Wang算法、Liu算法和Tsai算法的分割結果。由圖3可以看出,圖(b)和圖(c)的分割結果受斑點噪聲的干擾,無法得到準確的目標輪廓,檢測出許多虛假邊界。圖(d)的分割結果得到的曲線明顯比醫(yī)生勾畫的小。圖(e)得到了較接近醫(yī)生標注的輪廓,但在左下角處檢測出虛假輪廓。相比較而言,提出算法圖(f)更接近醫(yī)生的標注區(qū)域輪廓,而且分割結果沒有包括虛假邊緣。
圖4 對淋巴癌2分割算法產生的區(qū)域與醫(yī)生標注區(qū)域輪廓的比較
圖5是淋巴癌3的算法分割結果與醫(yī)生標注區(qū)域輪廓的比較。圖5(a)為醫(yī)生標注的淋巴癌輪廓。圖5(b)~圖5(e)分別為LBF算法、Li.Wang算法、Liu算法和Tsai算法的分割結果。由圖4可以看出,圖(b)和圖(c)的結果檢測出了許多的虛假邊界,無法得到準確的分割結果。圖(d)的分割結果與醫(yī)生標注的區(qū)域相差較大,圖(e)的輪廓較接近醫(yī)生標注的輪廓,但仍檢測出少許虛假輪廓。相比較而言,提出算法圖(f)提取的輪廓更接近醫(yī)生的標注區(qū)域輪廓,而且沒有包含多余的虛假邊緣。
圖5 對淋巴癌3分割算法產生的區(qū)域與醫(yī)生標注區(qū)域輪廓比較
為客觀比較提出分割算法的合理性,采用區(qū)域面積誤差參數(shù)和邊界誤差參數(shù)[9]作為指標,對本文算法與另外4種對照算法的分割結果與醫(yī)生標注的病變區(qū)域輪廓進行定量比較,判定分割結果的準確性。
3.3.1 區(qū)域面積誤差參數(shù)
劉等[9]提出用區(qū)域面積誤差評估有多少病變區(qū)域被生成的病變區(qū)域正確或是錯誤覆蓋。TP 表示真陽性的面積比,F(xiàn)P 表示假陽性的面積比,F(xiàn)N 表示假陰性面積比和SI表示相似性面積比。計算方法如下
式中:Am——真實輪廓或醫(yī)生手動畫出的病變區(qū)域的像素集合。Aa——分割算法生成的病變區(qū)域像素集合。當TP百分比較大時,更多的病變區(qū)域被生成的病變區(qū)域覆蓋。當FP 百分比較小時,較少的非病變區(qū)域被生成的病變區(qū)域覆蓋。當SI 百分比較大時,生成的區(qū)域更相似于醫(yī)生手工畫的區(qū)域。
3.3.2 邊界誤差參數(shù)
本文使用 Hausdorff 距離(HD)和絕對平均距離(MD)[9]兩個邊界誤差指標分析由分割算法自動得到的輪廓與放射科醫(yī)生手畫的輪廓差異對比,HD 是兩輪廓之間最大誤差,MD 是兩輪廓的平均誤差,HD 和MD 越小說明分割的結果越好。
設定真實的或醫(yī)生勾畫的邊界和算法分割的邊界分別為Q = {q1,q2,…,qr}和P = {p1,p2,…,pu},qr和pn分別是Q 和P 輪廓上的點。則輪廓P 上的點到輪廓Q 的最短距離為
式中:‖·‖ ——二維的歐氏距離。HD 和MD 定義為
式中:γ、μ——輪廓Q 和P 邊界像素的數(shù)量。HD 和MD 的一致性規(guī)范為
提出算法與4種對比算法的區(qū)域面積誤差參數(shù)和邊界誤差參數(shù)比較見表1。
表1 本文算法與另外4種算法的客觀指標比較
由表1知,本文的算法得到的TP 和SI 百分率比其它4種算法得到的更高,而Hausdorff距離和邊界誤差參數(shù)比對比算法得到的低,說明本文算法分割的輪廓比其它算法得到的輪廓更接近于真實的目標輪廓或醫(yī)生勾畫的輪廓。
本文提出了一種新的超聲圖像分割算法。提出算法根據(jù)不同的區(qū)域的均值和方程建立其能量函數(shù)。提出算法充分利用了全局變量驅使能量最小化,實現(xià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.