【作 者】陳詩燁，王成＊，鄭秀娟2,

1 上海交通大學(xué)醫(yī)學(xué)院生物醫(yī)學(xué)工程系，上海，200025

2 上海交通大學(xué)醫(yī)學(xué)院附屬仁濟(jì)醫(yī)院核醫(yī)學(xué)科，上海，200127

0 引言

要開展鼠腦Micro-PET/CT圖像自動定量分析，首要步驟是從頭部圖像中實(shí)現(xiàn)對全腦組織的分割和提取。這一過程通常是利用CT圖像提取大腦結(jié)構(gòu)，在此基礎(chǔ)上分析PET功能圖像。盡管在人類的CT或MRI圖像中分割并提取大腦組織的方法上已經(jīng)相當(dāng)成熟，然而與人類大腦的提取相較，在Micro-PET/CT圖像中的分割提取鼠腦的難度更大。其主要原因在于：①需要提取的鼠腦組織非常小，且結(jié)構(gòu)復(fù)雜，尤其是前腦部分，即使用人眼也很難分辨；② 雖然Micro-CT擁有極高的分辨率，但圖像邊緣依然比較模糊，且不連續(xù)，存在很多斷點(diǎn)，給分割帶來很大難度。因此，為了準(zhǔn)確分割鼠腦圖像，就不能簡單套用既有的相對成熟的人腦分割方法。本文提出一種基于水平集的方法，能有效地實(shí)現(xiàn)對鼠腦組織的三維自動化分割。

水平集方法最早是由Osher和Sethian[1]在1988年提出，現(xiàn)已在圖像處理和計算機(jī)視覺等領(lǐng)域得到了廣泛的應(yīng)用[2-3]。2008年，Zhang等人[4]提出了一種基于水平集的混合分割方法，其特點(diǎn)是同時結(jié)合了邊緣和區(qū)域信息，以達(dá)到準(zhǔn)確分割的效果。其中，邊緣信息可以幫助檢測目標(biāo)對象的準(zhǔn)確位置，而利用區(qū)域信息則可以防止邊緣泄漏問題。2009年，Uberti MG[5]提出了一種基于水平集方法并且針對小鼠MRI圖像提取大腦組織的半自動方法，其分割效果也較為顯著。但是，使用這些近年來所發(fā)展的水平集方法來分割鼠腦Micro-CT圖像仍然有較多不足：① 水平集初始表面輪廓仍需要手動設(shè)置；② 由于圖像邊緣模糊，存在邊緣泄露問題；③ 由于鼠腦組織結(jié)構(gòu)復(fù)雜，水平集在前腦等窄細(xì)結(jié)構(gòu)處的演化能力不夠；④ 停止條件（每次迭代的持續(xù)時間和迭代次數(shù)）仍然需要手動設(shè)置，存在諸多不便。

本文在水平集方法的基礎(chǔ)上做了三個方面的改進(jìn)，提出了一種針對大/小鼠Micro-CT圖像的腦組織分割方法，實(shí)現(xiàn)了圖像的自動化分割處理。

1 方法

1.1 改進(jìn)水平集方法的框架

本文基于改進(jìn)水平集方法，提出了一種針對大/小鼠Micro-CT圖像的腦組織分割框架，算法流程如圖1所示。首先使用模糊C均值聚類方法(Fuzzy-C-Means,FCM)，實(shí)現(xiàn)初始水平集表面輪廓的自動化設(shè)置；接著使用梯度矢量流(Gradient Vector Flow, GVF)增加外向演化能力，使得演化表面輪廓能更準(zhǔn)確地進(jìn)入窄細(xì)處等復(fù)雜結(jié)構(gòu)；最后提出一種判斷演化表面輪廓附近平均帶寬能量是否達(dá)到最大化(Average Bandwidth Energy Maximization, ABEM)的方法，作為停止條件，防止泄漏且能同時實(shí)現(xiàn)自動化停止。

圖1 算法總流程Fig.1 Algorithms of the proposed method

1.1.1 模糊C均值聚類方法的預(yù)分割

