基于導(dǎo)向濾波器的醫(yī)學(xué)圖像融合方法

郭 盼，何文超，梁龍凱，張 萌，呂緒浩，弓 馨

(1.東北師范大學(xué)人文學(xué)院 理工學(xué)院，吉林 長春 130117；2.吉林省高校汽車電子技術(shù)工程研究中心，吉林 長春 130011)

1 引 言

在臨床醫(yī)學(xué)中，通過不同的醫(yī)學(xué)成像設(shè)備得到各種不同的醫(yī)學(xué)圖像。單一成像設(shè)備所得圖像具有局限性且不利于醫(yī)生診斷，所以醫(yī)學(xué)圖像融合成為近年來的研究熱點。通過將不同醫(yī)學(xué)圖像進行融合，即得到細節(jié)信息量豐富、冗余信息少的有效融合圖像[1]。

圖像融合主要分為3大層次：像素級、特征級和決策級。本文主要針對像素級的圖像融合技術(shù)進行研究。早期的像素級圖像融合技術(shù)主要是在拉普拉斯金字塔變換和小波變換的基礎(chǔ)上進行的[2～5]，但是拉普拉斯金字塔變換的融合方法會產(chǎn)生大量的冗余信息；小波變換在分解方向上只有固定的幾個方向，具有局限性，所以近年來提出了許多多分辨率分析的多尺度變換方法，并已經(jīng)廣泛地應(yīng)用到圖像融合技術(shù)中。比較熱門的多尺度方法有2006年提出的NSCT(Non-Subsampled Contourlet Transform)[6-7]和2010年提出的NSST(Non-Subsampled Shearlet Transform)[8～11]。這兩種變換方法可以提取圖像更全面的方向信息和細節(jié)信息，減少冗余量，至今仍是很多學(xué)者研究的熱點。但是這類的變換方法在進行圖像處理的同時易產(chǎn)生人造紋理，而且運行時間相對較長，這是我們所不期望的。

2013年何凱明提出了導(dǎo)向濾波器GF(Guided Filter)[12]，該濾波器具有邊緣保持平滑作用，是一種線性濾波器，處理圖像速度相比多尺度變換要快。導(dǎo)向濾波器在對圖像進行濾波的同時，可以增強圖像的空間連續(xù)性，避免引入人造紋理[1]。隨后許多學(xué)者在GF濾波器的基礎(chǔ)上，提出了新的圖像融合算法，減少人造紋理的同時，得到較好的融合效果，如GF-CCT-SML[13]，CST-GF-FL[14]等。

本文將導(dǎo)向濾波器應(yīng)用到醫(yī)學(xué)圖像融合技術(shù)中，對醫(yī)學(xué)圖像進行處理，利用GF濾波器的特性提取原圖像的細節(jié)層信息；再采用基于含可變參數(shù)p的鄰域統(tǒng)計特性融合規(guī)則進行圖像融合。利用該方法可以有效提取圖像的細節(jié)信息，增強圖像的空間連續(xù)性，由于引入可變參數(shù)p，其融合效果具有可調(diào)節(jié)性，更易滿足不同類別的醫(yī)學(xué)圖像融合技術(shù)的要求。

2 基于GF濾波器的醫(yī)學(xué)圖像融合算法

2.1 導(dǎo)向濾波器

導(dǎo)向濾波器(GF)[12]是一種線性濾波器，具有邊緣保持平滑作用，主要包括輸入圖像I，導(dǎo)向圖像G及輸出圖像O。導(dǎo)向濾波器是基于局部線性模型得到的，假設(shè)在以一個像素點k為中心的局部窗口wk中，輸出圖像O是導(dǎo)向圖像G的線性變換，具體公式如下：

Oj=akGj+bk，?j∈wk

，

(1)

其中：wk是一個以半徑為r的方窗，其大小即為(2r+1)×(2r+1)。ak和bk是局部窗wk內(nèi)的常數(shù)，可由最小化下面代價函數(shù)估計出：

，

(2)

其中：ε為正則化參數(shù)，系數(shù)ak和bk可以由線性回歸求出：

(3)

，

(4)

由公式(1)可知，當窗wk在變化時，濾波輸出圖像Oj也隨之變化。所以就把系數(shù)ak和bk進行取均值來解決，具體公式如下：

，

(5)

導(dǎo)向濾波器由公式(3)、(4)和(5)構(gòu)成，當輸入圖像I和導(dǎo)向圖像G相同時，GF濾波器表現(xiàn)出與雙邊濾波相似的保邊特性。

將GF濾波器的輸入輸出關(guān)系表示為

O=GF(I,G)

.

(6)

2.2 融合框架

本文提出一種基于GF濾波器的醫(yī)學(xué)圖像融合算法，主要利用GF濾波器可以實現(xiàn)輸出圖像對導(dǎo)向圖像的線性變換。當輸入圖像I與導(dǎo)向圖像G不同時，濾波輸出的圖像在細節(jié)上與輸入圖像相似，在線性結(jié)構(gòu)和紋理上又與導(dǎo)向圖像G相似。具體的融合框圖如圖1所示。

