基于多尺度邊緣保持分解與PCNN 的醫(yī)學(xué)圖像融合

郭淑娟，高 媛，秦品樂，王麗芳

（中北大學(xué)大數(shù)據(jù)學(xué)院山西省生物醫(yī)學(xué)成像與影像大數(shù)據(jù)重點實驗室，太原 030051）

0 概述

圖像融合是指將不同模態(tài)的相關(guān)圖像合并為一幅新圖像。在此過程中，要求生成的融合圖像具有良好的清晰度，并使不同模態(tài)的圖像特征盡可能融合到一幅圖像中。圖像融合在多模態(tài)醫(yī)學(xué)圖像中的應(yīng)用是一個熱門的研究方向。隨著現(xiàn)代醫(yī)學(xué)成像技術(shù)的發(fā)展，結(jié)合不同的醫(yī)學(xué)成像方式能夠獲得器官和組織的全面信息，例如計算機(jī)斷層掃描（Computed Tomography，CT）圖像僅顯示高密度組織，磁共振（Magnetic Resonance，MR）圖像僅顯示高水分組織，正電子發(fā)射斷層掃描（Positron Emission Tomography，PET）圖像僅反映不同人體之間組織強(qiáng)度的差異［1］，而將這些醫(yī)學(xué)圖像的互補(bǔ)特征融合到一幅圖像中，能夠提高醫(yī)生的診斷效率和準(zhǔn)確性。

圖像融合技術(shù)可分為空間域方法和變換域方法兩類，其中，變換域方法可以更好地表征圖像特征，彌補(bǔ)空間域方法在細(xì)節(jié)提取方面的不足。在變換域方法中，多尺度分析方法通過將源圖像映射至不同的尺度空間提取潛在的重要信息，符合人眼視覺系統(tǒng)的生理機(jī)制，被認(rèn)為是一種主流的圖像融合方法［2］?，F(xiàn)有的多尺度分析方法主要有基于離散小波變換［3］、基于剪切波變換［4］、基于輪廓波變換［5］、基于非下采樣輪廓波變換（Non-Subsampled Contourlet Transform，NSCT）［6］和基于非下采樣剪切波變換（Non-Subsampled Shearlet Transform，NSST）［7］的算法，其中，基于NSCT 和NSST 的算法具有平移不變性，能夠很好地抑制偽吉布斯現(xiàn)象。

脈沖耦合神經(jīng)網(wǎng)絡(luò)（Pulse Coupled Neural Network，PCNN）是一種通過研究哺乳動物視覺皮層而獲得的具有仿生機(jī)制的單層神經(jīng)網(wǎng)絡(luò)模型［8］，其不需要學(xué)習(xí)或者訓(xùn)練，且具有全局耦合性和脈沖發(fā)射特性，能夠提取圖像的全局特征并且增強(qiáng)圖像細(xì)節(jié)，有利于圖像的實時處理。為提高融合圖像的視覺感知效果并充分利用圖像的空間信息，一些多尺度分析與PCNN 相結(jié)合的圖像融合方法相繼被提出［9］，如NSCT-SF-PCNN［10］和NSST-PAPCNN［11］。然而，傳統(tǒng)的多尺度分析方法大多采用線性濾波器，由于無法保留邊緣從而導(dǎo)致分解階段的強(qiáng)邊緣處出現(xiàn)模糊，使得邊緣處產(chǎn)生光暈［12］。

為解決該問題，研究人員將保邊濾波器引入圖像融合領(lǐng)域。保邊濾波器具有平移不變性、空間一致性和邊緣保護(hù)性能，而且計算效率較高。在此方面：文獻(xiàn)［13］提出一種多尺度雙邊濾波器對源圖像進(jìn)行分解；文獻(xiàn)［14］提出一種基于引導(dǎo)濾波的圖像融合方案，其對源圖像進(jìn)行兩尺度分解，通過引導(dǎo)濾波優(yōu)化權(quán)重實現(xiàn)了空間一致性；文獻(xiàn)［15］提出以多尺度極值濾波對源圖像進(jìn)行多尺度分解，采用對比度融合規(guī)則使融合圖像更符合人眼視覺感知特性；文獻(xiàn)［16］提出一種加權(quán)最小二乘濾波與引導(dǎo)濾波相結(jié)合的新方案，與只基于引導(dǎo)濾波的方法相比，該方案可以獲得更好的效果。盡管上述方法在一定程度上保護(hù)了邊緣信息，但是基礎(chǔ)層主要包含粗尺度結(jié)構(gòu)信息，直接對其融合可能在之后的融合過程中丟失或模糊保邊濾波保留的顯著的邊緣特征［17］。

