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

郭晉秀，張?jiān)路?，鄧紅霞，李海芳

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

1 引言

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

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

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

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

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

2 PCA-Levelset算法

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

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

刪除小面積連通區(qū)域可以達(dá)到降低計(jì)算量的目的。圖2為對(duì)32 127號(hào)獼猴腦組織連通區(qū)域的面積統(tǒng)計(jì)結(jié)果圖。為了展示小面積連通區(qū)域的統(tǒng)計(jì)結(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號(hào)獼猴腦組織小面積連通區(qū)域統(tǒng)計(jì)

2.2 分區(qū)

為了使初始輪廓包含較多目標(biāo)區(qū)域，提高實(shí)驗(yàn)結(jié)果的準(zhǔn)確度，本文將圖像I(包含N×M個(gè)像素點(diǎn))分成4個(gè)矩形區(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個(gè)矩形區(qū)域。假設(shè)像素點(diǎn)x(i,j)處于(1:N/2,1:M/2)內(nèi)，則將此像素點(diǎn)歸入I1矩形區(qū)域，以此類推。

2.3 融合Canny算子的能量泛函

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

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

(2)

(3)

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

(4)

(5)

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

(6)

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

(7)

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

(8)

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

(9)

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

3 實(shí)驗(yàn)結(jié)果及分析

3.1 實(shí)驗(yàn)過(guò)程

3.1.1 圖像二值化

使用Otus算法將獼猴大腦圖像進(jìn)行二值化處理，結(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個(gè)矩形區(qū)域，每一個(gè)矩形區(qū)域代表一個(gè)連通區(qū)域，圓點(diǎn)代表此連通區(qū)域的質(zhì)心點(diǎn)。

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

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

Figure 6 Initial contour圖6 初始輪廓

3.2 自對(duì)比實(shí)驗(yàn)

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

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

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

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

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

3.3 不同算法對(duì)比實(shí)驗(yàn)

本節(jié)采用BET算法、分水嶺算法、LBF算法、Unet算法以及本文算法在UC-Davis數(shù)據(jù)集上進(jìn)行實(shí)驗(yàn)，參數(shù)設(shè)置與3.2節(jié)相同。表3列出了32 127號(hào)獼猴第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算法無(wú)法正確剝離獼猴眼睛周圍的脂肪等非腦組織，并且對(duì)于腦邊緣部分提取效果較差。分水嶺算法無(wú)法正確提取腦組織。本文算法對(duì)邊緣部分提取較準(zhǔn)確，提取效果優(yōu)于其他算法。另外，表3中列出了標(biāo)準(zhǔn)mask圖像，可以發(fā)現(xiàn)本文算法提取結(jié)果與標(biāo)準(zhǔn)mask差異較小。對(duì)上述實(shí)驗(yàn)結(jié)果通過(guò)DSC和JS指標(biāo)進(jìn)行定量評(píng)估，結(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。

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

4 結(jié)束語(yǔ)

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