基于低秩表示和NSST的醫學圖像融合

孟令玉 聶仁燦 何 敏 周冬明 李華光

(云南大學信息學院 云南 昆明 650500)

0 引 言

隨著成像技術的發展，醫學圖像在臨床診斷等方面發揮著越來越重要的作用，已經成為醫生疾病診斷、疾病治療的得力助手，極大降低了誤診風險，提高了對癥治療的效率、速率。與普通的光學成像方式不同，醫學圖像成像機理比較復雜，涉及到物理、化學和生物等多個學科，因此，單模態醫學圖像的功能比較單一，很難全面顯示人體器官組織的信息。例如，核磁共振成像(Magnetic Resonance Imaging,MRI)可以充分顯示軟組織的信息，但其骨骼等信息的探測存在很大缺陷。與MRI不同的是，計算機斷層掃描成像(Computer Tomography,CT)可以很好地捕獲骨骼信息，但其對軟組織等結構信息缺乏探測能力[1]。此外，其他模態的醫學圖像還有具有很強病灶探測能力的正電子發射計算機斷層成像(Positron Emission Tomography，PET)以及具有顯示血液流動信息能力的單光子發射計算機斷層成像(Single Photon Emission Tomography，SPECT)等。全面了解器官組織的各項信息可以更有針對性地進行臨床診斷、疾病治療，因此，有必要運用一定的圖像融合算法將不同模態的醫學圖像根據需求進行融合得到一幅顯示信息更為全面的醫學圖像。

醫學圖像融合作為圖像融合的重要應用領域之一，是當今國內外研究的熱點問題，在醫學領域有著廣泛的應用前景和實用價值[2]。近年來，針對醫學圖像的融合方法已經得到很大發展，這些算法大致可以分為兩類：基于頻率域的算法和基于空間域的算法[3]。空間域類算法主要是在原始圖像像素上直接進行相關處理，對圖像主要特征信息捕捉并融合。此類方法實現起來比較簡單，但是對圖像的特性分析不夠充分，人為痕跡較明顯，噪聲等干擾信息較多。相較而言，頻率域類方法一般通過多尺度分析工具對圖像在不同方向進行分解，可以充分挖掘圖像的邊緣、細節、輪廓等信息。目前，國內外研究學者已經提出很多可行的算法。王健等[4]提出一種非下采樣輪廓波變換(NonSubsampled Contourlet Transform，NSCT)域內基于增補小波變換和脈沖耦合神經網絡(Pulse Coupled Neural Network，PCNN)的圖像融合方法，在一定程度上彌補了傳統NSCT的缺陷，充分利用了其對圖像的多分辨率、多方向性等分析特性，融合效果細節信息豐富，輪廓清晰，符合人眼視覺習慣。劉雯敏等[5]提出一種基于PCNN和非線性濾波萬有引力的醫學圖像融合方法，很好地改善了融合圖像的視覺效果，有效地保留了源始圖像的各項信息，邊緣和紋理信息突出。吳銀芳等[6]提出一種基于模糊變換耦合最大熵的醫學圖像融合方法，可以有效克服融合過程中的邊緣模糊和偽輪廓，對提高融合圖像的邊緣信息強度、峰值信噪比等效果顯著。Ramla等[7]提出一種非下采樣剪切波變換(NonSubsampled Shearlet Transform，NSST)域內基于PCNN和形態學梯度的醫學圖像融合方法，其將圖像的形態學梯度作為PCNN的激發項，以此作為高頻信息部分的融合規則，很好地保留了源始醫學圖像的特征信息。同樣是在NSST域內，Hermessi等[8]提出一種將NSST結合全卷積神經網絡(Fully Convolutional Neural Network，FCNN)學習框架的圖像融合算法。通過訓練好的CNN對NSST分解后的高頻信息進行特征提取，很好地保留了源始醫學圖像各部分的信息，細節清晰，對比度、亮度等恰當，極大提高了醫學融合圖像的視覺效果，但是該算法執行起來耗時較長，對設備性能要求較高。

為使不同模態醫學圖像的特征信息充分融合到一幅圖像當中，提高圖像融合質量，增強融合結果的抗噪能力，豐富梯度、細節等信息，本文提出一種基于低秩表示(Latent Low Rank Representation，LatentLRR)[9]和NSST的MRI和CT醫學圖像融合方法。低秩表示方法LatentLRR可以將原始圖像分解成基礎圖、特征圖像和噪聲圖三部分，具有很強的去噪和特征提取能力。對特征圖使用方向引導濾波[10]進行融合；對基礎圖進行NSST分解后的高頻部分采用本文的局部梯度能量方法進行融合，低頻部分使用本文的加權改進拉普拉斯能量和SML(Sum of Modified Laplacian)進行融合。

1 基本理論

1.1 低秩表示LatentLRR

