李智 陳業(yè)航 馮寶 張紹榮 李昌林 陳相猛 劉壯盛 龍晚生

(1.桂林電子科技大學 電子工程與自動化學院，廣西 桂林 541004；2.桂林航天工業(yè)學院 電子信息與自動化學院， 廣西 桂林 541004；3.中山大學附屬江門市中心醫(yī)院 醫(yī)學影像研究所，廣東 江門 529000)

腦梗死是目前世界范圍內(nèi)致殘率較高的疾病之一，大約有70%～80%的腦梗死患者在接受治療后仍然遺留不同程度的運動功能障礙[1]。腦梗死患者的早期發(fā)現(xiàn)和早期干預能夠明顯地改善患者預后。彌散加權(quán)成像(DWI)是一種檢測人體組織內(nèi)水分子布朗運動情況的核磁共振成像(MRI)序列[2]，它能夠反映早期腦梗死患者腦組織內(nèi)水分子變化情況，對腦梗死的早期臨床診斷具有很高的敏感性，是當前腦梗死患者影像診斷的主要技術(shù)之一。

針對急性期(發(fā)病一周內(nèi)為急性期)的腦梗死患者進行影像評估時，病灶區(qū)的位置和大小與腦梗死患者的功能損傷程度密切相關(guān)，是設(shè)計臨床治療方案的重要參考指標[3-4]。準確分割腦梗死病灶是診斷患者情況的一個重要的預處理步驟。然而，對于急性期腦梗死患者，由于梗死周圍水腫情況復雜多變[5]，使得整個腦梗死病灶在DWI影像上具有邊界模糊、形狀不規(guī)則且內(nèi)部亮度不均勻的特點，給病灶的準確分割帶來了挑戰(zhàn)。因此，如何在DWI影像上準確分割急性期患者的腦梗死病灶，對于影像輔助評估具有重要的臨床意義。

活動輪廓模型(ACM)是目前醫(yī)療影像分割的常用方法之一。它將圖像分割問題轉(zhuǎn)化為水平集函數(shù)求解問題，并可進行拓撲變化，受到了廣泛的關(guān)注[6-8]。張明慧等[9]提出了基于多圖譜的局部區(qū)域活動輪廓模型，該方法將多圖譜的形狀先驗信息引入活動輪廓模型中實現(xiàn)了對腦結(jié)構(gòu)的自動精確分割。欒紅霞等[10]利用目標輪廓的先驗知識對目標進行邊緣約束，然后將邊緣約束后的邊緣吸引力場進行正則化并融入邊界活動輪廓模型中，可以實現(xiàn)邊界模糊的腦MRI圖像中腦組織的分割。上述兩種方法都是基于目標先驗形狀信息的活動輪廓模型，但是腦梗死的發(fā)病部位和形狀都是隨機的，很難找到合適的先驗形狀信息。Sandhya等[11]提出一種自組織映射網(wǎng)絡(luò)與區(qū)域活動輪廓模型相結(jié)合的分割方法，通過自組織映射網(wǎng)絡(luò)訓練神經(jīng)元對圖像灰度分布進行全局建模，進而獲得圖像灰度分布的拓撲變化，然后將訓練好的神經(jīng)元引入活動輪廓模型中，驅(qū)動輪廓曲線的演化，但自組織映射的參數(shù)(如學習率、學習終止)需要人為控制，且神經(jīng)元的數(shù)目難以確定。馬超等[12]提出了一種整合級聯(lián)隨機森林與活動輪廓模型的腦MR圖像分割方法，該方法通過級聯(lián)隨機森林對演化曲線進行初始化，并利用活動輪廓模型進一步演化曲線。Ma等[13]提出基于隨機森林的活動輪廓模型，首先對隨機森林進行訓練，預測出待分割圖像的腦腫瘤結(jié)構(gòu)，然后將其作為局部活動輪廓模型的初始輪廓和空間先驗信息進行細分割。但由于腦梗死的DWI圖像噪聲較大，且發(fā)生腦梗死后，其水腫區(qū)域會以小時為單位不斷變化，隨機森林模型容易陷入過擬合。王醒策等[14]提出了一種結(jié)合圖像全局信息和局部信息的活動輪廓模型，使得該模型可以有效分割亮度不均勻的腦血管圖像。Charoensuk等[15]利用局部區(qū)域活動輪廓模型在DWI圖像上對腦梗死病灶進行分割。唐文杰等[16]提出一種雙水平集方法對腦部MR圖像進行分割。上述文獻[14]至文獻[16]都是基于局部區(qū)域的活動輪廓模型，原理是根據(jù)圖像的局部灰度統(tǒng)計信息驅(qū)動輪廓曲線的演變，可較好的分割亮度不均勻圖像，但由于腦梗死核心周圍有水腫，致使DWI圖像中梗死區(qū)域出現(xiàn)邊緣模糊的特點，基于局部灰度統(tǒng)計信息的活動輪廓模型無法準確收斂至目標邊界。近年來，一些學者也將深度學習應用到腦MRI圖像分割。Moeskops等[17]提出了一種多尺度的分段卷積神經(jīng)網(wǎng)絡(luò)來分割嬰兒和年輕人的大腦組織。Dou等[18]提出一種3D卷積神經(jīng)網(wǎng)絡(luò)來分割腦MRI圖像中腦出血病灶。但由于醫(yī)學領(lǐng)域的敏感性和特殊性，較少的訓練數(shù)據(jù)限制了深度學習方法的性能，深度學習模型龐大的模型規(guī)模所具有的優(yōu)勢此時并不明顯。

