姚倩　侯筱婷 張娟　王麗

摘 要：針對醫(yī)學(xué)圖像CT圖像像素不均勻?qū)D像局部分割算法影響較大的問題，提出一種基于Lagrangian粒子增強補種算法的混合水平集醫(yī)學(xué)CT圖像分割算法。首先，針對局部圖像的非均勻性，通過在計算水平集公式前先計算Lagrangian標記粒子來重建內(nèi)嵌交接界面，從而提高水平集算法的質(zhì)量守恒特性；其次，針對傳統(tǒng)粒子方法在處理界面奇異性和復(fù)雜幾何相關(guān)問題上的不確定性，通過增加速度矢量和單位法向量來促進奇異點和拓撲變化點速度場的收斂；最后，通過在合成數(shù)據(jù)測試集和真實CT圖像上的仿真測試表明，所提算法在邊緣分割收斂精度及運算速度上均要優(yōu)于對比算法。

關(guān)鍵詞：補種算法；奇異性；水平集； CT圖像；分割

中圖分類號：TP391 文獻標識碼：A

1 引 言

圖像處理方向?qū)W者Osher與Sethian等[1]聯(lián)合研究并提出了水平集（Level Set， LS）分割算法，主要思路是將低維閉合二維曲線映射到高維的水平集上進行處理，從而實現(xiàn)了圖像分割算法穩(wěn)定性的增強。盡管如此，圖像界面奇異性和復(fù)雜幾何相關(guān)性對水平集算法影響仍然很大，如何處理好上述問題，對于提高水平集分割算法在CT、X-ray等醫(yī)學(xué)成像中的應(yīng)用效果，至關(guān)重要，有利于增強算法在處理對比度差、像素不均勻、邊界模糊等問題并存的醫(yī)學(xué)影像成像分割問題上的表現(xiàn)[2-3]。

針對鋒利邊緣和極端變形等情況，不同學(xué)者提出各種改進方式，如文獻[4]基于歐拉-拉格朗日粒子的方法對水平集算法進行改進，該方法在保持分割界面質(zhì)量上表現(xiàn)不錯。此外，還有局部水平集方法[5-6]、間斷Galerkin方法[7-8]、梯度增強水平集方法[9-10]等。上述方法在處理鋒利邊緣和極端變形等情況有一定效果，但是當像素不均勻或交接界面存在奇異點或拓撲變化時，上述算法效果不佳。

對此，本文在文獻[4]基礎(chǔ)上，采用Lagrangian粒子增強補種算法與水平集算法相結(jié)合，通過增加速度矢量和單位法向量來促進奇異點速度場的收斂，提出了一種基于Lagrangian粒子增強補種算法的混合水平集醫(yī)學(xué)CT圖像分割算法（LPRLS）。

2 基于粒子的水平集算法

2.1 水平集算法

2.2 混合粒子水平集方法

該方法融合了歐拉水平集算法[16]和拉格朗日粒子算法[17]，采用拉格朗日無質(zhì)量粒子來糾正圖像界面待求解區(qū)域的水平集函數(shù)。兩組粒子被隨機放置在窄帶的界面附近，并吸引界面到正確的一邊，其中陽性粒子趨向>0的一邊，而陰性粒子趨向≤0的一邊，格朗日粒子將在給定速度下產(chǎn)生對流：

式中，xp為粒子位置；uxp為粒子的速度向量。由于演化方程的耗散缺失，因此流動的特征信息被完全保存。然后，利用三階TVD Runge Kutta方法來求解時間導(dǎo)數(shù)的演化方程。

每個粒子都有一個位置和范圍半徑，用來彌補因水平集函數(shù)的局限性所引起的的位置誤差。每個粒子的半徑由基于網(wǎng)格的最大值和最小值進行界定。根據(jù)公式（6）粒子半徑的最大和最小界定值，可由下式進行計算：

