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

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

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

0 引言

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

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

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

1 方法

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

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

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

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

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

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

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

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

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

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

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

以一幅二維采樣圖像為例（見(jiàn)圖2所示），需要分割的邊緣為像素高亮的腦骨骼組織，標(biāo)示為白色的演化曲線不斷向邊緣靠近，它周圍的黑色區(qū)域?yàn)樵撗莼€的窄帶。因?yàn)檫吘壍南袼赝哂懈叩幕叶鹊燃?jí)，且與周圍組織對(duì)比明顯，所以當(dāng)演化曲線不斷向圖像邊緣靠近時(shí)，窄帶中所包圍的平均能量也越大。當(dāng)其達(dá)到最大化時(shí)，窄帶中所包含的平均能量達(dá)到最大(圖3為該層圖像的平均帶寬能量曲線，最大值如五角星所指處)，即達(dá)到最接近圖像邊緣的程度。

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

1.2 算法分割性能評(píng)價(jià)

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

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

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

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

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

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

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

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

基于FCM預(yù)分割得到的初始水平集表面輪廓后，使用本文提出的GVF水平集演化，并通過(guò)判斷平均帶寬能量最大化作為停止的條件，分別獲得如圖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é)果評(píng)價(jià)

重疊指數(shù)的比較結(jié)果如圖7所示，OM_3D為對(duì)三維結(jié)果計(jì)算得到的重疊指數(shù)，OM_2D為分別對(duì)各層橫截面計(jì)算的平均重疊指數(shù)。從三維對(duì)比的結(jié)果看，使用本文提出的分割方法與手動(dòng)分割的法計(jì)算所得的大鼠和小鼠的重疊指數(shù)平均值分別為0.883和0.874，標(biāo)準(zhǔn)偏差分別為0.003和0.011。與之相對(duì)應(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平臺(tái)開(kāi)發(fā)的，所用的計(jì)算機(jī)配置為Intel Core 2 Quad CPU: 2.83GHz，RAM: 2 GB，Matlab版本為R2010a。圖像大小不同，其處理速度也不同，比如處理大小約為160×200×250像素的大鼠圖像，水平集演化過(guò)程需要8 min左右的時(shí)間。根據(jù)處理時(shí)間與圖像大小的比值得到大/小鼠的處理效率，詳細(xì)情況如表1所示。

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

3 討論

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

致謝

感謝上海交通大學(xué)醫(yī)學(xué)院附屬瑞金醫(yī)院小動(dòng)物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] 李傳富, 周康源, 黃丹, 等. 基于先驗(yàn)知識(shí)的顱腦CT圖像自動(dòng)化分割[J].中國(guó)醫(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.