圖1 圖像融合算法的實現(xiàn)框圖Fig.1 Realization block diagram of image fusion algorithm

具體步驟如下：

Step1：將兩幅醫(yī)學(xué)圖像X、Y分別當作輸入圖像I和導(dǎo)向圖像G，進行GF濾波后得到一幅基礎(chǔ)層信息圖像O1=GF(X,Y);

Step2：將輸入圖像X與輸出圖像O1作差，得到一幅細節(jié)層信息圖像，即為X-GF(X,Y);

Step3：將兩幅醫(yī)學(xué)圖像Y、X分別當做輸入圖像I和導(dǎo)向圖像G，進行GF濾波后得到一幅基礎(chǔ)層信息圖像O2=GF(Y,X);

Step4：將輸入圖像Y與輸出圖像O2作差，得到一幅細節(jié)層信息圖像，即為Y-GF(Y,X);

Step5：按照融合規(guī)則，分別對兩幅細節(jié)層信息圖像求解權(quán)值系數(shù)；

Step6：對兩幅原圖像按照所求得的權(quán)值系數(shù)進行圖像融合，得到最終的融合圖像。

2.3 融合規(guī)則

融合的目的就是最大程度地提取原圖像的有效信息，使得融合后的圖像細節(jié)信息豐富，對于醫(yī)學(xué)圖像來說，融合圖像不僅要含有豐富的細節(jié)信息，而且應(yīng)該便于醫(yī)生發(fā)現(xiàn)患者的病灶，根據(jù)MRI圖像和CT圖像的特點，提出一種增強融合圖像細節(jié)信息的融合規(guī)則，具體如下：

數(shù)學(xué)統(tǒng)計特性已經(jīng)被應(yīng)用到圖像融合規(guī)則中，用來求得融合系數(shù)的權(quán)值[15～17]。對于一個矩陣X，其協(xié)方差矩陣COV(X)定義為

COV(X)=E[(X-E[X])(X-E[X])T]

，

(7)

在數(shù)字圖像中，我們選取以某一像素點(x,y)為中心的局部方窗w，其大小為m×m。將局部方窗假設(shè)為矩陣Z，那么在像素點(x,y)處的協(xié)方差矩陣為

，

(8)

為了增強融合效果，在原有的求解協(xié)方差矩陣中加入可變參數(shù)p，具體公式如下：

，

(9)

其中，p≥1,當p=1時，即為原有的協(xié)方差矩陣。

，

(10)

，

(11)

再對αH(x,y)和αV(x,y)求和，即可得到權(quán)重系數(shù)W(x,y)：

W(x,y)=αH(x,y)+αV(x,y)

.

(12)

將窗口滑動經(jīng)過整幅圖像，既可以求得在每一個像素點處的權(quán)值系數(shù)。

在圖1所示的融合框圖中，分別對兩幅原圖像按照上述方法進行權(quán)值系數(shù)求解，分別求得WX(i,j)和WY(i,j),即可求得最終的融合圖像F(i,j)：

(13)

3 圖像融合仿真實驗及分析

實驗平臺為Matlab2014(a)，計算機主機配置為:Intel (R)處理器，CPU主頻3.4 GHz，內(nèi)存8 GB。

本文選擇正常人類頭部的CT和MRI圖像進行仿真實驗，圖像大小為256×256，如圖2所示。

為了驗證本文所提的圖像融合方法的有效性，選擇3組仿真實驗: (1)本文圖像融合方法的實現(xiàn)：將本文所提出的方法進行實驗仿真，按照融合框架，給出融合過程中的每一步實驗結(jié)果，實現(xiàn)醫(yī)學(xué)圖像融合；(2)圖像融合方法對比實驗：選擇不同的融合方法分別進行醫(yī)學(xué)圖像融合，給出融合圖像對比圖，再選擇合適的客觀評價指標進行對比分析；(3)客觀質(zhì)量評價指標隨參數(shù)p變化關(guān)系實驗：選擇不同參數(shù)p，根據(jù)本文方法進行仿真，得到不同的融合對比結(jié)果，并進行分析。

3.1 圖像融合方法的實現(xiàn)

將本文所提出的基于導(dǎo)向濾波器的醫(yī)學(xué)圖像融合方法通過仿真實驗展示其具體的運算過程結(jié)果。其中，GF濾波器采用matlab2014(a)中imguidedfilter函數(shù)進行濾波處理，參數(shù)選擇：鄰域半徑r=25，平滑系數(shù)ε=2.1[17]；在計算權(quán)重系數(shù)時選擇方窗w大小為5×5，參數(shù)p=2。按照圖1所示的融合框架進行圖像融合，具體結(jié)果見圖3。

圖2 原圖像Fig.2 Source images

圖3 圖像融合過程仿真結(jié)果Fig.3 Simulation results of image fusion process

3.2 圖像融合方法對比實驗

