融合分區(qū)與Canny泛函的水平集對猴腦提取的研究*

郭晉秀，張月芳，鄧紅霞，李海芳

(太原理工大學(xué)信息與計算機學(xué)院，山西 晉中 030600)

1 引言

腦部提取是將大腦的非腦組織剝離，提取出腦組織的過程。數(shù)據(jù)的采集使用磁共振成像技術(shù)，但磁場強度過大會對人腦造成損傷，因此研究人員選用了能承受較大磁場強度的獼猴進行研究。獼猴是研究人類的天然模型，它具有類似于人腦的組織結(jié)構(gòu)及其它諸多保守性組織特征[1]。但是，獼猴大腦比人腦體積小很多，眼睛周圍包裹著的脂肪組織與大腦組織距離較近，使得其很難分離，同時獼猴大腦額葉部分比較狹窄，這也在一定程度上增加了腦提取的難度。因此，猴腦提取方法大多還停留在由經(jīng)驗豐富的專家學(xué)者手動提取的水平上，耗時耗力。

目前大腦提取方法分為3類：基于區(qū)域的提取方法、基于圖譜的提取方法和基于深度學(xué)習(xí)的提取方法。基于區(qū)域[2,3]的提取方法依據(jù)目標(biāo)圖像的強度信息進行曲線演化，完成分割。但是，由于只考慮圖像的灰度信息，導(dǎo)致處理存在灰度不均勻問題的醫(yī)學(xué)圖像時不能達(dá)到準(zhǔn)確分割的目的。基于圖譜[4-6]的提取方法是指基于先驗信息制作針對不同獼猴種類的大腦概率圖譜，再通過特定的配準(zhǔn)方法進行大腦提取。這類方法需要投入大量的人力和時間制作不同獼猴種類的概率圖譜，時間開銷較大，無法達(dá)到快速提取的目的。基于深度學(xué)習(xí)[7-9]的提取方法主要是利用貝葉斯卷積神經(jīng)網(wǎng)絡(luò)、Unet算法等深度卷積模型實現(xiàn)對大腦的提取。然而，深度學(xué)習(xí)方法十分依賴于模型的訓(xùn)練和參數(shù)的選擇，需要大量時間和數(shù)據(jù)進行模型調(diào)試。

Figure 1 PCA-Levelset algorithm圖1 PCA-Levelset算法

綜上，面對目前研究成果中存在的問題，本文使用水平集算法來進行猴腦提取。水平集算法不需要大量的數(shù)據(jù)集，且對醫(yī)學(xué)圖像處理效果較好，但是傳統(tǒng)水平集算法對初始輪廓位置的選擇具有隨機性，在一定程度上降低了腦提取的速度，同時由于缺少對邊緣信息的處理算法，進行猴腦提取時不能實現(xiàn)對腦邊緣的精細(xì)提取。

針對上述問題，本文提出了融合分區(qū)與Canny泛函的水平集PCA-Levelset(Level set of fused Partition and CAnny functional)算法。首先，融合分區(qū)的思想，將目標(biāo)圖像分成4個矩形區(qū)域，結(jié)合各區(qū)域的形態(tài)信息在各個區(qū)域選定一個初始輪廓構(gòu)建點，保證構(gòu)建完成的初始輪廓包含較多的大腦組織，這樣可以解決傳統(tǒng)水平集算法中存在的初始輪廓位置隨機選擇的問題。其次，將Canny算子的能量泛函與傳統(tǒng)水平集的LBF(Local Binary Fitting)能量泛函進行融合以提高算法對于邊緣檢測的效果。實驗表明，本文算法進行腦提取的準(zhǔn)確度比傳統(tǒng)水平集算法提高了10%。

2 PCA-Levelset算法

水平集[10]算法是一種求解曲線演化的算法。其思想是：將低維曲線作為高維曲面的零水平集，通過曲線的演化迭代來完成分割。初始輪廓即水平集函數(shù)開始迭代的初始曲線，初始輪廓所包含的目標(biāo)區(qū)域越多，迭代時間越短，結(jié)果越準(zhǔn)確[11]。本文算法共4個步驟(如圖1所示)，包括圖像二值化、刪除小面積連通區(qū)域、分區(qū)以及通過融合Canny算子的能量泛函進行輪廓迭代完成猴腦提取。

2.1 刪除小面積連通區(qū)域

刪除小面積連通區(qū)域可以達(dá)到降低計算量的目的。圖2為對32 127號獼猴腦組織連通區(qū)域的面積統(tǒng)計結(jié)果圖。為了展示小面積連通區(qū)域的統(tǒng)計結(jié)果，圖中已刪除面積大于1 000像素的連通區(qū)域數(shù)據(jù)。從圖2中可以發(fā)現(xiàn)，小面積連通區(qū)域大多集中于200附近，因此設(shè)置Area閾值為200。