在界面清晰像素均勻的較容易解決圖像表面地區(qū)，水平集函數(shù)的解是足夠準確的，并且粒子不漂移距離交接界面區(qū)域過遠。然而在像素不均勻、界面模糊的不容易解決區(qū)域，水平集函數(shù)計算會產(chǎn)生質(zhì)量損失，粒子會漂移距離交接界面過遠。界面上逃脫的正（或負）粒子，基于網(wǎng)格點、粒子半徑和粒子符號可定義局部水平集函數(shù)。水平集函數(shù)集的糾正可通過對比網(wǎng)格水平集函數(shù)值和逃脫粒子的局部水平集函數(shù)值實現(xiàn)。基于網(wǎng)格水平集函數(shù)的符號，則逃脫陰陽離子的局部水平集函數(shù)可分別定義如下：

通過將>0和≤0的區(qū)域分別與網(wǎng)格水平集函數(shù)值+和-進行比較，然后基于誤差校正方法對水平集函數(shù)進行重建。上述水平集函數(shù)值+和-分別位于陽性和陰性區(qū)域，然后通過逃脫粒子計算定義的局部水平集函數(shù)值。x>0和x≤0區(qū)域的每個角落，可通過下式計算：

水平集函數(shù)值+和-可通過設(shè)置+和-的等價變量來合成為一個單一的水平集函數(shù)，該水平集都具有在每個網(wǎng)格點的最小值。變量可定義如下：

3 粒子補種方法及改進

3.1 粒子補種方法

在存在拉伸和撕裂的流動界面中，特別是不均勻區(qū)域，這樣的區(qū)域在執(zhí)行上述混合粒子水平集算法時會缺乏足夠數(shù)量的粒子產(chǎn)生。為了保持界面的準確性，在極度變形的交接界面附近進行粒子補種以保持適當?shù)牧Ｗ臃植际呛苤匾牟襟E。對此提出一種粒子補種算法，陰陽粒子補種范圍，如圖2所示：

如果粒子在距離交接界面3個網(wǎng)格點以內(nèi)的總數(shù)量大于最大值1.5Np，或者小于最小值0.5Np，式中Np取值可定義為：

Np=16， for 2D screen64， for 3D screen（15）

粒子被隨機刪除（或添加）直到顆粒的數(shù)目在規(guī)定的最大值和最小值之間，如圖2所示。補種過程不僅加強了區(qū)域特征信息實現(xiàn)對運動區(qū)域的完全描述，而且消除了不良粒子誘導(dǎo)。

在處理相關(guān)性問題中，在交接界面附近補種適量的粒子數(shù)量是很重要的。當基于水平集函數(shù)的幾何信息進行補種時，表面的數(shù)量（面積或體積）可由下式給出：

式中，Ω為操作區(qū)域；H為Heaviside輔助函數(shù)，可定義如下：

這種基于計算數(shù)或幾何信息補種程序方法最為簡單，但這種補種標準并不明確。如果特征被沖擊波或可壓縮流動等外界干擾破壞，則原有的粒子補種機制無法刪除含有錯誤信息的污染顆粒，從而影響到水平集切面的準確性。頻繁的粒子補種操作并不是解決問題的根本方案，補種操作并不能找出導(dǎo)致界面錯誤的非正常粒子，反而會刪除含有正確特征的粒子。

3.2 改進粒子補種算法

原始粒子補種算法對于處理界面奇異問題和復(fù)雜的幾何相關(guān)性問題時效果不佳，比如單獨或聯(lián)合區(qū)域劃分問題，因此需要一種更為嚴格和可靠的粒子補種方法。endprint