FCM方法是由Bezdek等[6]提出，該算法具有良好的局部收斂性和分割效果，被成功地用于圖像分割中[7-8]。因為腦組織和周圍組織存在差別，所以使用FCM方法，并結(jié)合閾值和形態(tài)學(xué)方法，可以獲取大致的腦部形態(tài)結(jié)構(gòu)，作為自動預(yù)分割。先通過設(shè)置一個參考閾值（該閾值由已經(jīng)發(fā)展成熟的人腦不同組織的CT值等先驗知識可以確定[9]），分割出部分頭骨組織；再通過膨脹、腐蝕和填充孔洞等形態(tài)學(xué)方法，提取包含整個腦組織的圖像區(qū)域，去除部分周圍組織的干擾；接著利用FCM方法對閾值和形態(tài)學(xué)方法分割得到的一次結(jié)果進(jìn)行再分割，獲得更為精確的腦組織區(qū)域。但是，經(jīng)過FCM分割的圖像邊緣不夠連續(xù)，周邊有部分噪聲，圖像上有很多孔洞。為了有效解決上述問題，再次利用形態(tài)學(xué)方法對圖像進(jìn)行處理，其結(jié)果作為下一步水平集演化的初始圖像。

1.1.2 梯度矢量流的水平集演化

經(jīng)典水平集演化是以邊界上點(diǎn)的演化速度作為動力的。在圖像處理上，常以灰度梯度作為計算速度的要素，據(jù)此水平集方法的演化方程可以簡化為：

其中，I是需要分割的圖像，μ表示灰度級下限的預(yù)定義參數(shù)，φ代表活動輪廓曲線C={x│φ(x)=0}的初始水平集，H(φ)是Heaviside 方程，其定義為是圖像梯度圖，Ω表示圖像域，α和β是用來調(diào)節(jié)方程兩個因子比重所預(yù)定義的參數(shù)。

1.1.3 平均帶寬能量最大化 (ABEM)

為了能自動獲得最佳演化結(jié)果，本文提出了一種判斷演化表面輪廓周圍一定帶寬區(qū)域的平均能量是否達(dá)到最大化，作為停止水平集繼續(xù)演化的條件。平均帶寬能量(Average Bandwidth Energy，ABE) 被定義為

其中，I(i)代表第i個像素的灰度級，n代表帶寬中的像素數(shù)。帶的寬度可根據(jù)圖像來設(shè)置，以期獲得最佳效果。通常是帶寬越寬，帶寬中所包含的圖像信息也越大，隨之計算量也會更大，但是會更加準(zhǔn)確。

以一幅二維采樣圖像為例（見圖2所示），需要分割的邊緣為像素高亮的腦骨骼組織，標(biāo)示為白色的演化曲線不斷向邊緣靠近，它周圍的黑色區(qū)域為該演化曲線的窄帶。因為邊緣的像素往往具有更高的灰度等級，且與周圍組織對比明顯，所以當(dāng)演化曲線不斷向圖像邊緣靠近時，窄帶中所包圍的平均能量也越大。當(dāng)其達(dá)到最大化時，窄帶中所包含的平均能量達(dá)到最大(圖3為該層圖像的平均帶寬能量曲線，最大值如五角星所指處)，即達(dá)到最接近圖像邊緣的程度。

圖2 平均帶寬能量最大化的原理Fig.2 Principle of ABEM

1.2 算法分割性能評價

本文采用基于區(qū)域法的評價方法，即計算兩個對應(yīng)區(qū)域的重疊比例，這也是國際上公認(rèn)的評價方法之一（Uberti MG[5]）。重疊指數(shù)定義為

圖3 平均帶寬能量曲線Fig.3 Energy curve of ABEM

其中 V1和 V2為兩幅需要對照的分割結(jié)果圖像，Vo是兩幅圖像的重疊區(qū)域。 OM 的范圍在0到1之間，0代表兩幅圖像之間沒有重疊部分，而1則代表兩幅圖像完全重疊。

為了評價本文所提出方法的分割效果，將手動分割的結(jié)果作為評價本文所提出方法是否分割準(zhǔn)確的標(biāo)準(zhǔn)，公式(5)中的V2由手動分割獲得，采用了Analyze?軟件[12]中勾畫感興趣區(qū)域的功能進(jìn)行一層層的分割并整合而得。V1則分別采用本文所提出的方法和Zhang等人在其論文中所提到的混合分割方法(Hybrid Level-Set，HLS)而得。2D/3D image segmentation toolbox?是Zhang等人在其論文中所提到的基于Matlab所開發(fā)的工具箱，用該方法獲得的圖像作為本文所提出的圖像處理方法的對照。為了使用的方便性，其中初始水平集表面輪廓的獲取也采用了利用FCM方法，其他步驟和方法基本參照原方法，沒有變動。方法中所設(shè)定的參數(shù)也通過反復(fù)試驗調(diào)節(jié)所得。