文中針對急性期腦梗死病灶在DWI圖像中存在的邊界模糊、形狀不規(guī)則和亮度不均勻等問題，提出一種模糊速度函數(shù)驅(qū)動下的活動輪廓模型，來完成腦梗死病灶的準確分割：(1)在模型初始化方面，利用小波能量圖下的貝葉斯概率獲取模型初始輪廓，使得初始輪廓快速定位于腦梗死病灶的真實邊界附近，增強模型的魯棒性和準確性。(2)將圖像局部熵引入活動輪廓模型中。圖像局部熵反映了該局部圖像亮度信息的有序性，可表征腦梗死病灶和正常腦組織之間水分子分布的差異性，在一定程度上解決圖像的亮度不均勻性和噪聲的干擾問題。(3)根據(jù)腦梗死病灶區(qū)的邊界模糊特性，提出結(jié)合圖像局部熵和灰度的模糊聚類算法計算模糊隸屬度，進一步加強腦梗死病灶與正常組織的區(qū)分。(4)將基于模糊隸屬度的模糊速度函數(shù)引入到模型中，使輪廓曲線在腦梗死病灶的模糊邊界處停止演變，完成腦梗死病灶的分割。(5)最后利用臨床數(shù)據(jù)驗證本文方法的有效性。

1 模糊速度函數(shù)驅(qū)動下的活動輪廓模型

如圖1所示，文中提出的分割算法主要包括4個步驟：1)在小波變換域下計算圖像的貝葉斯概率，進行活動輪廓模型的輪廓曲線初始化，并得到基于初始輪廓的局部區(qū)域；2)結(jié)合圖像局部熵和灰度的模糊聚類算法計算模糊隸屬度；3)將基于模糊隸屬度的模糊速度函數(shù)引入到活動輪廓模型中，構(gòu)建基于模糊速度函數(shù)的活動輪廓模型；4)活動輪廓模型能量泛函演化。

1.1 基于貝葉斯概率的活動輪廓模型的初始化

模型輪廓曲線的初始位置對于模型的分割精度有很大的影響。為了解決對初始輪廓敏感的問題，在小波域下構(gòu)建貝葉斯概率對模型進行初始化。首先，計算圖像的小波能量，小波能量可以增強腦梗死病灶區(qū)與正常組織的區(qū)分。然后，在小波能量圖上計算每個像素的貝葉斯概率，獲取活動輪廓模型的初始輪廓，使得初始輪廓位于腦梗死病灶的真實邊界附近，提高分割的準確性。

對DWI圖像的像素(x，y)進行二維離散小波變換后，小波能量可以定義為[19]

W(x，y)=EA(x，y)+EH(x，y)+EV(x，y)+

ED(x，y)

(1)

式中：EA(x，y)、EH(x，y)、EV(x，y)和ED(x，y)分別為低頻小波能量和3個高頻小波能量，定義如下：

圖1 文中算法流程圖

(2)

式中：DA(x，y)為低頻小波系數(shù)；DH(x，y)、DV(x，y)和DD(x，y)分別為水平、垂直、對角線方向的高頻小波系數(shù)；Q是整數(shù)的集合；r和s分別為鄰域Q內(nèi)的像素的橫坐標和縱坐標；K(·)為高斯核函數(shù)。采用小波框架方法，將小波頻域下的圖像大小擴展到原始圖像的大小。

假設(shè)小波能量圖符合混合高斯模型[22]，則將小波能量圖像分為腦梗死區(qū)域(r=1)與背景區(qū)域(r=2)兩部分?；旌细咚鼓Ｐ涂杀硎緸?/p>

(3)

L(θ)=logP(w|θ)

(4)