速度場和幾何流會引起界面的奇異性，比如角落或邊緣。通過分析速度場和界面之間的相關(guān)性可知包含錯誤信息的特定區(qū)域。當速度場（或特征線）收斂于界面上的一點，由于拉格朗日粒子特性保持的完全性，那么逃逸粒子將會沿該奇異點堆積。這些逃逸的受污染粒子將會通過補種算法清除污染粒子前的校正階段影響交接界面劃分的準確性。由于粒子方法無法區(qū)分拓撲事件，例如，合并或分割，因此拓撲結(jié)構(gòu)的變化是影響水平集交接界面準確性的另一個因素。

圖3 拓撲變化示意圖

本文求解思路是在處理奇異性和拓撲合并或破裂問題時（如圖3所示），采用求解雙曲型方程的柔性解決方案。采取的策略是：當檢測到區(qū)域中含有奇異點或拓撲結(jié)構(gòu)變化時，粒子的校正階段停止執(zhí)行，水平集函數(shù)值保持原始柔性計算值。奇異點是由速度場收斂而出現(xiàn)，并根據(jù)相應(yīng)的交接界面和檢測的速度向量和單位法向量的角度值進行評測，如下：

式中，θ表示交接界面附近節(jié)點特征聚集程度，如果所有的θ都接近于0，則說明此處存在奇異點，在奇異點處的特征被融合并且傳統(tǒng)的粒子算法不能正確計算。在這里假定θ∈-10o，10o時，該處為奇異點。

在圖3中，拓撲合并和分裂時發(fā)生在兩界面彼此非常接近時，這種拓撲結(jié)構(gòu)的變化可通過檢查周圍環(huán)境和水平集值進行檢測。如果像素點的水平集值在0，0.1Δx范圍內(nèi)，且水平集值為負值，則該位置附近將發(fā)生合并或者分裂的拓撲結(jié)構(gòu)變化。例如，在圖3中，節(jié)點Ei，j為發(fā)生拓撲結(jié)構(gòu)變化的可疑點。

綜上所述，增強粒子補種算法的計算步驟如下：

Step1：（初始化）讀取水平集待分割圖像，并初始化算法參數(shù)；

Step2：（奇異點檢測）根據(jù)3.3節(jié)檢測當前區(qū)域是否屬于奇異點或拓撲結(jié)構(gòu)變化點，若是轉(zhuǎn)Step4，否則轉(zhuǎn)Step3；

Step3：（粒子補種）根據(jù)3.1節(jié)，執(zhí)行粒子補種算法；

Step4：（粒子狀態(tài)檢測）檢查環(huán)繞點Ai-1，j-1、Bi+1，j-1、Ci+1，j和Di-1，j的粒子符號狀態(tài)；

Step5：（粒子刪除）判斷是否滿足刪除條件，環(huán)繞點符號是否為“-”，若是則刪除節(jié)點Ei，j附近的粒子。否則轉(zhuǎn)Step6；

Step6：（水平集分割）根據(jù)2.1～2.2節(jié)執(zhí)行混合粒子水平集分割算法。判斷是否滿足終止條件，是則終止算法，并輸出結(jié)果，否則轉(zhuǎn)Step2，繼續(xù)執(zhí)行粒子的補種和刪除操作。

4 仿真實驗與分析

仿真設(shè)備參數(shù)：處理器為i7 2.4GHz，內(nèi)存為 4G ddr1333。仿真平臺選取：matlab2012a。實驗對象：選取常用醫(yī)學(xué)用CT圖像分別為血管CT、腫瘤CT、乳腺CT、腦CT以及胸腔CT五幅圖像。評價指標選取圖像分割算法的識別率、運算時間及方差。

仿真參數(shù)設(shè)置：角度θ∈-10o，10o、公式（7）參數(shù)a=0.1，b=0.5，粒子補種參數(shù)Np=16。所選取的仿真對比算法為標準LS算法以及本文所提出的LPRLS算法。仿真結(jié)果如圖4（a）～（e）所示，通過對matlab水平集工具箱函數(shù)進行研究，并且結(jié)合本文所提算法理論，經(jīng)過編寫接口程序?qū)崿F(xiàn)LPRLS算法。除上述參數(shù)設(shè)置外，其他水平集參數(shù)設(shè)置同工具箱內(nèi)水平集算法參數(shù)設(shè)置。水平集算法的分割精度、時間消耗以及方差等評測指標仿真結(jié)果如表1所示，表1中的仿真結(jié)果為上述兩種算法分別獨立運行20次所求取的評價指標均值。