低秩表示LRR[11]是在已經確定學習字典之后，一種從稀疏噪聲中恢復低秩矩陣結構的有效模型。用最低秩數矩陣的線性組合表示所有的數據向量，近幾年一直是圖像領域的研究熱點。相比于稀疏表示，低秩表示可以更好地表示數據的整體結構。LRR的一般模型可以表示表示為：

(1)

s.t.XO=AZ

LatentLRR是在LRR基礎上發展而來，其不僅將已知的觀測矩陣作為學習字典，并考慮了隱藏的數據信息對學習字典的影響，具有很強的特征信息提取能力和去噪能力。為實現秩數最小化表示，可將LatentLRR看作帶有核范數的凸優化問題：

(2)

s.t.XO=XOZ+LXO+E

(a) 原圖 (b) 基礎 (c) 特征圖 (d) 噪聲圖圖1 LatentLRR分解結果

1.2 非下采樣剪切波變換NSST

小波變換只能對圖像在三個高通子帶進行高頻部分信息分析，方向性信息有限，而通過經典仿射系統理論將幾何分析與多分辨率分析結合發展而來的剪切波變換可以將圖像分解到各種方向上，可對圖像的細節、邊緣、輪廓等特征信息進行更加詳細的分析處理。NSST[13-14]不僅具有平移不變性等優良特性，并且由于分解過程沒有使用采樣運算，NSST不會產生吉布斯現象。在二維圖像空間中，帶有合成膨脹的仿射系統可以表示為：

MXY(ψ)={ψl,j,k(x)=|detXj/2|ψ(YlXjx-k)}

(3)

(4)

(5)

為實現函數ψi,j,k(x)需要以下頻率范圍作為支撐集:

(6)

(7)

為實現NSST，先對源圖像進行非下采樣金字塔(Nosubsampled Pyramid,NSP)得到高頻信息部分和低頻信息部分。接著用以上方法對高頻信息部分進行局部方向性分析。然后再對低頻信息部分繼續使用NSP分解，依次迭代操作完成對源始圖像的分解過程。相應的對各子帶信息處理后，通過逆NSST運算便可以得到重構圖像。

2 融合規則及融合方法

圖像融合一般包括三個研究方向：圖像分解方法、圖像融合規則和融合結果評價方法。當圖像經過恰當的多尺度分解方法分解為高頻部分和低頻部分之后，需要一定的融合規則將對應的信息進行整合。融合規則作為圖像融合過程中的重要組成部分，對源始圖像信息的正確選取和最終融合結果的優良具有決定性的作用。

2.1 特征圖融合規則

為充分將源始圖像經過LatLRR分解得到的特征信息進行有效融合處理，更加高效地保留相關重要特征信息，本文將引導濾波[15]與一組方向矩陣結合作為特征圖的融合規則。引導濾波是在局部線性模型的基礎上發展而來，在濾波過程中，算法實施起來簡單，使用一幅引導圖指導整個濾波過程，效率高，對圖像的邊緣等信息有著很好的平滑特性，是一種邊緣保持效果最好的邊緣保持濾波器[16-17]。本文中用GFF表示引導濾波器算法。

設輸入圖像為Pn，則通過以下方式使用平均濾波器AVA將源始圖像進行分解：

Bn=Pn*AVA

Dn=Pn-B

(8)

式中：*為卷積運算符；Bn、Dn為分解后得到的圖像。為得到Bn、Dn的顯著性成分圖，并為提取更加豐富的細節信息，在計算權重圖時使用一組方向性矩陣。設拉普拉斯濾波器為LAPF，wi(i=1,2,…,8)為方向性矩陣，則有:

(9)

式中：Sn為顯著性成分圖。8個方向性矩陣分別為：

(10)

(11)

通過引導濾波分別得到Bn、Dn的權重圖：

(12)

式中：r1、ε1和r2、ε2為引導濾波相關參數。最后，可以通過以下方法得到特征圖的融合結果Ff：

(13)

2.2 NSST低頻部分融合規則

為提高融合圖像的對比度、進一步提升對原始圖像細節等信息的豐富程度，本文采用窗口大小為3×3的加權改進拉普拉斯能量和作為低頻融合規則。設經過NSST分后的低頻圖為DBn(n=1,2,…,N)，則本文中在(x,y)位置處的加權拉普拉斯能量和MLn(x,y)由以下方式計算得到：

MLn(x,y)=

|1.6DBn(x,y)-0.8DBn(x-1,y)-0.8DBn(x+1,y)|+

|1.6DBn(x,y)-0.8DBn(x,y-1)-0.8DBn(x,y+1)|+

|2DBn(x,y)-DBn(x-1,y-1)-DBn(x+1,y+1)|+

|2DBn(x,y)-DBn(x+1,y-1)-DBn(x-1,y+1)|

(14)