2 實(shí)驗結(jié)果

2.1 初始水平集表面輪廓的獲取

從圖4可以看出，通過使用基于FCM方法自動獲取的初始水平集表面輪廓（圖4中以深色曲線表示），已經(jīng)非常接近腦組織邊界。但是，前腦部分以及腦部和頸部連接處的分割效果還是不夠。以大鼠一為例，經(jīng)過FCM預(yù)分割的初始水平集表面輪廓與手工分割的結(jié)果作對照，其重疊指數(shù)為0.628。此與重疊指數(shù)為0.881的最終分割結(jié)果相比，還有一定距離，所以需要利用水平集方法作進(jìn)一步的演化處理。

2.2 水平集演化的最終結(jié)果

基于FCM預(yù)分割得到的初始水平集表面輪廓后，使用本文提出的GVF水平集演化，并通過判斷平均帶寬能量最大化作為停止的條件，分別獲得如圖5和圖6的三維和二維圖像。

圖4 初始水平集分割結(jié)果的二維采樣例圖Fig.4 Sampled axial slices of the obtained initial level-set surface

圖5 最終分割結(jié)果的三維例圖Fig.5 Segmentation result in 3D

圖6 最終分割結(jié)果的二維例圖Fig.6 Sampled axial slices of segmentation result

2.3 分割結(jié)果評價

重疊指數(shù)的比較結(jié)果如圖7所示，OM_3D為對三維結(jié)果計算得到的重疊指數(shù)，OM_2D為分別對各層橫截面計算的平均重疊指數(shù)。從三維對比的結(jié)果看，使用本文提出的分割方法與手動分割的法計算所得的大鼠和小鼠的重疊指數(shù)平均值分別為0.883和0.874，標(biāo)準(zhǔn)偏差分別為0.003和0.011。與之相對應(yīng)的，使用Zhang等所提出的方法分割的結(jié)果，重疊指數(shù)平均值分別為0.663和0.819，標(biāo)準(zhǔn)偏差分別為0.008和0.009。比較而言，大鼠的分割改善效果小鼠較好。相比Zhang的方法，大鼠的分割準(zhǔn)確率提高了33%，小鼠提高了6.7%?？梢钥吹剑疚乃岢龅姆椒ㄔ跍?zhǔn)確性上有顯著改善。

圖7 兩種分割方法重疊指數(shù)比較Fig.7 Comparison of the proposed method versus HLS by using OM

2.4 算法的處理速度

本文提出的算法是基于Matlab平臺開發(fā)的，所用的計算機(jī)配置為Intel Core 2 Quad CPU: 2.83GHz，RAM: 2 GB，Matlab版本為R2010a。圖像大小不同，其處理速度也不同，比如處理大小約為160×200×250像素的大鼠圖像，水平集演化過程需要8 min左右的時間。根據(jù)處理時間與圖像大小的比值得到大/小鼠的處理效率，詳細(xì)情況如表1所示。

表1 本文算法的處理效率Tab.1 The processing ef fi ciency of the algorithm

3 討論

隨著成像技術(shù)的發(fā)展，利用Micro-PET/CT技術(shù)來對動物模型做研究正在被越來越多研究者所采用。然而，就我們所知，還沒有專門針對大/小鼠Micro-CT圖像的腦組織自動化提取方法。因此，本文提出了基于改進(jìn)水平集方法的算法框架，可以實(shí)現(xiàn)對鼠腦組織的Micro-CT圖像自動化分割，而且在準(zhǔn)確性上有較大改善，可以作為后續(xù)進(jìn)一步量化分析的基礎(chǔ)。

3.1 關(guān)于框架流程的考慮