圖4（a）～（e）顯示的分別為標準LS算法和LPRLS算法在血管CT、腫瘤CT、乳腺CT、腦CT以及胸腔CT五組圖像上的分割效果。（a）圖血管CT圖像分割中，LPRLS算法比LS算法在主干血管的識別率上要明顯更高，但是存在的問題是，由于靈敏度過高，導(dǎo)致存在一定的識別噪聲點，如圖中散落的小圓圈；（b）圖腫瘤CT圖像分割中，LPRLS算法比LS算法更全面的對腫瘤輪廓進行分割，并且腫瘤內(nèi)部較精細的輪廓也能識別出來；（c）圖乳腺CT圖像分割中，LPRLS算法對整個乳腺的識別全面性要明顯高于LS算法，并且乳腺內(nèi)部也可較為精確的進行分割識別；（d）圖腦CT圖像分割中，腦CT圖像的紋理結(jié)構(gòu)更加復(fù)雜，LS算法的分割效果不是很理想，從分割的精細程度上要明顯不如LPRLS算法；（e）圖胸腔CT圖像分割中，LS算法在左胸中部輪廓識別中存在雜亂現(xiàn)象，整體輪廓也較雜亂，而LPRLS算法對于胸腔輪廓識別要更圓潤，識別效果更好。圖4從視覺直觀性上給出了標準LS算法和LPRLS算法的圖像分割效果對比，直觀上LPRLS算法要明顯優(yōu)于標準LS算法。表1分別給出標準LS算法和LPRLS算法在血管CT、腫瘤CT、乳腺CT、腦CT以及胸腔CT五組圖像上的分割數(shù)值對比結(jié)果。仿真對比指標

（a）血管CT圖像

（b）腫瘤CT圖像

（c）乳腺CT圖像

（d）腦CT圖像

（e）胸腔CT圖像

選取識別率、運算時間和識別率方差三種，算法各允許20次求均值。在算法分割識別率指標上，LPRLS算法要明顯高于標準LS算法，并且這種優(yōu)勢非常明顯，在該指標上LPRLS算法要平均高出20%～30%，如在血管CT、腫瘤CT、乳腺CT三組圖像中，標準LS算法的識別率近50%～60%，而LPRLS算法在這三幅圖像上的分割識別率均達到85%以上；在腦CT和胸腔CT兩組圖像分割中，標準LS算法略高達到70%以上，但是LPRLS算法在這兩組分割圖像的分割識別率均達到90%以上。在計算時間指標對比中，標準LS算法在這五組圖像上的運算時間在20s～35s之間，而LPRLS算法在這五組圖像上的分割時間則分布在16s～20s之間，可看出在計算時間方面LPRLS算法也要優(yōu)于標準LS算法。在代表算法穩(wěn)定性方面的算法方差指標中，兩者相差不大，LPRLS算法略優(yōu)于標準LS算法。

5 結(jié)束語endprint

利用Lagrangian粒子增強補種算法對水平集算法進行改進，提高算法應(yīng)對醫(yī)學(xué)圖像CT圖像像素不均勻所帶來的分割精度降低問題。通過在計算水平集公式前先計算Lagrangian標記粒子來重建內(nèi)嵌交接界面，從而提高水平集算法的質(zhì)量守恒特性；通過增加速度矢量和單位法向量來促進奇異點和拓撲變化點速度場的收斂，來處理界面奇異性和復(fù)雜幾何相關(guān)問題。算法仿真結(jié)果顯示，所提算法比原始算法性能大為提升。從仿真結(jié)果可看出，存在的不足是，算法存在誤識別率，下一步將主要針對此問題進行進一步分析和研究。

