鄭 偉，李 涵，安曉林，劉帥奇，張曉丹，馬澤鵬*

(1.河北大學(xué) 電子信息工程學(xué)院， 保定 071002；2.河北省數(shù)字醫(yī)療工程重點(diǎn)實(shí)驗(yàn)室，保定 071002；3.河北省機(jī)器視覺工程技術(shù)研究中心，保定071002；4.河北大學(xué)附屬醫(yī)院，保定 071000)

引 言

中國已進(jìn)入老齡化階段，預(yù)計老齡化進(jìn)程在未來十幾年內(nèi)將不斷加速[1]。阿爾茨海默病(Alzheimer disease,AD)在老年群體中發(fā)病率不斷增加，并給社會、家庭帶來不可估計的經(jīng)濟(jì)壓力[2]。影像學(xué)是AD診斷的重要方式,包括計算機(jī)斷層掃描技術(shù)(computed tomography，CT)、核磁共振成像(magnetic resonance imaging,MRI)、正電子發(fā)射斷層成像(positron emission tomography,PET)等，近年來PET/CT,PET/MRI兩種影像結(jié)合的形式也成為重要的診斷方式[3]。其中CT檢查中的X射線對活體組織有一定損傷，且對血管、神經(jīng)等軟組織分辨率有限。MRI 的磁場與人體的水分子包含的氫原子核的旋轉(zhuǎn)作用釋放能量產(chǎn)生非常清晰的解剖結(jié)構(gòu)。醫(yī)生通過AD患者的MRI腦影像清晰地看到腦溝回加深加寬、腦室擴(kuò)大程度、海馬萎縮等，但由于很多疾病可以導(dǎo)致腦萎縮，僅僅依靠MRI幾乎不能確診。PET通過在人體中注射放射性物質(zhì)，可以無創(chuàng)生理地研究人體內(nèi)生物分子狀態(tài)。AD患者表現(xiàn)為神經(jīng)突觸的丟失或機(jī)能的降低，PET示蹤劑18氟-氟脫氧葡萄糖(18F-fluorodeoxyglucose，18F-FDG)通過葡萄糖代謝高低反映神經(jīng)突觸的活性，表現(xiàn)在皮層區(qū)域葡萄糖代謝下降，但圖像空間分辨率低。因此,采用圖像融合技術(shù)生成PET/MRI圖像，利用MRI清晰的結(jié)構(gòu)精準(zhǔn)判斷PET葡萄糖代謝低的位置，提高AD的特異性和準(zhǔn)確性。

PET/MRI圖像融合可在空域和變換域中實(shí)現(xiàn)，其中空域算法如加權(quán)平均、像素取大等原理簡單、且耗時短，但忽略了紋理等信息?；谧儞Q域的多尺度幾何變換[4-5]的結(jié)果與人眼視覺效果相近被廣泛使用。ZHANG等人[6]提出基于小波變換的PET和MRI圖像融合，小波不能很好地表達(dá)具有方向特征的線、輪廓等，導(dǎo)致融合結(jié)果出現(xiàn)塊效應(yīng)。因此提出了具有更好表達(dá)方向信息的變換，如輪廓波[7]、雙樹小波[8]、剪切波變換[9]。由于剪切波同時具有數(shù)學(xué)表達(dá)式簡單、通過一個或幾個函數(shù)仿射變換產(chǎn)生、最佳表達(dá)高維信號的各向異性、可實(shí)現(xiàn)空域緊支、可擴(kuò)展到3維甚至更高維的優(yōu)點(diǎn)是其它變換不能同時擁有的，掀起了廣大學(xué)者的研究熱潮。LABATE等人[10-11]提出3維形式的剪切波變換。DUAN等人[12]基于3-D帶限剪切波變換的磁共振圖像融合，證明了3維剪切波彌補(bǔ)了小波變換對于方向特征的不足，但高低頻所用的融合規(guī)則未顧及到體素點(diǎn)所在位置結(jié)構(gòu)特征。WANG等人[13]基于3-D剪切波變換分解并用廣義高斯密度建模，不同子帶間的相關(guān)性用相對熵統(tǒng)計，產(chǎn)生較大的冗余運(yùn)行速度較慢。YANG[14]基于3維醫(yī)學(xué)圖像體素點(diǎn)的物理特性和3-D剪切波融合，但是未顧及低頻子帶中的細(xì)節(jié)信息，使融合結(jié)果缺少細(xì)節(jié)特征。

