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

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

引 言

中國已進入老齡化階段，預計老齡化進程在未來十幾年內(nèi)將不斷加速[1]。阿爾茨海默病(Alzheimer disease,AD)在老年群體中發(fā)病率不斷增加，并給社會、家庭帶來不可估計的經(jīng)濟壓力[2]。影像學是AD診斷的重要方式,包括計算機斷層掃描技術(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腦影像清晰地看到腦溝回加深加寬、腦室擴大程度、海馬萎縮等，但由于很多疾病可以導致腦萎縮，僅僅依靠MRI幾乎不能確診。PET通過在人體中注射放射性物質(zhì)，可以無創(chuàng)生理地研究人體內(nèi)生物分子狀態(tài)。AD患者表現(xiàn)為神經(jīng)突觸的丟失或機能的降低，PET示蹤劑18氟-氟脫氧葡萄糖(18F-fluorodeoxyglucose，18F-FDG)通過葡萄糖代謝高低反映神經(jīng)突觸的活性，表現(xiàn)在皮層區(qū)域葡萄糖代謝下降，但圖像空間分辨率低。因此,采用圖像融合技術(shù)生成PET/MRI圖像，利用MRI清晰的結(jié)構(gòu)精準判斷PET葡萄糖代謝低的位置，提高AD的特異性和準確性。

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

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

1 ShearLab 3D理論

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

(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)過頭動校正、配準、歸一化處理得到兩個待處理對象圖像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),其中高頻子帶以方差為閾值獲得強邊緣子帶為中頻子帶Cm,A(x,y,z)和留下的高頻子帶Ch,1,A(x,y,z),Ch,1,B(x,y,z)。本文中提出將低中高頻融合規(guī)則擴展為3維，以加權(quán)局部能量與拉普拉斯算子加權(quán)和融合低頻系數(shù)并以銳化矩陣作為權(quán)重參量捕捉圖像中灰度跳變區(qū)域，以絕對值為活動度量融合中頻子帶，以3個3維擴展的底層特征加權(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圖像的腦溝回等細節(jié)特征不足。

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

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

(10)

(11)

W3×3×3=

(12)

設LWSTML,s(x,y,z)是低頻子帶位置(x,y,z)處二十六鄰域的改進拉普拉斯算子的加權(quán)和(weight sum of twenty six-neighborhood based modified laplacian，WSTML)。LTML,s[15]為改進的拉普拉斯算子(twenty six-neighborhood based modified laplacian，TML)，將相鄰點的權(quán)重設為歐氏距離。為了充分利用相鄰體素間的關系,將LTML,s擴展到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)，為了更好地融合圖像的邊緣信息在高頻子帶中提取包含源圖像強邊緣信息中頻子帶。由于方差可以判斷圖像的紋理、邊緣、平滑區(qū)域，其中紋理區(qū)方差最大，邊緣次之，平滑最小[16]。因此以方差為閾值得到中頻子帶，且方差的選取與圖像特點有關。

(16)

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

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

G(x,y,z)=

(17)

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

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

(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(x,y,z)=1時,PET圖像的中頻系數(shù)為融合圖像的中頻子帶系數(shù)。

(20)

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

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

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

(21)

式中,θk是k的方向角，An,θk為第n個傅里葉分量的幅值和角度，ε=0.001是正數(shù)去除圖像中的DC分量，Eθk(x,y,z)可見參考文獻[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)的體素點。

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(x,y,z)=1時,PET圖像系數(shù)為融合圖像的高頻系數(shù)。

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

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

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

為了驗證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患者MRI腦圖像;圖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)平均相對于其它方法融合圖像對比度略有降低導致邊緣輪廓模糊。像素取大法由于PET和MRI呈現(xiàn)原理不同導致融合結(jié)果丟失MRI圖像紋理信息，變換域的小波變換融合圖像的邊緣輪廓紋理等信息都受到了塊效應的影響。ShearLab 3D融合圖像從很大程度上彌補了它們的不足，可以觀察到腦萎縮程度及相應組織的代謝情況。

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

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

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

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

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

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

4 結(jié) 論

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