參考文獻

[1] OSHER S，SETHIAN J A. Fronts propagating with curvaturedependent speed： algorithms based on HamiltonJacobi formulations[J]. Journal of Computational Physics，1988，79（1）：12-49.

[2] 王斌， 李潔， 高新波. 一種基于邊緣與區(qū)域信息的先驗水平集圖像分割方法[J]. 計算機學(xué)報， 2012， 35（5）： 1067-1071.

[3] 薛維琴， 周志勇， 張濤. 灰度不均的弱邊緣血管影像的水平集分割方法[J]. 軟件學(xué)報， 2012， 23（9）： 2489-2499.

[4] GIBOU F，MIN C，F(xiàn)EDKIW R. High resolution sharp computational methods for elliptic and parabolic problems in complex geometries[J]. Journal of Science Computing， 2013， 54（2）： 369-413.

[5] BOUMEHED M，ALSHAQAQI B，OUAMRI A. Moving Objects Localization by Local Regions Based Level Set： Application on Urban Traffic[J]. Journal of Mathematical Imaging and Vision， 2013， 46（2）： 258-274.

[6] JEFFREY C L，ZACHARY M. Level sets of the Takagi function： local level sets[J]. Monatshefte für Mathematik， 2012， 166（2）： 201-238.

[7] MARCHANDISE E，REMACLE J F，CHEVAUGEON N. A quadraturefree discontinuous Galerkin method for the level set equation[J]. Journal of Computational Physics， 2006， 212（1）： 338–357.

[8] MU L，WANG J P，WANG Y P. A computational study of the weak Galerkin method for secondorder elliptic equations[J]. Numerical Algorithms， 2013， 63（4）： 753-777.

[9] NAVE J C，ROSALES R R，SEIBOLD B. A gradientaugmented level set method with an optimally local coherent advection scheme[J]. Journal of Computational Physics， 2010， 229（10）： 3802–3827.

[10]LEE C，DOLBOW J，MUCHA P J. A narrowband gradientaugmented level set method for multiphase incompressible flow[J]. Journal of Computational Physics， 2014， 273（5）： 12–37.

[11]楊誼， 喻德曠. 基于改進水平集圖像分割方法的乳腺超聲病灶提取[J]. 計算機應(yīng)用與軟件， 2014， 31（11）： 217-221.

[12]Li Changyang， Wang Xiuying， Eberl S. A Likelihood and Local Constraint Level Set Model for Liver Tumor Segmentation from CT Volumes[J]. IEEE Transactions on Biomedical Engineering， 2013， 60（10）： 2967-2977.

[13]KACHROO P，RATLIFF L，SASTRY S. Analysis of the GodunovBased Hybrid Model for Ramp Metering and Robust Feedback Control Design[J]. IEEE Transactions on Intelligent Transportation Systems， 2014， 15（5）： 2132-2142.

[14]Zhuang Chijie， Zeng Rong， Zhang Bo. A WENO Scheme for Simulating Streamer Discharge With Photoionizations[J]. IEEE Transactions on Magnetics， 2014， 50（2）： #7007904.

[15]劉世興， 宋端， 賈林. 辛RungeKutta 方法在求解LagrangeMaxwell方程中的應(yīng)用研究[J]. 物理學(xué)報， 2013， 62（3）： #034501.

[16]OSHERA S，Cheng Litien， KANG M J. Geometric Optics in a PhaseSpaceBased Level Set and Eulerian Framework[J]. Journal of Computational Physics， 2002， 179（2）： 622-648.

[17]王金華， 沈永明， 石峰. 基于拉格朗日粒子追蹤的渤海冬季與夏季環(huán)流及影響因素[J]. 水利學(xué)報， 2011， 42（5）： 544-553.endprint