為此,本文中充分考慮PET和MRI圖像的結(jié)構(gòu)特征和相同位置強(qiáng)度差異較大問題，提出了一種新的PET/MRI醫(yī)學(xué)圖像融合方法。采用ShearLab 3D將PET和MRI圖像分解，并以方差為閾值獲得高頻子帶中的強(qiáng)邊緣子帶，融合規(guī)則充分考慮到3維體素的空間位置信息，將低中高頻融合規(guī)則由2維擴(kuò)展到3維形式。低頻子帶以銳化矩陣增強(qiáng)圖像邊緣信息，并以3維擴(kuò)展的加權(quán)局部能量與拉普拉斯算子加權(quán)和保留圖像能量和細(xì)節(jié)信息。中頻子帶以絕對值為活動度量捕捉邊緣信息。高頻子帶結(jié)合3維擴(kuò)展的3個底層特征加權(quán)增強(qiáng)圖像的細(xì)節(jié)特征提高融合圖像亮度和對比度。結(jié)果表明，與空域算法、小波變換、不同的融合規(guī)則相比本文中的算法融合圖像邊緣清晰、細(xì)節(jié)豐富符合人眼視覺效果，更適用于PET/MRI圖像融合。

1 ShearLab 3D理論

3-D剪切波系統(tǒng)是由各向異性膨脹矩陣和剪切矩陣組成[10]，它們分別控制尺度和不同尺度的方向，設(shè)各向異性膨脹矩陣為:

(1)

剪切矩陣為:

(2)

式中,j∈Z為尺度參量，α∈(0,2)為各向異性程度，k=(k1,k2)∈Z2為剪切方向，Z為整數(shù)集。移位采樣晶格為下式:

(3)

式中，c1>0,c2>0為平移量。

(4)

式中，

Φ(φ;c1)={φm=φ(·-c1m),m∈Z2)}

(5)

(j≥0,|k|≤「2j(αj-1)/2?,m∈Z2)}

(6)

(j≥0,|k|≤「2j(αj-1)/2?,m∈Z2)}

(7)

(j≥0,|k|≤「2j(αj-1)/2?,m∈Z2)}

(8)

(9)

式中，函數(shù)φ1和ψ1是緊密支持的函數(shù)，2維扇形濾波函數(shù)P都滿足一定條件[11]。

2 PET/MRI圖像融合

首先從ADNI(Alzheimer’s Disease Neuroimaging Initiative)數(shù)據(jù)庫中下載AD患者的PET和MRI腦功能影像，分別在統(tǒng)計參量圖(statistical parametric mapping,SPM)中對源圖像經(jīng)過頭動校正、配準(zhǔn)、歸一化處理得到兩個待處理對象圖像A,B。然后利用ShearLab 3D分解分別得到低頻子帶Cl,A(x,y,z)和Cl,B(x,y,z)高頻子帶Ch,0,A(x,y,z)，Ch,0,B(x,y,z),其中高頻子帶以方差為閾值獲得強(qiáng)邊緣子帶為中頻子帶Cm,A(x,y,z)和留下的高頻子帶Ch,1,A(x,y,z),Ch,1,B(x,y,z)。本文中提出將低中高頻融合規(guī)則擴(kuò)展為3維，以加權(quán)局部能量與拉普拉斯算子加權(quán)和融合低頻系數(shù)并以銳化矩陣作為權(quán)重參量捕捉圖像中灰度跳變區(qū)域，以絕對值為活動度量融合中頻子帶，以3個3維擴(kuò)展的底層特征加權(quán)融合高頻系數(shù)。最后，經(jīng)過ShearLab 3D逆變換得到融合結(jié)果PET/ MRI圖像F，流程如圖1所示。

Fig.1 PET/MRI image fusion flowchart

2.1 低頻子帶融合規(guī)則

ShearLab 3D分解生成的低頻子帶不但包含圖像絕大多數(shù)的能量，也含有圖像的紋理特征。因此低頻子帶融合規(guī)則的選擇會影響融合結(jié)果的好壞，常用的加權(quán)平均法使融合結(jié)果輪廓不完整清晰、對比度降低。像素取大法的融合結(jié)果中MRI圖像的腦溝回等細(xì)節(jié)特征不足。