為保留圖像完整清晰的邊緣信息，提高融合圖像的視覺感知效果，本文提出一種新的醫(yī)學(xué)圖像融合方法。將多尺度邊緣保持分解方法引入圖像分解過程，對已配準(zhǔn)的源圖像進(jìn)行分解得到低頻層和高頻層。在此基礎(chǔ)上，利用改進(jìn)的空間頻率以及區(qū)域能量激勵PCNN 對兩者進(jìn)行融合，并通過逆變換獲得最終的融合圖像。

1 相關(guān)理論

1.1 加權(quán)最小二乘濾波與高斯濾波

文獻(xiàn)［18］提出一種基于加權(quán)最小二乘（Weighted Least Square，WLS）優(yōu)化框架的非線性濾波器。給定一個輸入圖像g，WLS 濾波器的目的是使濾波后的平滑圖像u盡可能與輸入圖像g近似，則濾波后的圖像u可以表示為以下目標(biāo)函數(shù)：

其中，下標(biāo)p表示像素的空間位置，（up-gp）2用以保證濾波后的平滑圖像u更接近輸入圖像g，是一個正則化項，通過求u的偏導(dǎo)數(shù)來實現(xiàn)平滑，λ是正則化參數(shù)，用以調(diào)整兩項之間的平衡，ax，p和ay，p代表平滑因子。ax，p和ay，p的計算公式如下：

其中，l是輸入圖像g的對數(shù)亮度通道，指數(shù)α確定對g的梯度的敏感度，而ε是一個非常小的常數(shù)，通常為0.000 1，可防止在g恒定的區(qū)域被零除。

高斯濾波是一種線性濾波，由于沒有考慮相鄰像素的影響，其對圖像濾波時去除高頻信息的同時也模糊了邊緣，因此在圖像的平滑和去噪中被廣泛應(yīng)用。高斯濾波公式可表示為：

1.2 非下采樣方向濾波器組

文獻(xiàn)［19］提出一種方向濾波器組（Directional Filter Bank，DFB）。DFB 將圖像劃分為具有方向性的子帶，捕捉圖像幾何特征，但是其中的下采樣操作使其缺乏平移不變性。非下采樣方向濾波器組（Non-subsampled Directional Filter Bank，NDFB）［20］與DBF 類似，可以對圖像進(jìn)行多方向分解，提取圖像方向細(xì)節(jié)特征，但是其進(jìn)行上采樣操作，因而能夠獲得位移不變的方向擴(kuò)展。NDFB 是由扇形濾波器和棋盤濾波器構(gòu)成的四方向濾波器組，基本結(jié)構(gòu)如圖1 所示。

圖1 NDFB 基本結(jié)構(gòu)Fig.1 Basic structure of NDFB

1.3 脈沖耦合神經(jīng)網(wǎng)絡(luò)

圖2 PCNN 的簡化模型Fig.2 Simplified model of PCNN

脈沖耦合神經(jīng)網(wǎng)絡(luò)（PCNN）是一個M×N的二維網(wǎng)絡(luò)，其中每一個神經(jīng)元對應(yīng)于圖像的一個特定的像素。每個神經(jīng)元由三部分組成，分別是接收域、調(diào)制域和脈沖產(chǎn)生器。一種PCNN 的簡化模型和數(shù)學(xué)表達(dá)式如圖2 和式（5）所示。其中：下標(biāo)L表示圖像解等級；Fij(n)和Lij(n)分別是PCNN 的饋送輸入和鏈接輸入，處于PCNN 模型的接收域，F(xiàn)ij(n)通常接收的是圖像像素的歸一化灰度值，Lij(n)則通過突觸權(quán)重與8 個相鄰神經(jīng)元的先前放電狀態(tài)相關(guān)；aL和VL是Lij(n)的時間衰減常數(shù)和振幅增益；Uij(n)代表內(nèi)部活動，處于PCNN 模型的調(diào)制域，通過Fij(n)和Lij(n)進(jìn)行非線性調(diào)制得到，其中參數(shù)β是連接強(qiáng)度。