本文所提出的框架是在水平集演化之前加入了基于FCM方法的預(yù)分割，雖然這個部分需要耗費(fèi)額外的處理時間，但是對于整個結(jié)果的準(zhǔn)確性和便捷性都是很有用的。使用FCM方法并結(jié)合閾值和形態(tài)學(xué)的方法來獲取初始水平集表面輪廓，不僅避免了手動設(shè)置，而且獲取的初始水平集表面輪廓已經(jīng)非常接近目標(biāo)邊緣，因此為后續(xù)演化過程節(jié)約了大量時間，并且實(shí)現(xiàn)了全自動化獲取，不需要手動勾畫。

3.2 使用梯度矢量流的改善效果

引入梯度矢量流，增強(qiáng)了水平集演化的外部張力，使在圖像的窄細(xì)等復(fù)雜結(jié)構(gòu)處的演化能力大大增強(qiáng)。通過測試比較發(fā)現(xiàn)，在同樣使用FCM預(yù)分割得到的結(jié)果作為初始水平集表面，并且其他參數(shù)均設(shè)置相同的條件下，以水平集表面輪廓成功演化到前腦部分以及腦部和頸部等窄細(xì)結(jié)構(gòu)處為參照標(biāo)準(zhǔn)，使用GVF作為梯度圖通常比使用一般梯度方程要快大約10個迭代過程。從分割的整體效果看，使用一般梯度方程在演化表面輪廓成功到達(dá)前腦等區(qū)域的同時，其橫截面的其他大部分層面均出現(xiàn)了較為嚴(yán)重的泄漏，但是觀察使用GVF的結(jié)果圖，泄漏情況大為改善，且各個層面的分割進(jìn)度較為統(tǒng)一。所以實(shí)驗結(jié)果證明了GVF擴(kuò)張?zhí)荻葓鰪埩Φ淖饔?，可以加快水平集表面輪廓更快地演化到遠(yuǎn)端，同時加強(qiáng)了在窄細(xì)結(jié)構(gòu)處的演化能力。

3.3 平均帶寬能量最大化 的改善效果

通過使用平均帶寬能量最大化作為停止條件，省去了人工設(shè)置的麻煩，較好地解決了部分層面泄漏的問題。這個過程是在二維層面完成的，因為在水平集演化過程中，不同層面的演化進(jìn)度是不一樣的，雖然通過使用GVF擴(kuò)張了梯度場張力，使得水平集表面輪廓可以更快地演化到遠(yuǎn)端，例如前腦等窄細(xì)結(jié)構(gòu)處，但是還是有部分中間層面出現(xiàn)了泄漏的情況。為了使水平集表面輪廓在成功演化到遠(yuǎn)端的同時，控制中間層面的演化進(jìn)度，就需要計算各個橫截面的平均帶寬能量，并分別判斷達(dá)到最大化的時間和結(jié)果，最終得到優(yōu)化的結(jié)果。雖然在二維層面處理增加了處理的時間，但是與水平集演化所需要的時間相比，這部分時間可以忽略。

3.4 參數(shù)的優(yōu)化選擇

在水平集演化方程中，主要涉及參數(shù)有μ、α和β。如上文所介紹的，μ為全局閾值，用來表示目標(biāo)對象邊界的低灰度級。演化方程(1)的第一部利用圖像區(qū)域信息，推動曲面輪廓向大于 的方向演化。μ的選擇不是一件容易的事，如果選擇太小，則可能會導(dǎo)致演化曲面輪廓越過邊界從而過度分割；如果選擇太大，則可能會導(dǎo)致演化曲面輪廓始終到達(dá)不了目標(biāo)邊緣。此外，由于圖像不均勻，很難確定一個特定的值賦值給。本文采用LBF(local binary fi tting energy)[13-14]，—— 一種局部二值擬合能量泛函，根據(jù)局部圖像信息設(shè)定該閾值。與全局閾值法相比，使用LBF能夠提取更準(zhǔn)確的圖像信息。

α和β是用來調(diào)節(jié)方程兩部分比重所預(yù)定義的參數(shù)。α為演進(jìn)權(quán)重指數(shù)，α越大，則水平集演化所利用的區(qū)域信息比重越大，且受μ值的影響也大。演化方程（1）的第二部分利用梯度信息，使演化曲面輪廓向梯度值大的圖像邊緣靠近，β為該部分的權(quán)重指數(shù)。關(guān)于α和β的選取，以大鼠為例，通過窮盡法選取正確的數(shù)量級區(qū)間，觀測3例圖像分割結(jié)果的平均重疊指數(shù)效果變化，比較好的結(jié)果出現(xiàn)在β=0.001，α在0.001到0.1之間，針對這個結(jié)果，再進(jìn)行一定范圍的測試，觀察重疊的效果變化，得到最優(yōu)近似值為α=0.009，β=0.001 。利用同樣的方法對小鼠的和 也做了選取，得到最優(yōu)近似值為α=0.009，β=0.001。