本文中考慮到每個體素點(diǎn)的空間信息將2維局部拉普拉斯能量加權(quán)[15]擴(kuò)展到3維形式。設(shè)EWLE,s(x,y,z)為低頻子帶位置(x,y,z)的加權(quán)局部能量(weight local energy，WLE)，s∈{A,B}分別表示待融合的PET和MRI圖像。為了突出了MRI圖像邊緣信息，本文中將W設(shè)置為(2r+1)×(2r+1)×(2r+1)的銳化矩陣，r是W的半徑，即:

Cl,s(x+m,y+n,z+p)2

(10)

(11)

W3×3×3=

(12)

設(shè)LWSTML,s(x,y,z)是低頻子帶位置(x,y,z)處二十六鄰域的改進(jìn)拉普拉斯算子的加權(quán)和(weight sum of twenty six-neighborhood based modified laplacian，WSTML)。LTML,s[15]為改進(jìn)的拉普拉斯算子(twenty six-neighborhood based modified laplacian，TML)，將相鄰點(diǎn)的權(quán)重設(shè)為歐氏距離。為了充分利用相鄰體素間的關(guān)系,將LTML,s擴(kuò)展到3維，即:

LTML,s(x+m,y+n,z+p)2

(13)

[y-(y+g)]2+[z-(z+e)]2}-1/2×

|Cl,s(x,y,z)-Cl,s(x+i,y+q,z+e)|

(14)

式中,i，q，e都為整數(shù)且不能同時為0；s,W,r同上。

位置(x,y,z)的低頻融合系數(shù)Fl(x,y,z)為下式，選取局部拉普拉斯能量大的系數(shù)為融合圖像的低頻系數(shù):

Fl(x,y,z)=

(15)

2.2 中頻子帶融合規(guī)則

ShearLab 3D分解后高頻子帶包含圖像的邊緣和紋理結(jié)構(gòu)，為了更好地融合圖像的邊緣信息在高頻子帶中提取包含源圖像強(qiáng)邊緣信息中頻子帶。由于方差可以判斷圖像的紋理、邊緣、平滑區(qū)域，其中紋理區(qū)方差最大，邊緣次之，平滑最小[16]。因此以方差為閾值得到中頻子帶，且方差的選取與圖像特點(diǎn)有關(guān)。

(16)

中頻子帶采取基于絕對值的活動度量[17]保留中頻子帶中的邊緣信息，并充分考慮到空間結(jié)構(gòu)信息將其擴(kuò)展到3維情況。

設(shè)Cm,A(x,y,z)與Cm,B(x,y,z)為PET,MRI圖像的中頻子帶，以下式二值化處理：

G(x,y,z)=

(17)

式中,G(x,y,z)為二值化矩陣，當(dāng)PET圖像的中頻子帶系數(shù)絕對值大于MRI圖像的中頻子帶系數(shù)絕對值時,G(x,y,z)=1，否則為0。

d(x,y,z)為中頻子帶融合的決策圖，g(x,y,z)為活動指標(biāo)，通過以下等式計算:

(18)

g(x,y,z)={(x0,y0,z0)∈Ω|G(x,y,z)|}

(19)

式中,Ω是以(x,y,z)為中心的M×N×P的滑動窗口。

中頻子帶的融合系數(shù)為Fm(x,y,z)，當(dāng)決策圖d(x,y,z)=1時,PET圖像的中頻系數(shù)為融合圖像的中頻子帶系數(shù)。

(20)

2.3 高頻子帶融合規(guī)則

高頻子帶包含了MRI源圖像腦溝回等細(xì)節(jié)特征，對AD病情的診斷非常重要。因此高頻子帶融合規(guī)則應(yīng)有效的保留細(xì)節(jié)信息，局部相位一致的大部分特征為邊緣狀和角狀特征[18]。局部對比度、局部能量可以反映對比度和亮度信息。本文中采取這3種底層特征加權(quán)[19]的融合規(guī)則，考慮到每個體素點(diǎn)的空間信息將其擴(kuò)展到3維形式。