Figure 2 Statistics of small area connected area of No.32127 macaque’s brain tissue圖2 32 127號獼猴腦組織小面積連通區(qū)域統(tǒng)計

2.2 分區(qū)

為了使初始輪廓包含較多目標(biāo)區(qū)域，提高實驗結(jié)果的準(zhǔn)確度，本文將圖像I(包含N×M個像素點)分成4個矩形區(qū)域,分別作為構(gòu)建初始輪廓的種子區(qū)域。將目標(biāo)圖像I(N×M)按照式(1)分區(qū)：

I=I1+I2+I3+I4

(1)

其中，I1∈(1:N/2,1:M/2),I2∈(N/2+1:N,1:M/2),I3∈(1:N/2,M/2+1:M),I4∈(N/2+1:N,M/2+1:M)。

圖像I被分成I1,I2,I3,I4共4個矩形區(qū)域。假設(shè)像素點x(i,j)處于(1:N/2,1:M/2)內(nèi)，則將此像素點歸入I1矩形區(qū)域，以此類推。

2.3 融合Canny算子的能量泛函

為克服傳統(tǒng)水平集算法在處理灰度不均勻圖像時對于圖像邊緣檢測不準(zhǔn)確的問題，本文在傳統(tǒng)水平集算法的LBF能量泛函中引入了Canny算子。對于目標(biāo)圖像I，閉合曲線C將圖像域Ω分割為2部分：外部區(qū)域Ω1=outside(C)和內(nèi)部區(qū)域Ω2=inside(C)。

在圖像I中，定義如下LBF能量函數(shù)：

(2)

(3)

為了提高傳統(tǒng)水平集算法對腦邊緣的檢測能力，本文在LBF能量泛函中融合了Canny算子的能量泛函。對于任意一幅圖像I，定義其Canny能量泛函為：

(4)

(5)

增加距離懲罰項函數(shù)[12]保證水平集函數(shù)始終為符號距離函數(shù)，保證演化過程順利完成:

(6)

因此，本文算法得到的能量泛函為：

(7)

其中，μ和v都是常數(shù)。

(8)

固定f1(x)和f2(x)可以得到使能量函數(shù)最小化的梯度下降流為：

(9)

Canny算子對邊緣點的檢測不僅依賴于梯度運算，獨特的方向性使其邊緣檢測和定位的效果優(yōu)于其他算子，具有更好的邊緣強度估計效果，融合到能量泛函中可以準(zhǔn)確檢測腦組織邊緣。

3 實驗結(jié)果及分析

3.1 實驗過程

3.1.1 圖像二值化

使用Otus算法將獼猴大腦圖像進行二值化處理，結(jié)果如圖3所示。

Figure 3 Image binarization圖3 圖像二值化

3.1.2 刪除小面積連通區(qū)域

刪除面積小于200像素的連通區(qū)域后結(jié)果如圖4所示，從圖4中可以看出，獼猴眼睛周圍的小面積肌肉等非腦組織被刪除。

Figure 4 Delete small connected areas圖4 刪除小面積連通區(qū)域

3.1.3 初始輪廓構(gòu)建

如圖5所示，將一幅圖像分成4個矩形區(qū)域，每一個矩形區(qū)域代表一個連通區(qū)域，圓點代表此連通區(qū)域的質(zhì)心點。

Figure 5 Centroids of connected areas in each rectangular area圖5 各矩形區(qū)域中各連通區(qū)域的質(zhì)心

計算分區(qū)點(即點(N/2,M/2),其中,N和M分別為圖像的長和寬)與各區(qū)域質(zhì)心點的歐氏距離，在4個矩形區(qū)域中分別選取歐氏距離最小的質(zhì)心點作為構(gòu)建初始輪廓的種子點，完成初始輪廓的構(gòu)建。圖6中多邊形為構(gòu)建完成的初始輪廓。

Figure 6 Initial contour圖6 初始輪廓

3.2 自對比實驗

將本文算法應(yīng)用于UC-Davis、Princeton和Mountsinai-P 3個不同的數(shù)據(jù)集進行實驗，使用DSC(Dice Similarity Coefficient)和JS(Jaccard Similarity)相似性系數(shù)指標(biāo)來衡量實驗結(jié)果的準(zhǔn)確性，相似性系數(shù)值越接近1，表示相似性越高，結(jié)果越準(zhǔn)確。