(5)

則模型的初始輪廓可定義為

(6)

圖2顯示了小波域下采用貝葉斯概率獲取的初始輪廓。圖2(a)是腦梗死病灶區(qū)的DWI原圖，圖2(b)是小波能量分布圖，圖2(c)是小波域下基于貝葉斯概率獲得的初始輪廓，可看出該初始輪廓位于梗死病灶真實邊界附近。圖2(d)是水平集函數(shù)的初始化圖，中間白色曲線為水平集初始化曲線對應的輪廓曲線。

圖2 基于貝葉斯概率初始輪廓

1.2 計算圖像模糊隸屬度

文中根據(jù)腦梗死病灶邊界的模糊特性計算模糊速度函數(shù)，即采用結(jié)合局部熵和灰度的模糊聚類算法，在初始輪廓確定的局部區(qū)域中，計算模糊隸屬度，從而進一步加強腦梗死病灶和正常組織的區(qū)分。

由于急性腦梗死患者梗死中心區(qū)的神經(jīng)細胞已經(jīng)壞死，無法正常新陳代謝，且周圍伴隨缺血半暗帶和水腫情況的出現(xiàn)，因此梗死區(qū)的水分子分布失去固定規(guī)律，與正常腦組織有顯著差異。而圖像的局部熵反映了該局部圖像亮度信息的有序性，可表征腦梗死病灶和正常腦組織之間水分子分布的差異性，在一定程度上解決圖像的亮度不均勻性和噪聲的干擾問題，從而加強腦梗死病灶和正常腦組織的區(qū)分。對于一幅大小為M×N的圖像，設(shè)I(x，y)為圖像中像素點(x，y)處的灰度值，則局部熵可定義為[20]

(7)

其中，

(8)

式中，pxy為像素點(x，y)處的灰度分布概率，如果M×N是以(x，y)為中心的一個局部圖像窗口，則H為圖像的局部熵。文中的局部圖像窗口大小設(shè)置為9×9。從圖3(e)和圖3(f)可以看出正常組織的局部熵值小于腦梗死病灶局部熵值，這與腦梗死病灶的水分子分布失去固定規(guī)律，比正常組織具有更多信息的臨床背景相吻合。但從圖3(c)、圖3(d)、圖3(e)和圖3(f)可以看出，單一使用圖像的灰度和圖像局部熵很難將腦梗死病灶從正常組織中分離出來。所以考慮構(gòu)建結(jié)合局部熵特征和灰度特征的模糊聚類算法，計算圖像的模糊隸屬度，進一步加強腦梗死病灶和正常腦組織的區(qū)分。

模糊聚類算法通過目標函數(shù)迭代尋找最優(yōu)參數(shù)，目標函數(shù)定義為

(9)

式中：n表示圖片像素點的總數(shù)；l表示第l個像素點，Ol為圖像局部熵和灰度構(gòu)成的特征向量；U為2×n的隸屬度矩陣；ulh是向量Ol隸屬h類的程度；G為聚類中心向量；m為一個模糊控制常數(shù)；d(Ol，Gh)為Ol和距離中心Gh之間的距離，該距離函數(shù)定義為

d(Ol，Gh)=‖Ol-Gh‖

(10)

重復計算ulh和Gh，當所有的特征向量被正確分類時，J(U，G)目標函數(shù)達到最小值。ulh和Gh的計算公式如下：

(11)

(12)

根據(jù)ulh構(gòu)造隸屬腦梗死病灶的模糊隸屬度矩陣Z，則像素點(x，y)的模糊隸屬度可表示為Z(x，y)∈[0，1]。從圖3(g)和圖3(h)中可以看出，文中結(jié)合灰度值和局部熵的模糊聚類算法，將基于灰度的圖像空間區(qū)域轉(zhuǎn)變到基于模糊隸屬度的圖像空間域，使得腦梗死病灶和正常組織的模糊隸屬度以0.5為界，從而加強了腦梗死病灶與正常組織的區(qū)分。

(a)原始圖像 (b)病灶區(qū)域放大圖

(c)正常組織的灰度值 (d)病灶組織的灰度值

(e)正常組織的局部熵值 (f)病灶組織的局部熵值

(g)正常組織的模糊隸屬度(h)病灶組織的模糊隸屬度

Fig.3 Gray value，local entropy and fuzzy membership of DWI image

1.3 構(gòu)建模糊速度函數(shù)