選擇基于小波融合[5]、基于NSCT融合[6-7]、NSST+NMF[11]、GF+COVW[17]及本文所提算法分別進行醫(yī)學(xué)圖像融合。其中，小波變換選擇‘Db4’小波，分解層數(shù)為5；NSCT變換分解層數(shù)選為2；NSST變換分解層數(shù)設(shè)置為3，方向數(shù)分別選為4,8,16；GF+COVW融合方法中GF的參數(shù)：鄰域半徑r=25，平滑系數(shù)ε=2.1，w=5；本文方法參數(shù)：GF濾波器的鄰域半徑r=25，平滑系數(shù)ε=2.1，w=5，參數(shù)p=5。給出各融合方法的對比圖，如圖4所示；再選擇參數(shù)熵值(E)、平均梯度(AG)、標準差(SD)、互信息(MI)和空間頻率(SF)和均值(μ)共6個參數(shù)作為客觀評價指標進行對比分析，如表1所示。其中，空間頻率的大小可以反映融合圖像紋理的細致與粗糙程度，SF的值越大，說明圖像的紋理越細致，融合效果越好[18]。

通過觀察圖4，從視覺效果的角度進行對比，可以發(fā)現(xiàn)，采用GF濾波器的方法明顯優(yōu)于其他方法，而GF+COVW與本文方法相比較，后者在邊緣處有視覺效果增強，這樣更便于醫(yī)生明顯地發(fā)現(xiàn)患者的病灶，使得醫(yī)學(xué)圖像融合的視覺效果更佳；通過表1的客觀評價指標的對比，NSST+NMF的均值最大，NSCT的標準差最大，NSST+NMF的熵值最大，本文的平均梯度、空間頻率和互信息都是最大的，綜合對比發(fā)現(xiàn)，本文的方法在綜合性能上要明顯優(yōu)于其他方法。

圖4 圖像融合實驗對比圖Fig.4 Experimental comparison of image fusion

表1 客觀評價指標Tab.1 Objective quality evaluation of different fusion methods

3.3 客觀質(zhì)量評價參數(shù)與參數(shù)p變化關(guān)系實驗

為了得到融合效果隨參數(shù)p變化的關(guān)系，選擇不同的參數(shù)p按照本文所給出的融合框架進行圖像融合，并將其融合圖像的客觀質(zhì)量評價指標進行對比分析，如表2所示。為了清楚地觀察到各評價指標與參數(shù)p之間的變化趨勢，分別畫出各個參數(shù)隨p變化的曲線，詳見表2和圖5。

通過觀察表2和圖5可知，隨著參數(shù)p的增大，各評價參數(shù)值整體呈現(xiàn)增加的趨勢。其中均值μ和標準差SD在p>3后曲線變緩，p>10后基本趨于穩(wěn)定；熵值MSG是先變大后變小，而后又趨于平穩(wěn)的趨勢，在p=3左右時達到最大；空間頻率SF曲線隨著p的增大，一直呈現(xiàn)上升趨勢；平均梯度AG和互信息MI在p<3時，曲線較陡，上升幅度明顯，在p>3后曲線上升趨勢變緩；在p>10時后基本趨于穩(wěn)定。綜上所述，本文提出的醫(yī)學(xué)圖像融合方法含有可變參數(shù)p，融合效果也隨著p的變化而變化。隨著p的增大，融合效果明顯變好，而p值的最優(yōu)解問題將在后續(xù)的工作中繼續(xù)進行研究。

表2 隨參數(shù)p變化的客觀評價參數(shù)值Tab.2 Objective evaluate parameter values changing with the parameter p

圖5 客觀評價指標隨p變化的曲線圖Fig.5 Graph of objective evaluation index changing withp

4 結(jié) 論

導(dǎo)向濾波器是一種線性平滑濾波器，具有良好的邊緣保持和細節(jié)增強的特性，將GF應(yīng)用到醫(yī)學(xué)圖像融合技術(shù)中，可以有效增強融合圖像的空間連續(xù)性，減少人造紋理的產(chǎn)生，對圖像的融合效果有明顯的提高。本文給出了GF濾波器的理論基礎(chǔ)和融合框架，首先，利用GF濾波器的細節(jié)增強特性，提取圖像的基礎(chǔ)層信息，再經(jīng)過與原圖像的作差運算，即得細節(jié)層信息；再通過融合規(guī)則進行圖像融合。本文提出一種基于可變參數(shù)p的鄰域數(shù)學(xué)統(tǒng)計量的權(quán)值計算方法，在原有的協(xié)方差矩陣求解中，加入?yún)?shù)p,當p>1時，明顯增強了融合圖像的細節(jié)信息，經(jīng)過實驗仿真結(jié)果對比和分析，隨著參數(shù)p變大，融合效果明顯提高；當p>10后，除空間頻率SF外，其他客觀評價指標都已趨于穩(wěn)定。該方法可以有效實現(xiàn)醫(yī)學(xué)圖像融合，滿足醫(yī)學(xué)圖像融合的技術(shù)指標。