脈沖發(fā)生域控制脈沖的發(fā)生：將Uij與動態(tài)閾值θij進(jìn)行比較，若Uij＞θij，Yij=1 即神經(jīng)元點火，產(chǎn)生脈沖輸出；反之Yij=0，aθ和Vθ是θij的時間衰減常數(shù)和振幅增益。將上述過程多次迭代直到滿足設(shè)定的迭代條件，把此時由神經(jīng)元點火形成的點火映射圖作為輸出結(jié)果。

2 本文醫(yī)學(xué)圖像融合方法

對已配準(zhǔn)的源圖像進(jìn)行圖像融合的主要步驟包括圖像分解、圖像融合和圖像重構(gòu)。本文在圖像分解過程中引入多尺度邊緣保持分解方法，對已配準(zhǔn)的源圖像進(jìn)行分解得到低頻層和高頻層，而在圖像融合過程中針對得到的低頻層和高頻層不同的特點，分別采用不同的融合規(guī)則進(jìn)行融合，在圖像重構(gòu)中獲得最終的融合的圖像。融合過程如圖3所示。

圖3 醫(yī)學(xué)圖像融合過程Fig.3 Process of medical image fusion

2.1 多尺度邊緣保持分解

WLS 濾波器可以有效地在模糊和銳化之間進(jìn)行最佳折衷，更好地保留邊緣，解決雙邊濾波在邊緣處容易產(chǎn)生梯度反轉(zhuǎn)的問題，彌補(bǔ)雙邊濾波無法很好地在多尺度上提取到細(xì)節(jié)信息的不足，從而更準(zhǔn)確地提取多尺度信息。由于高斯濾波器只考慮像素的空間分布而沒有考慮鄰近像素的影響，因此具有對圖像進(jìn)行濾波時能夠剔除細(xì)節(jié)信息、模糊邊緣的特點。結(jié)合兩種濾波的特性，本文對源圖像進(jìn)行多尺度分解。分解過程由以下步驟組成：

1）對輸入圖像進(jìn)行WLS 濾波分解得到細(xì)節(jié)層和基礎(chǔ)層。使用WLS 濾波對輸入圖像進(jìn)行保邊平滑處理得到基礎(chǔ)層，細(xì)節(jié)層則由輸入圖像與基礎(chǔ)圖像做差得到?；A(chǔ)層按照式（1）計算。設(shè)Cinput是輸入圖像，則細(xì)節(jié)層表示為：

其中，CD_layer和CB_layer是經(jīng)過WLS 濾波分解得到的細(xì)節(jié)層和基礎(chǔ)層。由于WLS 濾波的保邊特性，因此基礎(chǔ)層含有強(qiáng)邊緣信息。

2）為避免光譜損失，使用高斯濾波器再次分解基礎(chǔ)層提取邊緣信息，得到低頻層和邊緣層。低頻層圖像依據(jù)式（4）計算，則邊緣層計算公式為：

其中，CL_layer和CE_layer是基礎(chǔ)層經(jīng)過高斯濾波分解得到的低頻層和邊緣層?；赪LS 濾波和高斯濾波的分解單元結(jié)構(gòu)如圖4 所示。

圖4 基于WLS 濾波和高斯濾波的分解單元結(jié)構(gòu)Fig.4 Decomposition unit structure based on WLS filtering and Gaussian filtering