為了使模型在腦梗死病灶邊界處停止演變，根據(jù)模糊隸屬度矩陣，分別定義使輪廓曲線向內(nèi)演變的模糊速度函數(shù)V1、輪廓曲線向外演變的模糊速度函數(shù)V2和反映輪廓邊界的模糊速度函數(shù)V3，并有：1)當輪廓曲線遠離邊界時，V1、V2和V3越來越大，驅(qū)動輪廓曲線快速演變；2)當輪廓曲線趨近邊界時，V1、V2和V3逐漸變小，且在目標邊界處時V1≈V2≈V3≈0，輪廓曲線停止演變。

文中所提出的模型中，V1、V2和V3分別定義為：

V1(x，y)=|e-t1[Z(x，y)-0.5]-1|

(13)

V2(x，y)=|e-t2[0.5-Z(x，y)]-1|

(14)

V3(x，y)=et3[Z(x，y)-0.5]2-1

(15)

式中，t1、t2和t3為整數(shù)，Z(x，y)為像素點(x，y)隸屬腦梗死病灶的模糊隸屬度。

圖4展示了一個腦梗死病灶的模糊速度函數(shù)。從圖4(b)和圖4(c)中可以看出，在腦梗死病灶邊界點A、點B、點C和點D處，模糊速度函數(shù)Vk(A)、Vk(B)、Vk(C)和Vk(D)，k∈(1，2，3)分別趨近0。

(a)原始圖像

(b)實型直線的模糊速度

(c)點型直線的模糊速度

1.4 構(gòu)建活動輪廓模型的能量泛函

假設(shè)Ω是DWI圖像的定義域，輪廓曲線C將Ω分為輪廓曲線的內(nèi)部區(qū)域Ω1和外部區(qū)域Ω2。根據(jù)模糊速度函數(shù)V1和V2可定義活動輪廓模型的區(qū)域項：

V2(x，y)|I(x，y)-c2|2(1-H(φ))]dxdy

(16)

式中，φ為水平集函數(shù)，I(x，y)為像素點(x，y)處的灰度值，c1和c2分別表示輪廓曲線內(nèi)部灰度平均值和輪廓曲線外部灰度平均值。H(φ)為正則化Heaviside函數(shù)，參數(shù)ε用以控制δ(φ)的有效寬度。

根據(jù)反映輪廓邊界的模糊速度函數(shù)V3可定義活動輪廓模型的邊界檢測項：

(17)

結(jié)合式(16)和式(17)，活動輪廓模型的能量泛函可定義為

V2(x，y)|I(x，y)-c2|2(1-H(φ))]dxdy+

(18)

式中：等號右邊第1項是控制輪廓曲線向目標邊界演變的區(qū)域項；第2項是控制輪廓曲線在目標邊界處停止演化的邊界檢測項；第3項是為了避免輪廓曲線在演變過程中對水平集函數(shù)φ進行重新初始化[21]。μ、λ和γ分別是控制區(qū)域項、邊界檢測項和正則項的權(quán)重參數(shù)，在文中μ、λ、γ和t1、t2、t3的取值都為1[21]。

1.5 能量泛函演化

活動輪廓模型輪廓曲線演化過程可歸結(jié)為最小化能量泛函E(φ)。首先，固定水平集函數(shù)φ，相對c1和c2最小化E(φ)，可得到：

(19)

為了簡化表示，書寫上省略變量(x，y)，E(φ)可簡化為

F(φ，φx，φy)=λe1H(φ)+λe2(1-H(φ))+

(20)

然后，固定c1和c2，最小化E(φ)對應于求解如下偏微分方程：

(21)

其中F對φ、φx、φy的偏導數(shù)可以表示為

(22)

(23)

(24)

進一步對式(23)和(24)求導：

(25)

(26)

同時結(jié)合

(27)

將式(22)、(25)-(27)代入式(21)，可以得到F關(guān)于水平集函數(shù)φ的歐拉-拉格朗日方程：

(28)

由變分原理，可得φ的梯度下降流：

(29)

2 分割結(jié)果評價指標

在評估算法有效性的過程中，分別定義陽性率(RTP)、假陽性率(RFP)和相似度(DS)3個指標來評價算法的分割結(jié)果。

(30)

式中:Am為醫(yī)生手動分割結(jié)果；Aa為分割算法結(jié)果；RTP的值越大，說明Aa對Am的覆蓋率越好；RFP的值越小，說明Aa包含越少的背景；DS的值越大，說明Aa與Am越相似。

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

實驗數(shù)據(jù)采用江門市中心醫(yī)院的臨床數(shù)據(jù)。由影像科醫(yī)生從臨床數(shù)據(jù)中隨機選擇11 個邊界模糊、亮度不均勻的腦梗死患者數(shù)據(jù)進行分割結(jié)果的對比展示。將文中提出的方法與基于區(qū)域的ACM、基于邊界的ACM、結(jié)合二值水平集和形態(tài)學的ACM、基于局部高斯模型的ACM、基于模糊水平集的ACM和基于圖像局部熵的ACM[20]進行分割結(jié)果的對比。