然后，通過比較各個圖像在(x,y)處的能量和得到融合圖像Ffdb：

Ffdb(x,y)=DBn(x,y)

ifMLn(x,y)=max(ML1,ML2,…,MLN)

(15)

2.3 NSST高頻部分融合規則

經過多層次、多方向性分解后，得益于NSST的優良特性，圖像的特征信息已經被充分表示出來。因此，本文中對NSST分解后的高頻信息部分采用局部梯度能量的方法進行融合。

首先，分別計算高頻信息部分DGn在橫向和縱向的梯度DGngx、DGngy：

DGngx=gradientx(DGn)

DGngy=gradienty(DGn)

(16)

式中：gradientx、gradienty分別為橫向和縱向的梯度算子。

然后，取得綜合性梯度信息DGng：

(17)

通過以下方式計算DGng在(x,y)處的能量：

(18)

式中：w=3為低頻圖像;We為大小為3×3的局部窗，

(19)

設能量最大值為Fn(x,y):

Fn(x,y)=max{En(x+i,y+j)|1≤i,j≤3}

(20)

通過比較不同高頻信息的局部梯度能量得到權重:

(21)

從而，高頻部分融合后的圖像系數可以表示為：

(22)

通過對Ffdb和DGf使用逆NSST得到基礎圖融合結果。

2.4 融合方法

首先使用LatLRR將源始圖像分解為特征圖和基礎圖，對相應的特征圖使用引導濾波進行融合得到特征融合圖。然后將基礎圖進行NSST分解得到高頻圖、低頻圖，對低頻圖使用加權改進拉普拉斯能量和進行融合，對高頻圖運用局部梯度能量進行融合，對低頻融合結果和高頻融合結果使用逆NSST進行圖像重構得到基礎圖的融合結果。由于特征圖和基礎圖為互補的關系，因此最后通過加法運算將基礎融合圖和特征融合圖結合得到最終融合結果。融合框架如圖2所示。

圖2 本文算法融合框架

3 實驗結果及分析

為驗證本文算法的有效性，選取三組已經配準好的MRI-CT圖作為實驗對象，并與幾種經典算法DWT(Discrete Wavelet Transform)[18]、FPDE(Fourth Order Partial Differential Equations)[19]、MSVD(Multi-resolution Singular Value Decomposition)[20]、NSCT-SF-PCNN[21]、CUR(Curvelet Transform)[22]作比較。其中，DWT、CUR的高頻部分、低頻部分分別采用取最大值、取平均值的融合方法。實驗結果如圖3-圖5所示。

(a) 源MRI1 (b) 源CT2 (c) 文獻[18] (d) 文獻[19]

(a) 源MRI2 (b) 源CT2 (c) 文獻[18] (d) 文獻[19]

(e) 文獻[20] (f) 文獻[21] (g) 文獻[22] (h) 本文算法圖4 第二組融合結果

(a) 源MRI3 (b) 源CT3 (c) 文獻[18] (d) 文獻[19]

(e) 文獻[20] (f) 文獻[21] (g) 文獻[22] (h) 本文算法圖5 第三組融合結果

為方便客觀評價各算法的融合結果，選取PSNR、VIF、STD、Qabf作為評價指標。PSNR可以反映算法在融合過程中對噪聲的抑制能力；VIF是判斷融合圖像結果對人眼視覺習慣符合程度的重要指標；STD可以對圖像對比度信息進行評測；Qabf是反映融合結果中對源始圖像邊緣信息保留情況的有效指標。以上各項評測方法均為值越大，融合效果越好。圖3-圖5的各項評測指標分別如表1-表3所示。

表1 圖3各算法的評測指標

表2 圖4各算法的評測指標

表3 圖5各算法的評測指標

從圖3-圖5各算法的融合結果可以看出，文獻[19-20]邊緣、細節信息不夠清晰，有很強的涂抹感。文獻[18，21-22]對CT圖的信息保留較少，并帶有噪點。文獻[18]CT圖信息明顯融合不均，文獻[22]細節信息豐富，但亮度、對比度不及本文算法。表1-表3更加客觀地反映了各算法融合結果的差異性，這與圖3-圖5的視覺感受相符。顯然，本文算法在各項指標中占有優勢，具有很強的抗噪聲能力，原始圖像的細節信息保留最好，對比度、梯度信息豐富，更加符合人眼視覺習慣，更能滿足臨床診斷的需求。

4 結 語

醫學圖像融合是提高臨床診斷效率的重要方法。本文運用方向引導濾波、加權改進拉普拉斯能量和等作為融合規則，提出一種基于低秩表示和NSST的MRI/CT醫學圖像融合方法，其實現簡單，在一定程度上提高了融合圖像的抗噪能力、細節梯度信息，并擁有不錯的視覺效果。今后將努力完善本文算法，挖掘本文算法在多通道醫學圖像中的應用價值。