3）將邊緣層CE_layer和細(xì)節(jié)層CD_layer進(jìn)行疊加，構(gòu)建圖像的高頻層CH_layer。同時為提取高頻層系數(shù)中存在的大量幾何特征，對高頻層使用NSDFB 進(jìn)行方向分析，捕獲高頻層系數(shù)不同方向的特征。邊緣保持分解框架如圖5 所示。

圖5 邊緣保持分解框架Fig.5 Frame of edge-preserving decomposition

將得到的低頻層CL_layer作為下一級分解的輸入圖像，經(jīng)過多級分解實現(xiàn)多尺度分解。通過上述分解模式最終將源圖像分解為一個低頻層和一系列不同方向的高頻層方向子帶。將源圖像A 和源圖像B通過多尺度分解得到的低頻層表示為和，不同方向的高頻層方向子帶表示為和，分別表示源圖像A 和源圖像B 在第i級別第d方向上的高頻層。

2.2 圖像融合規(guī)則

2.2.1 低頻層融合規(guī)則

低頻層圖像是源圖像去除高頻信息之后的近似圖像，包含圖像大量的能量信息，其中子帶相鄰系數(shù)之間存在區(qū)域相關(guān)性［21］。區(qū)域能量將鄰域內(nèi)每個中心像素的能量表征為其自身和鄰域內(nèi)像素值的平方和，考慮了圖像局部特征，減少了融合圖像中的不連續(xù)性，邊緣信息在該融合模式下能夠很好地得以保留。因此，本文將區(qū)域能量作為PCNN的外部激勵融合低頻子帶。首先計算低頻層的區(qū)域能量作為PCNN 的外部激勵，然后使用表征圖像清晰度的平均梯度作為PCNN 的連接強(qiáng)度，最后通過PCNN 輸出的點火圖進(jìn)行融合。低頻層融合過程如下：

其中，W是一個3×3 的滑動窗口。為體現(xiàn)低頻層子帶的平滑特性，選取W為：

其中，η是一個大于0 的常數(shù)，用于調(diào)節(jié)鏈接強(qiáng)度的值，本文設(shè)置為0.2(i，j)是計算所得的平均梯度。的計算公式如下：

3）將上述計算的區(qū)域能量作為PCNN 模型的外部激勵，平均梯度構(gòu)建的鏈接系數(shù)作為PCNN 模型的鏈接強(qiáng)度，通過PCNN 模型得到低頻子帶相應(yīng)的點火圖

2.2.2 高頻層融合規(guī)則

高頻層反映了圖像的紋理細(xì)節(jié)和邊緣信息，對最終的融合結(jié)果有著至關(guān)重要的影響。改進(jìn)的空間頻率能夠反映圖像的邊緣和細(xì)節(jié)信息，同時表示圖像灰度值的變化，因此，本文采用改進(jìn)的空間頻率激PCNN 對高頻層進(jìn)行融合。首先計算每一級不同方向的高頻層方向子帶改進(jìn)的空間頻率作為PCNN 的外部刺激，然后使用表征圖像清晰度的平均梯度作為PCNN 的連接強(qiáng)度，最后通過PCNN 輸出的點火圖進(jìn)行融合。高頻層融合過程如下：

3）將上述計算的改進(jìn)的空間頻率作為PCNN 模型的外部激勵，將平均梯度構(gòu)建的鏈接系數(shù)作為PCNN 模型的鏈接強(qiáng)度，通過PCNN 模型得到高頻方向子帶相應(yīng)的點火圖

通過NDFB 逆變換，可以得到第i層的高頻層融合結(jié)果

2.3 圖像重構(gòu)

將得到的高頻層融合結(jié)果和低頻層融合結(jié)果重構(gòu)為最終的融合結(jié)果：

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

3.1 實驗參數(shù)設(shè)置

本文實驗數(shù)據(jù)選取自哈佛大學(xué)醫(yī)學(xué)院的公開數(shù)據(jù)，配準(zhǔn)的醫(yī)學(xué)圖像是大小為256 像素×256 像素、顏色深度為8 bit 的灰度圖。