4 結(jié)論

本文對大/小鼠腦組織提取框架是基于水平集方法的改進(jìn)，并被驗證具有較好的準(zhǔn)確性和便捷性。通過與已有經(jīng)典方法的對比，本文所提出的方法被證明具有更好的效果。利用該方法，研究人員可以很容易地對大/小鼠Micro-CT圖像進(jìn)行自動化處理，為進(jìn)一步與圖譜（例如Digimouse[15]）配準(zhǔn)，獲取豐富的腦區(qū)信息或不同個體腦組織之間的比對，打下了基礎(chǔ)。

致謝

感謝上海交通大學(xué)醫(yī)學(xué)院附屬瑞金醫(yī)院小動物PET中心和上海交通大學(xué)MED- X研究院在本研究中所提供的幫助。

[1] Osher S, Sethian JA. Fronts propagating with curvature dependent speed: algorithms based on the Hamilton-Jacobi formulation [J].Journal of Computational Physics, 1988, 79(1): 12-49.

[2] Caselles V, Catte F, Coll T, Dibos F. A geometric model for active contours in image processing [J]. Numeric Math, 1993, 66(1): 1-31.

[3] Malladi R, Sethian JA, Vemuri BC. Shape modeling with front propagation: a level set approach [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1995, 17(2): 158-175.

[4] Zhang Y, Matuszewski BJ, Shark L, et al. Medical image segmentation using new hybrid level-set method [A]. In:Proceedings of the 2008 Fifth International Conference BioMedical Visualization: Information Visualization in Medical and Biomedical Informatics[C], 2008, 71-76.

[5] Uberti MG, Boska MD, Liu Y. A semi-automatic image segmentation method for extraction of brain volume from in vivo mouse head magnetic resonance imaging using Constraint Level Sets [J]. J Neurosci Methods, 2009, 179(2): 338-344.

[6] Bezdek JC. Pattern recognition with fuzzy objective function algorithms [M]. New York: Plenum Press, 1981, 95-107.

[7] Pham D, Prince J. An adaptative fuzzy C-means algorithm for image segmentation in the presence of intensity inhomogeneities [J].Pattern Recognition Letter, 1999, 20(1): 57-68.

[8] W.J. Chen, M.L. Giger, U. Bick. A fuzzy c-means (FCM)-based approach for computerized segmentation of breast lesions in dynamic contrast enhanced MRI images [J]. Acad. Radiol, 2006,13(1): 63-72.

[9] 李傳富, 周康源, 黃丹, 等. 基于先驗知識的顱腦CT圖像自動化分割[J].中國醫(yī)療器械雜志, 2004, 28(3):168-171.

[10] Chenyang Xu and Jerry L. Prince. Gradient Vector Flow: A New External Force for Snakes [J]. IEEE Transactions on Image Processing, 1998, 7(3): 359-369.

[11] Chenyang Xu, Prince, JL. Snakes, shapes, and gradient vector fl ow[J]. IEEE Transactions on Image Processing, 1998, 7(3):359-369.

[12] Robb RA. The biomedical imaging resource at mayo clinic [J].IEEE Transactions on Medical Imaging, 2001, 20(9):854-867.

[13] C M Li, C Kao, J Gore, et al. Implicit active contours driven by local binary fi tting energy [A]. In: IEEE Conference on Computer Vision and Pattern Recognition [C]. Minneapolis, 2007, 1-7.

[14] Qingqi Hong, Qingde Li, Jie Tian. Local hybrid level-set method for mra image segmentation[A]. In: 2010 10th IEEE International Conference on Computer and Information Technology (CIT 2010).Bradford [C], 2010, 1397- 1402.

[15] Dogdas B, Stout D, Chatziioannou AF, et al. Digimouse: a 3D whole body mouse atlas from CT and cryosection data [J]. Physics in Medicine and Biology, 2007, 52(3), 577-587.