郭 盼,何文超,梁龍凱,張 萌,呂緒浩,弓 馨
(1.東北師范大學(xué)人文學(xué)院 理工學(xué)院,吉林 長春 130117;2.吉林省高校汽車電子技術(shù)工程研究中心,吉林 長春 130011)
在臨床醫(yī)學(xué)中,通過不同的醫(yī)學(xué)成像設(shè)備得到各種不同的醫(yī)學(xué)圖像。單一成像設(shè)備所得圖像具有局限性且不利于醫(yī)生診斷,所以醫(yī)學(xué)圖像融合成為近年來的研究熱點(diǎn)。通過將不同醫(yī)學(xué)圖像進(jìn)行融合,即得到細(xì)節(jié)信息量豐富、冗余信息少的有效融合圖像[1]。
圖像融合主要分為3大層次:像素級、特征級和決策級。本文主要針對像素級的圖像融合技術(shù)進(jìn)行研究。早期的像素級圖像融合技術(shù)主要是在拉普拉斯金字塔變換和小波變換的基礎(chǔ)上進(jìn)行的[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]。這兩種變換方法可以提取圖像更全面的方向信息和細(xì)節(jié)信息,減少冗余量,至今仍是很多學(xué)者研究的熱點(diǎn)。但是這類的變換方法在進(jìn)行圖像處理的同時易產(chǎn)生人造紋理,而且運(yùn)行時間相對較長,這是我們所不期望的。
2013年何凱明提出了導(dǎo)向?yàn)V波器GF(Guided Filter)[12],該濾波器具有邊緣保持平滑作用,是一種線性濾波器,處理圖像速度相比多尺度變換要快。導(dǎo)向?yàn)V波器在對圖像進(jìn)行濾波的同時,可以增強(qiáng)圖像的空間連續(xù)性,避免引入人造紋理[1]。隨后許多學(xué)者在GF濾波器的基礎(chǔ)上,提出了新的圖像融合算法,減少人造紋理的同時,得到較好的融合效果,如GF-CCT-SML[13],CST-GF-FL[14]等。
本文將導(dǎo)向?yàn)V波器應(yīng)用到醫(yī)學(xué)圖像融合技術(shù)中,對醫(yī)學(xué)圖像進(jìn)行處理,利用GF濾波器的特性提取原圖像的細(xì)節(jié)層信息;再采用基于含可變參數(shù)p的鄰域統(tǒng)計特性融合規(guī)則進(jìn)行圖像融合。利用該方法可以有效提取圖像的細(xì)節(jié)信息,增強(qiáng)圖像的空間連續(xù)性,由于引入可變參數(shù)p,其融合效果具有可調(diào)節(jié)性,更易滿足不同類別的醫(yī)學(xué)圖像融合技術(shù)的要求。
導(dǎo)向?yàn)V波器(GF)[12]是一種線性濾波器,具有邊緣保持平滑作用,主要包括輸入圖像I,導(dǎo)向圖像G及輸出圖像O。導(dǎo)向?yàn)V波器是基于局部線性模型得到的,假設(shè)在以一個像素點(diǎn)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)可知,當(dāng)窗wk在變化時,濾波輸出圖像Oj也隨之變化。所以就把系數(shù)ak和bk進(jìn)行取均值來解決,具體公式如下:
,
(5)
導(dǎo)向?yàn)V波器由公式(3)、(4)和(5)構(gòu)成,當(dāng)輸入圖像I和導(dǎo)向圖像G相同時,GF濾波器表現(xiàn)出與雙邊濾波相似的保邊特性。
將GF濾波器的輸入輸出關(guān)系表示為
O=GF(I,G)
.
(6)
本文提出一種基于GF濾波器的醫(yī)學(xué)圖像融合算法,主要利用GF濾波器可以實(shí)現(xiàn)輸出圖像對導(dǎo)向圖像的線性變換。當(dāng)輸入圖像I與導(dǎo)向圖像G不同時,濾波輸出的圖像在細(xì)節(jié)上與輸入圖像相似,在線性結(jié)構(gòu)和紋理上又與導(dǎo)向圖像G相似。具體的融合框圖如圖1所示。
圖1 圖像融合算法的實(shí)現(xiàn)框圖Fig.1 Realization block diagram of image fusion algorithm
具體步驟如下:
Step1:將兩幅醫(yī)學(xué)圖像X、Y分別當(dāng)作輸入圖像I和導(dǎo)向圖像G,進(jìn)行GF濾波后得到一幅基礎(chǔ)層信息圖像O1=GF(X,Y);
Step2:將輸入圖像X與輸出圖像O1作差,得到一幅細(xì)節(jié)層信息圖像,即為X-GF(X,Y);
Step3:將兩幅醫(yī)學(xué)圖像Y、X分別當(dāng)做輸入圖像I和導(dǎo)向圖像G,進(jìn)行GF濾波后得到一幅基礎(chǔ)層信息圖像O2=GF(Y,X);
Step4:將輸入圖像Y與輸出圖像O2作差,得到一幅細(xì)節(jié)層信息圖像,即為Y-GF(Y,X);
Step5:按照融合規(guī)則,分別對兩幅細(xì)節(jié)層信息圖像求解權(quán)值系數(shù);
Step6:對兩幅原圖像按照所求得的權(quán)值系數(shù)進(jìn)行圖像融合,得到最終的融合圖像。
融合的目的就是最大程度地提取原圖像的有效信息,使得融合后的圖像細(xì)節(jié)信息豐富,對于醫(yī)學(xué)圖像來說,融合圖像不僅要含有豐富的細(xì)節(jié)信息,而且應(yīng)該便于醫(yī)生發(fā)現(xiàn)患者的病灶,根據(jù)MRI圖像和CT圖像的特點(diǎn),提出一種增強(qiáng)融合圖像細(xì)節(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ù)字圖像中,我們選取以某一像素點(diǎn)(x,y)為中心的局部方窗w,其大小為m×m。將局部方窗假設(shè)為矩陣Z,那么在像素點(diǎn)(x,y)處的協(xié)方差矩陣為
,
(8)
為了增強(qiáng)融合效果,在原有的求解協(xié)方差矩陣中加入可變參數(shù)p,具體公式如下:
,
(9)
其中,p≥1,當(dāng)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)過整幅圖像,既可以求得在每一個像素點(diǎn)處的權(quán)值系數(shù)。
在圖1所示的融合框圖中,分別對兩幅原圖像按照上述方法進(jìn)行權(quán)值系數(shù)求解,分別求得WX(i,j)和WY(i,j),即可求得最終的融合圖像F(i,j):
(13)
實(shí)驗(yàn)平臺為Matlab2014(a),計算機(jī)主機(jī)配置為:Intel (R)處理器,CPU主頻3.4 GHz,內(nèi)存8 GB。
本文選擇正常人類頭部的CT和MRI圖像進(jìn)行仿真實(shí)驗(yàn),圖像大小為256×256,如圖2所示。
為了驗(yàn)證本文所提的圖像融合方法的有效性,選擇3組仿真實(shí)驗(yàn): (1)本文圖像融合方法的實(shí)現(xiàn):將本文所提出的方法進(jìn)行實(shí)驗(yàn)仿真,按照融合框架,給出融合過程中的每一步實(shí)驗(yàn)結(jié)果,實(shí)現(xiàn)醫(yī)學(xué)圖像融合;(2)圖像融合方法對比實(shí)驗(yàn):選擇不同的融合方法分別進(jìn)行醫(yī)學(xué)圖像融合,給出融合圖像對比圖,再選擇合適的客觀評價指標(biāo)進(jìn)行對比分析;(3)客觀質(zhì)量評價指標(biāo)隨參數(shù)p變化關(guān)系實(shí)驗(yàn):選擇不同參數(shù)p,根據(jù)本文方法進(jìn)行仿真,得到不同的融合對比結(jié)果,并進(jìn)行分析。
將本文所提出的基于導(dǎo)向?yàn)V波器的醫(yī)學(xué)圖像融合方法通過仿真實(shí)驗(yàn)展示其具體的運(yùn)算過程結(jié)果。其中,GF濾波器采用matlab2014(a)中imguidedfilter函數(shù)進(jìn)行濾波處理,參數(shù)選擇:鄰域半徑r=25,平滑系數(shù)ε=2.1[17];在計算權(quán)重系數(shù)時選擇方窗w大小為5×5,參數(shù)p=2。按照圖1所示的融合框架進(jìn)行圖像融合,具體結(jié)果見圖3。
圖2 原圖像Fig.2 Source images
圖3 圖像融合過程仿真結(jié)果Fig.3 Simulation results of image fusion process
選擇基于小波融合[5]、基于NSCT融合[6-7]、NSST+NMF[11]、GF+COVW[17]及本文所提算法分別進(jìn)行醫(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)、標(biāo)準(zhǔn)差(SD)、互信息(MI)和空間頻率(SF)和均值(μ)共6個參數(shù)作為客觀評價指標(biāo)進(jìn)行對比分析,如表1所示。其中,空間頻率的大小可以反映融合圖像紋理的細(xì)致與粗糙程度,SF的值越大,說明圖像的紋理越細(xì)致,融合效果越好[18]。
通過觀察圖4,從視覺效果的角度進(jìn)行對比,可以發(fā)現(xiàn),采用GF濾波器的方法明顯優(yōu)于其他方法,而GF+COVW與本文方法相比較,后者在邊緣處有視覺效果增強(qiáng),這樣更便于醫(yī)生明顯地發(fā)現(xiàn)患者的病灶,使得醫(yī)學(xué)圖像融合的視覺效果更佳;通過表1的客觀評價指標(biāo)的對比,NSST+NMF的均值最大,NSCT的標(biāo)準(zhǔn)差最大,NSST+NMF的熵值最大,本文的平均梯度、空間頻率和互信息都是最大的,綜合對比發(fā)現(xiàn),本文的方法在綜合性能上要明顯優(yōu)于其他方法。
圖4 圖像融合實(shí)驗(yàn)對比圖Fig.4 Experimental comparison of image fusion
表1 客觀評價指標(biāo)Tab.1 Objective quality evaluation of different fusion methods
為了得到融合效果隨參數(shù)p變化的關(guān)系,選擇不同的參數(shù)p按照本文所給出的融合框架進(jìn)行圖像融合,并將其融合圖像的客觀質(zhì)量評價指標(biāo)進(jìn)行對比分析,如表2所示。為了清楚地觀察到各評價指標(biāo)與參數(shù)p之間的變化趨勢,分別畫出各個參數(shù)隨p變化的曲線,詳見表2和圖5。
通過觀察表2和圖5可知,隨著參數(shù)p的增大,各評價參數(shù)值整體呈現(xiàn)增加的趨勢。其中均值μ和標(biāo)準(zhǔn)差SD在p>3后曲線變緩,p>10后基本趨于穩(wěn)定;熵值MSG是先變大后變小,而后又趨于平穩(wěn)的趨勢,在p=3左右時達(dá)到最大;空間頻率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ù)進(jìn)行研究。
表2 隨參數(shù)p變化的客觀評價參數(shù)值Tab.2 Objective evaluate parameter values changing with the parameter p
圖5 客觀評價指標(biāo)隨p變化的曲線圖Fig.5 Graph of objective evaluation index changing withp
導(dǎo)向?yàn)V波器是一種線性平滑濾波器,具有良好的邊緣保持和細(xì)節(jié)增強(qiáng)的特性,將GF應(yīng)用到醫(yī)學(xué)圖像融合技術(shù)中,可以有效增強(qiáng)融合圖像的空間連續(xù)性,減少人造紋理的產(chǎn)生,對圖像的融合效果有明顯的提高。本文給出了GF濾波器的理論基礎(chǔ)和融合框架,首先,利用GF濾波器的細(xì)節(jié)增強(qiáng)特性,提取圖像的基礎(chǔ)層信息,再經(jīng)過與原圖像的作差運(yùn)算,即得細(xì)節(jié)層信息;再通過融合規(guī)則進(jìn)行圖像融合。本文提出一種基于可變參數(shù)p的鄰域數(shù)學(xué)統(tǒng)計量的權(quán)值計算方法,在原有的協(xié)方差矩陣求解中,加入?yún)?shù)p,當(dāng)p>1時,明顯增強(qiáng)了融合圖像的細(xì)節(jié)信息,經(jīng)過實(shí)驗(yàn)仿真結(jié)果對比和分析,隨著參數(shù)p變大,融合效果明顯提高;當(dāng)p>10后,除空間頻率SF外,其他客觀評價指標(biāo)都已趨于穩(wěn)定。該方法可以有效實(shí)現(xiàn)醫(yī)學(xué)圖像融合,滿足醫(yī)學(xué)圖像融合的技術(shù)指標(biāo)。