設(shè)PPC,s(x,y,z)為高頻子帶(x,y,z)處的相位一致(phase congruency，PC)的點(diǎn)，s∈(A,B)，即:

(21)

式中,θk是k的方向角，An,θk為第n個傅里葉分量的幅值和角度，ε=0.001是正數(shù)去除圖像中的DC分量，Eθk(x,y,z)可見參考文獻(xiàn)[19]。

SLSCM,s(x,y,z)為(x,y,z)處的局部對比度(local measure of sharpness change,LSCM)，(2M+1)×(2N+1)×(2P+1)是鄰域大小，SSCM,s(x,y,z)為高頻子帶(x,y,z)處的對比度(measure of sharpness change,SCM)的值。

y+b,z+t)

(22)

Ch,1,s(x0,y0,z0)]

(23)

式中,取M=N=P=1，Ω0是(x,y,z)的3×3×3的局部區(qū)域，(x0,y0,z0)表示Ω0內(nèi)的體素點(diǎn)。

ELE,s(x,y,z)為高頻子帶(x,y,z)處的局部能量(local energy,LE)，且M=N=P=1。

ELE,s(x,y,z)=

(24)

用計算出來的PPC,s，SLSCM,s，ELE,s組成活躍度測量Ns;α1，β1和γ1是用來調(diào)整PPC,s，SLSCM,s，ELE,s的參量，分別為1,2,2[19]:

Ns(x,y,z)=[PPC,s(x,y,z)]α1·

[SLSCM,s(x,y,z)]β1·[ELE,s(x,y,z)]γ1

(25)

d(x,y,z)為高頻子帶融合的決策圖，通過(18)式計算，其中g(shù)(x,y,z)由等下式計算:

g(x,y,z)={(x0,y0,z0)∈Ω1|Ns≥

max(NA(x0,y0,z0),NB(x0,y0,z0))}

(26)

高頻子帶的融合系數(shù)Fh(x,y,z)與(20)式相同，當(dāng)決策圖d(x,y,z)=1時,PET圖像系數(shù)為融合圖像的高頻系數(shù)。

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

本文中所有實(shí)驗(yàn)均采用Window 8操作系統(tǒng)，使用MATLAB 2018a平臺進(jìn)行仿真實(shí)驗(yàn),且PET和MRI數(shù)據(jù)均為同一AD患者同一時期的腦影像。首先在空域和變換域中采取相同的融合規(guī)則進(jìn)行比較，然后在ShearLab 3D域中將本文中的融合規(guī)則與現(xiàn)有的PET/MRI算法[17]進(jìn)行比較。選取信息熵(entropy,IE)、平均梯度(average gradient,AG)、空間頻率(spatial frequency,SF)、邊緣強(qiáng)度(edge intensity,EI)、綜合熵(comprehensive entropy,CE)作為客觀評價標(biāo)準(zhǔn),它們的值越大融合效果越好。

3.1 基于不同變換域的融合算法的比較

為了驗(yàn)證ShearLab 3D相較于空域算法和其它變換域算法的優(yōu)越性，在加權(quán)平均、像素取大、小波域和ShearLab 3D域內(nèi)融合PET/MRI圖像，其中小波域和ShearLab 3D域低頻系數(shù)使用加權(quán)平均法，高頻子帶系數(shù)使用體素絕對值取大法，且都為一層分解。圖2所示為PET和MRI圖像和不同算法產(chǎn)生的結(jié)果圖像。圖2a是AD患者注射18F-FDG的 PET圖像，其中左上是冠狀面，右上是矢狀面，左下是軸狀面;圖2b表示AD患者M(jìn)RI腦圖像;圖2c是空域加權(quán)平均結(jié)果;圖2d為空域像素取大結(jié)果;圖2e為變換域小波變換結(jié)果;圖2f為變換域ShearLab 3D的融合結(jié)果。

Fig.2 PET/MRI images in different domainsa—PET image b—MRI image c—weighted average d—take the bigger voxel e—wavelet f—ShearLab 3D

從主觀角度分析，空域的加權(quán)平均相對于其它方法融合圖像對比度略有降低導(dǎo)致邊緣輪廓模糊。像素取大法由于PET和MRI呈現(xiàn)原理不同導(dǎo)致融合結(jié)果丟失MRI圖像紋理信息，變換域的小波變換融合圖像的邊緣輪廓紋理等信息都受到了塊效應(yīng)的影響。ShearLab 3D融合圖像從很大程度上彌補(bǔ)了它們的不足，可以觀察到腦萎縮程度及相應(yīng)組織的代謝情況。