在圖像分解階段，基于WLS 分解的參數(shù)主要包括平滑參數(shù)λ和邊緣保留參數(shù)α。λ是一個正則化參數(shù)，用于保持?jǐn)?shù)據(jù)項和平滑度項之間的平衡，增加λ可產(chǎn)生更平滑的圖像。α通過非線性縮放梯度來確定對梯度的敏感程度，增加α可保留更清晰的邊緣。本文參考文獻(xiàn)［16］，設(shè)置λi=[0.075，0.6，4.8]，α=0.6，其中，i=1，2，…，L?；诟咚篂V波分解的主要參數(shù)是用于控制圖像平滑效果的高斯分布的標(biāo)準(zhǔn)差σ，第1 級設(shè)置為σ1=1.1，其他級設(shè)置為σi+1=kσi，其中，k是相鄰級之間的分解比例因子，設(shè)置為k=2。圖像分解等級設(shè)置為L=1，2，3，方向分析的分解水平分別設(shè)置為第1 層方向數(shù)為16 方向，第2 層和第3 層分解分別為8 方向和4 方向，采用的方向濾波器為“vk”，將第L級的得到的低頻圖像作為下一級分解的輸入圖像以實現(xiàn)多尺度分解。

在圖像融合階段，本文參考文獻(xiàn)［9］，將用于融合圖像的PCNN 的參數(shù)設(shè)置為：aL=10，aθ=0.815，VL=1.0，Vθ=10，n=200

3.2 實驗結(jié)果分析

3.2.1 不同融合方法的比較

本文采用多發(fā)性腦梗塞MR-T1 和MR-T2 圖像、腦弓形蟲CT 和MR 圖像、正常腦部腦CT 和MR 圖像進(jìn)行3 組仿真實驗，在64 位Windows10 操作系統(tǒng)、Matlab2014b 平臺上進(jìn)行仿真。選取以下5 種融合方法作為對比方法：1）基于NSCT 的圖像融合方法（NSCT）；2）基于多尺度變換和稀疏表示的圖像融合方法（MST-SR）［22］；3）基于NSCT域下空間頻率激勵PCNN的圖像融合方法（NSCT-SF-PCNN）［10］；4）基于NSST 與自適應(yīng)PCNN 的圖像融合方法（NSST-PAPCNN）［11］；5）基于引導(dǎo)濾波的圖像融合方法（GFF）［14］。NSCT采用高頻子帶區(qū)域能量取大、低頻子帶加權(quán)平均的融合規(guī)則，MST-SR 的多尺度分解方法選取為LP，分解等級為4。6 種方法的實驗結(jié)果如圖6～圖8 所示。

圖6 不同方法的多發(fā)性腦梗塞MR-T1/MR-T2圖像融合結(jié)果Fig.6 Fusion results for MR-T1/MR-T2 image of multiple cerebral infarction by different methods

圖7 不同方法的腦弓形蟲CT/MR 圖像融合結(jié)果Fig.7 Fusion results for CT/MR image of toxoplasma gondii by different methods

圖8 不同方法的正常腦部CT/MR 圖像融合結(jié)果Fig.8 Fusion results for CT/MR image of normal brain by different methods

由實驗結(jié)果可以看出：NSCT 和NSCT-SF-PCNN的融合結(jié)果丟失了骨質(zhì)信息，同時NSCT 的融合圖像整體亮度偏低；NSCT-SF-PCNN 在第3 組實驗中的融合結(jié)果MR 信息嚴(yán)重缺失并且出現(xiàn)邊緣模糊的現(xiàn)象；NSST-PAPCNN 的融合結(jié)果雖然較好地保存了骨質(zhì)信息，但是在邊緣處產(chǎn)生了邊緣模糊和光暈現(xiàn)象；MST-SR 能夠較好地保存骨質(zhì)信息和邊緣信息，但是在第2 組實驗中MST-SR 融合圖像MR信息有所缺失，存在對比度下降和細(xì)節(jié)丟失的問題；GFF 的方法能夠有效保護(hù)邊緣信息，但是第1 組和第2 組實驗結(jié)果存在明顯的細(xì)節(jié)丟失現(xiàn)象。

