基于多尺度邊緣保持分解與PCNN 的醫學圖像融合

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

（中北大學大數據學院山西省生物醫學成像與影像大數據重點實驗室，太原 030051）

0 概述

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

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

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

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

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

1 相關理論

1.1 加權最小二乘濾波與高斯濾波

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

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

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

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

1.2 非下采樣方向濾波器組

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

圖1 NDFB 基本結構Fig.1 Basic structure of NDFB

1.3 脈沖耦合神經網絡

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

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

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

2 本文醫學圖像融合方法

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

圖3 醫學圖像融合過程Fig.3 Process of medical image fusion

2.1 多尺度邊緣保持分解

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

1）對輸入圖像進行WLS 濾波分解得到細節層和基礎層。使用WLS 濾波對輸入圖像進行保邊平滑處理得到基礎層，細節層則由輸入圖像與基礎圖像做差得到。基礎層按照式（1）計算。設Cinput是輸入圖像，則細節層表示為：

其中，CD_layer和CB_layer是經過WLS 濾波分解得到的細節層和基礎層。由于WLS 濾波的保邊特性，因此基礎層含有強邊緣信息。

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

其中，CL_layer和CE_layer是基礎層經過高斯濾波分解得到的低頻層和邊緣層。基于WLS 濾波和高斯濾波的分解單元結構如圖4 所示。

圖4 基于WLS 濾波和高斯濾波的分解單元結構Fig.4 Decomposition unit structure based on WLS filtering and Gaussian filtering

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

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

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

2.2 圖像融合規則

2.2.1 低頻層融合規則

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

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

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

3）將上述計算的區域能量作為PCNN 模型的外部激勵，平均梯度構建的鏈接系數作為PCNN 模型的鏈接強度，通過PCNN 模型得到低頻子帶相應的點火圖

2.2.2 高頻層融合規則

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

3）將上述計算的改進的空間頻率作為PCNN 模型的外部激勵，將平均梯度構建的鏈接系數作為PCNN 模型的鏈接強度，通過PCNN 模型得到高頻方向子帶相應的點火圖

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

2.3 圖像重構

將得到的高頻層融合結果和低頻層融合結果重構為最終的融合結果：

3 實驗結果與分析

3.1 實驗參數設置

本文實驗數據選取自哈佛大學醫學院的公開數據，配準的醫學圖像是大小為256 像素×256 像素、顏色深度為8 bit 的灰度圖。

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

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

3.2 實驗結果分析

3.2.1 不同融合方法的比較

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

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

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

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

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

為更客觀地評價圖像的融合結果，選取邊緣評價因子QAB/F、平均梯度、互信息、標準差和空間頻率5 種指標進行客觀評價。3 組實驗的客觀評價指標如表1～表3 所示，其中加粗數據表示最優值。可以看出：在5 種評價指標中，各組實驗中本文方法的QAB/F均達到最大值，驗證了其在邊緣保持方面的可行性和有效性；此外，除第1 組實驗的空間頻率和第2 組、第3 組實驗的標準差之外，本文方法其余各實驗指標均達到最大值，表明其能較好地保留原圖像信息，體現圖像的紋理細節。

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

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

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

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

3.2.2 不同濾波器的比較

在本文方法中使用不同的保邊濾波器進行實驗，選取引導濾波器（Guided Filter，GF）和雙邊濾波器（Bilateral Filter，BF）與上文使用的WLS 濾波器進行對比。雙邊濾波器參數設置為δr=2，4，8，δs=0.1，0.2，0.4，引導濾波器參數設置為r=2，4，8，ε=0.12，0.22，0.42。表4～表6 列出了使用不同濾波器的客觀定量指標，其中，加粗數據表示最優值。可以看出，使用WLS 濾波器效果較使用其他濾波器更好。

表4 本文方法不同濾波器下的多發性腦梗塞MR-T1/MR-T2圖像融合評價指標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 圖像融合評價指標Table 5 Evaluation indexes of CT/MR image fusion oftoxoplasma gondii by the proposed method using different filters

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

4 結束語

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