實驗中參數(shù)設(shè)置為：標(biāo)準(zhǔn)方差σ=4，μ=1，λ1=λ2=1，v=0.001×255×255，時間步長Δt=0.2，空間步長h=1，Heaviside函數(shù)和Dirac函數(shù)中參數(shù)ε=1。

不同數(shù)據(jù)集上猴腦提取結(jié)果如表1所示。從表1中可以發(fā)現(xiàn)，對于不同的獼猴數(shù)據(jù)集，本文算法都可以實現(xiàn)對腦邊緣的準(zhǔn)確檢測，實現(xiàn)獼猴腦的準(zhǔn)確提取。

Table 1 Comparison of extraction results of our algorithm on different data sets表1 不同數(shù)據(jù)集上本文算法的提取結(jié)果對比

本文算法在不同數(shù)據(jù)集上猴腦提取后的DSC和JS平均值如表2所示。從表2可以發(fā)現(xiàn)，本文算法應(yīng)用于不同數(shù)據(jù)集得到的腦組織與標(biāo)準(zhǔn)腦組織的DSC和JS相似度均可達(dá)到0.74，根據(jù)文獻[13]，DSC值和JS值大于0.7表示算法分割效果較好，本文算法針對不同數(shù)據(jù)集的實驗結(jié)果均達(dá)到此標(biāo)準(zhǔn)，表明本文算法具有較好性能。

3.3 不同算法對比實驗

本節(jié)采用BET算法、分水嶺算法、LBF算法、Unet算法以及本文算法在UC-Davis數(shù)據(jù)集上進行實驗，參數(shù)設(shè)置與3.2節(jié)相同。表3列出了32 127號獼猴第300,320,340的腦提取結(jié)果。

Table 2 Average values of DSC and JS of the brain extracted by our algorithm on different data sets表2 本文算法在不同數(shù)據(jù)集上猴腦提取后的DSC和JS平均值

Table 3 Comparison of extraction results of different algorithms表3 不同算法的提取結(jié)果比較

如表3所示，BET算法、LBF算法和Unet算法無法正確剝離獼猴眼睛周圍的脂肪等非腦組織，并且對于腦邊緣部分提取效果較差。分水嶺算法無法正確提取腦組織。本文算法對邊緣部分提取較準(zhǔn)確，提取效果優(yōu)于其他算法。另外，表3中列出了標(biāo)準(zhǔn)mask圖像，可以發(fā)現(xiàn)本文算法提取結(jié)果與標(biāo)準(zhǔn)mask差異較小。對上述實驗結(jié)果通過DSC和JS指標(biāo)進行定量評估，結(jié)果分別如圖7和圖8所示。

Figure 7 Dice similarity coefficient圖7 DSC相似性系數(shù)

Figure 8 Jaccard similarity coefficient圖8 JS相似性系數(shù)

圖7為DSC相似性系數(shù)的比較結(jié)果。從圖7中可以看出，本文算法相似度集中在0.77左右，最高可到0.86，LBF算法集中在0.65左右，最高可達(dá)0.8,而Watershed集中在0.58左右，最高可達(dá)到0.72，Unet算法集中在0.7左右，最高可達(dá)到0.8，BET算法的DSC集中在0.63左右，最高可到達(dá)0.84。

圖8為Jaccard相似性系數(shù)的比較結(jié)果。從圖8中可以看出，本文算法相似度的JS集中在0.78左右，最高可到0.86，LBF算法集中在0.55左右，最高可達(dá)0.72,而Watershed集中在0.4左右，最高可達(dá)到0.74，Unet算法集中在0.72，最高可達(dá)0.79，BET算法集中在0.62左右，最高可達(dá)到0.87。

綜上，采用本文算法進行腦提取時可以實現(xiàn)獼猴腦組織邊緣的準(zhǔn)確提取，DSC和JS相似性系數(shù)值最高，提取效果最好，且算法穩(wěn)定。

4 結(jié)束語

獼猴大腦組織的自動化提取是獼猴大腦研究中遇到的首要問題，目前對于此問題的解決方法大多是專家手動標(biāo)記，耗時耗力。因此，本文提出了PCA-Levelset算法，首先融合分區(qū)的思想，結(jié)合各區(qū)域形態(tài)信息構(gòu)建初始輪廓，經(jīng)過二值化、刪除小面積連通區(qū)域、分區(qū)域計算質(zhì)心等步驟完成初始輪廓的構(gòu)建，提高了腦提取效率。其次，針對邊緣提取準(zhǔn)確度較低的問題，將Canny算子融合到LBF能量泛函中，提高了本文算法對腦組織邊緣的提取準(zhǔn)確度。大量實驗表明，本文算法可以實現(xiàn)對獼猴腦組織的準(zhǔn)確提取。