分割結(jié)果如圖5所示，藍色曲線包圍的區(qū)域為分割結(jié)果。圖5(d)和圖5(f)分別展示了基于區(qū)域的ACM和結(jié)合二值水平集和形態(tài)學的ACM的分割結(jié)果，可以看出在圖像邊界模糊、亮度不均勻的情況下這兩種方法的分割結(jié)果并不理想，如第K個數(shù)據(jù)分割結(jié)果出現(xiàn)了邊界泄漏。圖5(e)展示了基于邊界的ACM的分割結(jié)果，在模糊邊界處輪廓曲線無法收斂到目標邊界且在凹陷邊界處都存在邊界泄漏的問題，見第A、C、I和J個數(shù)據(jù)。圖5(g)是基于局部高斯模型的ACM的分割結(jié)果，第A、J、H、I和K個數(shù)據(jù)均發(fā)生了邊界泄漏問題。圖5(h)展示了基于模糊水平集的ACM分割結(jié)果，可以看到在多病灶的情況下分割不完整(如第G、H和K個數(shù)據(jù))，且在邊界模糊和亮度不均勻的情況下存在著過分割(如第J個數(shù)據(jù))。圖5(i)是基于圖像局部熵的分割結(jié)果，當腦梗死病灶存在邊界模糊的問題時，該方法不能準確的分割腦梗死病灶，如第A，F(xiàn)和G個數(shù)據(jù)。圖5(j)是文中提出方法的分割結(jié)果，較其他幾種分割算法更接近醫(yī)生手動分割結(jié)果。

圖5 腦梗死病灶分割結(jié)果

對上述分割結(jié)果分別計算RTP、RFP和DS3個指標，結(jié)果如表1所示。結(jié)合二值水平集和形態(tài)學的ACM雖然得到最好的DS值，但存在欠分割問題；基于邊界的ACM雖然得到最好的RTP值，但存在過分割問題；而文中方法DS為0.887，RTP為0.908，RFP為0.025，說明文中方法的結(jié)果與醫(yī)生手動分割的結(jié)果最為相近。

表1 11個臨床病例的分割結(jié)果

為了進一步驗證本文方法的準確性，又隨機選擇58個病例進行分割實驗，結(jié)果如表2所示。綜合3個指標，可以看出文中方法的準確性更好。

從以上結(jié)果可知，由于腦梗死病灶存在模糊特性和對比度低等問題，基于邊界的、區(qū)域的、結(jié)合二值水平集和形態(tài)學的、基于局部高斯模型的和基于模糊水平集的ACM無法準確分割病灶。同樣的，基于圖像局部熵特征的ACM是一個結(jié)合全局信息和局部信息的ACM，該方法利用全局或局部的灰度統(tǒng)計信息驅(qū)動活動輪廓模型進行演變，但由于腦梗死病灶具有邊緣模糊的特點，基于全局或局部的灰度統(tǒng)計信息無法有效處理邊界模糊的問題。

表2 58個臨床病例的分割結(jié)果

文中方法利用結(jié)合局部熵特征和灰度特征的模糊聚類算法計算圖像的模糊隸屬度，加強腦梗死病灶與正常組織的區(qū)分。然后分別計算3個基于模糊隸屬度的模糊速度函數(shù)構(gòu)建模型，使得模型在腦梗死病灶邊界處能夠停止演變。最后選取一個腦梗死病例的所有DWI圖像進行分割，該病灶有8個連續(xù)層面。結(jié)果如圖6所示，對比圖6(b)和圖6(c)可以看出，文中提出方法的分割結(jié)果與醫(yī)生手動分割結(jié)果相似。

圖6 文中提出的方法對腦梗死病灶DWI序列圖像分割

4 結(jié)論

為了準確分割DWI圖像中的腦梗死病灶，文中提出用模糊速度函數(shù)驅(qū)動下的活動輪廓模型對其進行分割。該模型一方面提出結(jié)合圖像局部熵特征和灰度特征的模糊聚類算法計算模糊隸屬度，以加強腦梗死病灶與正常組織的區(qū)分，并利用基于模糊隸屬度的模糊速度函數(shù)構(gòu)建模型；另一方面利用小波變換域下的貝葉斯概率獲取模型初始輪廓，增強模型的魯棒性和準確性。實驗結(jié)果表明，本文的模型可有效分割DWI圖像中的腦梗死病灶。