為更客觀地評價圖像的融合結(jié)果，選取邊緣評價因子QAB/F、平均梯度、互信息、標(biāo)準(zhǔn)差和空間頻率5 種指標(biāo)進(jìn)行客觀評價。3 組實驗的客觀評價指標(biāo)如表1～表3 所示，其中加粗?jǐn)?shù)據(jù)表示最優(yōu)值?？梢钥闯觯涸? 種評價指標(biāo)中，各組實驗中本文方法的QAB/F均達(dá)到最大值，驗證了其在邊緣保持方面的可行性和有效性；此外，除第1 組實驗的空間頻率和第2 組、第3 組實驗的標(biāo)準(zhǔn)差之外，本文方法其余各實驗指標(biāo)均達(dá)到最大值，表明其能較好地保留原圖像信息，體現(xiàn)圖像的紋理細(xì)節(jié)。

表1 不同方法的多發(fā)性腦梗塞MR-T1/MR-T2 圖像融合評價指標(biāo)Table 1 Evaluation indexes of MR-T1/MR-T2 image fusion of multiple cerebral infarction by different methods

表2 不同方法的腦弓形蟲CT/MRI 圖像融合評價指標(biāo)Table 2 Evaluation indexes of CT/MRI image fusion of toxoplasma gondii by different methods

表3 不同方法的正常腦部CT/MRI 圖像融合評價指標(biāo)Table 3 Evaluation indexes of CT/MRI image fusion of normal brain by different methods

通過對以上主觀視覺和客觀指標(biāo)的分析可知，本文方法能夠很好地保留邊緣信息，所得圖像具有較高的對比度和豐富的細(xì)節(jié)信息，較對比方法能夠得到更好的視覺效果。

3.2.2 不同濾波器的比較

在本文方法中使用不同的保邊濾波器進(jìn)行實驗，選取引導(dǎo)濾波器（Guided Filter，GF）和雙邊濾波器（Bilateral Filter，BF）與上文使用的WLS 濾波器進(jìn)行對比。雙邊濾波器參數(shù)設(shè)置為δr=2，4，8，δs=0.1，0.2，0.4，引導(dǎo)濾波器參數(shù)設(shè)置為r=2，4，8，ε=0.12，0.22，0.42。表4～表6 列出了使用不同濾波器的客觀定量指標(biāo)，其中，加粗?jǐn)?shù)據(jù)表示最優(yōu)值?？梢钥闯?，使用WLS 濾波器效果較使用其他濾波器更好。

表4 本文方法不同濾波器下的多發(fā)性腦梗塞MR-T1/MR-T2圖像融合評價指標(biāo)Table 4 Evaluation indexes of MR-T1/MR-T2 image fusion of multiple cerebral infarction by the proposed method using different filters

表5 本文方法不同濾波器下的腦弓形蟲CT/MR 圖像融合評價指標(biāo)Table 5 Evaluation indexes of CT/MR image fusion oftoxoplasma gondii by the proposed method using different filters

表6 本文方法不同濾波器下的正常腦部CT/MR 圖像融合評價指標(biāo)Table 6 Evaluation indexes of CT/MR image fusion of normal brain by the proposed method using different filters

4 結(jié)束語

傳統(tǒng)基于多尺度分析的圖像融合方法因無法保留邊緣特征而導(dǎo)致邊緣處產(chǎn)生光暈。針對該問題，本文提出一種多尺度邊緣保持分解與PCNN 結(jié)合的醫(yī)學(xué)圖像融合方法。在分解過程保護(hù)邊緣的同時進(jìn)行多尺度多方向分解，并采用PCNN 作為融合規(guī)則，從而改善融合圖像的視覺感知效果。實驗結(jié)果表明，該方法能夠突出醫(yī)學(xué)圖像的邊緣輪廓并增強(qiáng)圖像細(xì)節(jié)，有利于將更多的顯著特征從源圖像分離并轉(zhuǎn)移到融合圖像中。本文方法使用了大量參數(shù)，而參數(shù)設(shè)置將影響融合結(jié)果的質(zhì)量，因此，下一步將從參數(shù)自適應(yīng)角度著手對其進(jìn)行優(yōu)化。