從客觀角度分析，觀察表1、圖2c和圖2e，其各項(xiàng)評價指標(biāo)數(shù)值普遍較低，尤其是圖2c中的AG,EI,SF指數(shù)過低，表示融合圖像紋理特征不夠清晰，灰度變化率低對比度差。圖2e中IE和CE較低融合圖像信息丟失嚴(yán)重。圖2d中的IE, EI,CE都較好，但融合圖像中只含有極少部分MRI圖像的細(xì)節(jié)紋理，缺乏實(shí)用性。圖2f的AG最好，其它指標(biāo)與最優(yōu)相比相差甚少，相對于其它算法，ShearLab 3D輪廓清晰、細(xì)節(jié)明顯更具有實(shí)用性。

Table 1 Objective evaluation of PET/MRI images in different domains

3.2 基于ShearLab 3D不同拉普拉斯算子的比較

在ShearLab 3D域?qū)⒈疚闹刑岢龅娜诤弦?guī)則與參考文獻(xiàn)[17]中提出的PET/MRI的融合規(guī)則進(jìn)行對比。ShearLab 3D為一層分解，圖3a是參考文獻(xiàn)[17]中融合規(guī)則的結(jié)果，圖3b是本文中融合算法結(jié)果。從視覺角度看,兩幅圖像均保留了原圖像的主要信息，但是圖3a中的腦溝回處邊緣輪廓不明顯且對比度較低，無法判斷腦萎縮狀況，降低了醫(yī)療診斷的實(shí)用性。本文中算法在腦溝回處邊緣輪清晰，并準(zhǔn)確地判斷海馬體葡萄糖代謝降低萎縮嚴(yán)重，右側(cè)額葉與右側(cè)顳葉葡萄糖代謝降低發(fā)生病變。

Fig.3 PET/MRI images with different fusion rules in ShearLab 3D domaina—reference [17] fusion results b—algorithm fusion results in this paper

從客觀指標(biāo)分析觀察表2，本文中算法融合結(jié)果的IE值略低于參考文獻(xiàn)[17]中的算法，但SF,AG,CE的評價指數(shù)都高于對比算法且EI有明顯提高，表明融合結(jié)果的體素活躍程度高，對微小細(xì)節(jié)的表示度高，圖像清晰。綜合主客觀評價,本文中算法的融合結(jié)果邊緣紋理特征豐富，清晰水平更好。

Table 2 Objective evaluation of PET/MRI images with different fusion rules in ShearLab 3D domain

4 結(jié) 論

針對醫(yī)學(xué)圖像融合要求融合圖像含有豐富清晰的紋理信息準(zhǔn)確的判斷病灶，本文中基于ShearLab 3D的方向敏感對圖像進(jìn)行分解，并根據(jù)方差識別強(qiáng)邊緣特征的特點(diǎn)將高頻子帶分為中高頻子帶。低中高頻子帶的融合規(guī)則由2維擴(kuò)展到3維，低頻子帶采用加權(quán)局部能量與拉普拉斯算子加權(quán)融合規(guī)則，并以銳化矩陣作為權(quán)重參量，改善了融合結(jié)果輪廓模糊的問題。中頻子帶以絕對值為活動度量增強(qiáng)融合圖像的邊緣信息。高頻子帶結(jié)合局部相位一致性、局部對比度、局部能量3個低層特征加權(quán)的融合規(guī)則，增強(qiáng)圖像的細(xì)節(jié)特征且增加了融合結(jié)果的亮度和對比度。通過實(shí)驗(yàn)結(jié)果證明，ShearLab 3D相對于空域和小波相比處理高維數(shù)據(jù)更有優(yōu)勢；相對融合規(guī)則，本文中提出的算法僅有信息熵略低，其它指標(biāo)都高于對比算法，其中平均梯度和邊緣強(qiáng)度表明,本文中融合結(jié)果含有更豐富且清晰的輪廓和紋理特征，空間頻率和綜合熵表明,融合圖像整體活躍度和融合質(zhì)量更高符合人類視覺效果。