王雷 李彬 田聯(lián)房

(華南理工大學(xué)自動化科學(xué)與工程學(xué)院，廣東廣州510640)

基于多尺度分解的醫(yī)學(xué)圖像融合方法的核心問題之一是多尺度分解工具的選擇.小波變換是最典型的多尺度分解工具［1］.然而，二維離散小波變換只能將圖像分解為水平、垂直、對角3個高頻子帶，這使其只能捕獲有限的方向信息，因此，基于小波變換的融合方法由于小波不能很好地表示圖像中的線、邊緣、輪廓等高維信息而會導(dǎo)致偽影現(xiàn)象.剪切波是為克服小波的缺點而提出的一種新穎的多尺度幾何分析工具.相比圖像融合中常用的其它多尺度分解工具，如曲線波變換［2］、輪廓波變換［3］、非下采樣輪廓波變換 (NSCT)［4］等，剪切波變換具有其自身的優(yōu)勢.從近似理論的觀點來看，剪切波形成一個具有小波形式的多尺度多方向的緊支集結(jié)構(gòu)，它是圖像信號(如邊緣)的真正的二維稀疏表示，目前只有曲線波具有類似的性質(zhì).然而，曲線波的實現(xiàn)過程過于復(fù)雜，并且不能在圖像的多分辨分析框架內(nèi)實現(xiàn).相比輪廓波和NSCT，雖然剪切波變換具有相似的實現(xiàn)過程，但是剪切波在剪切過程中沒有方向數(shù)目和支撐基尺寸大小的限制.此外，剪切波逆變換只需對正向變換中的剪切濾波器進(jìn)行加和處理，而不需要像輪廓波變換和NSCT那樣對方向濾波器逆合成，這說明剪切波的實現(xiàn)過程具有更高的計算效率.目前，剪切波變換已經(jīng)被應(yīng)用于遙感圖像的融合中，并取得了良好的效果［5-6］.然而，在上述文獻(xiàn)中，剪切波變換的離散化過程通過下采樣策略實現(xiàn)，由于其不具備平移不變的性質(zhì)而易產(chǎn)生偽吉布斯現(xiàn)象.

基于多尺度分解的醫(yī)學(xué)圖像融合方法的另一個核心問題是頻域系數(shù)的融合規(guī)則.待融合圖像進(jìn)行多尺度分解后，低頻系數(shù)表示圖像的近似逼近，高頻系數(shù)代表圖像的細(xì)節(jié)信息.傳統(tǒng)的低頻取平均策略易降低融合圖像的對比度.對于高頻融合，基于人工神經(jīng)網(wǎng)絡(luò)的融合策略近年來越來越受到學(xué)者們的重視，如脈沖耦合神經(jīng)網(wǎng)絡(luò)［7-8］.但是在這些應(yīng)用中，神經(jīng)網(wǎng)絡(luò)包含許多參數(shù)，需要專家依靠經(jīng)驗去設(shè)置，這限制了它們的使用范圍.自生成神經(jīng)網(wǎng)絡(luò)(SGNN)是一類特點鮮明的自組織神經(jīng)網(wǎng)絡(luò)，它不需要外在參數(shù)的設(shè)定，因而具有自主性高的優(yōu)點［9-10］.但是上述文獻(xiàn)中的基于SGNN的圖像融合方法是直接對圖像的像素進(jìn)行聚類，只考慮了圖像灰度值的變化，其本質(zhì)是將待融合圖像像素的灰度值聚類后進(jìn)行加權(quán)處理，因而可能導(dǎo)致融合后的圖像產(chǎn)生塊狀效應(yīng).基于上述分析，文中提出一種基于平移不變剪切波變換的醫(yī)學(xué)圖像融合新方法.

1 基于平移不變剪切波變換的融合方法

文中提出的基于平移不變剪切波變換的醫(yī)學(xué)圖像融合方法主要包括3步:(1)將圖像進(jìn)行多尺度的平移不變剪切波分解，獲得圖像的頻率子帶;(2)對低頻系數(shù)采用基于區(qū)域系數(shù)絕對值和權(quán)重的融合規(guī)則進(jìn)行處理;對于高頻系數(shù)，在每個高頻子帶，計算其支持向量值，將支持向量值輸入到SGNN中進(jìn)行聚類分析，根據(jù)聚類結(jié)果確定融合后圖像的高頻系數(shù);(3)通過逆平移不變剪切波變換重建融合后的圖像.文中待融合的圖像記為A和B，融合后的圖像記為F，具體過程如圖1所示.

圖1 基于平移不變剪切波變換的醫(yī)學(xué)圖像融合算法示意圖Fig.1 Schematic diagram of medical image fusion algorithm based on shift-invariant shearlet transformation

1.1 圖像的平移不變剪切波變換

式中:ψ∈L2(R2)，L表示可積空間，det表示矩陣的行列式;A稱為各向異性膨脹矩陣，A=

平移不變剪切波變換［11］的離散化過程主要分為兩步:多尺度剖分和方向局部化.多尺度剖分是通過非下采樣的金字塔濾波器組［12］(NSPF)實現(xiàn)的，每一級需要對上級中所采用的濾波器按矩陣D=進(jìn)行上采樣，圖像經(jīng)過d級非下采樣塔式濾波，可以得到d+1個與源圖像大小相同的子帶圖像.

方向局部化是通過剪切濾波器(SF)完成的.標(biāo)準(zhǔn)的剪切波變換中使用的剪切濾波器是在偽極化網(wǎng)格中通過窗函數(shù)的平移操作實現(xiàn)的，其過程需要下采樣，因而不具有平移不變性.平移不變的剪切變換把標(biāo)準(zhǔn)的剪切濾波器從偽極化網(wǎng)格系統(tǒng)映射回到笛卡爾坐標(biāo)系統(tǒng).通過逆傅里葉變換，可以證明其操作可以直接通過二維卷積完成(參看文獻(xiàn)［11］定理3.1)，避免了下采樣操作因而具有平移不變性.其具體過程又稱為Meyer窗口的構(gòu)建.

圖2 離散的平移不變剪切波變換示意圖Fig.2 Schematic diagram of discrete shift-invariant shearlet transformation

令j表示圖像分解的尺度，j=1，2，…，M，平移不變剪切波變換的具體過程如圖2所示，可描述為:

(1)利用非下采樣金字塔策略將圖像fj分解成為低通圖像fj+1和細(xì)節(jié)圖像gj+1;

(2)針對細(xì)節(jié)圖像gj+1，構(gòu)建Meyer窗口進(jìn)行多尺度剖分;

(3)對各個方向子帶進(jìn)行逆傅里葉變換得到平移不變剪切波系數(shù)cj+1.

剪切波滿足拋物線尺度，對某尺度子帶圖像進(jìn)行l(wèi)級平移不變剪切波方向分解，可得到2l+2個與原始輸入圖像尺寸大小相同的子帶圖像.因此，對圖像j級分解后，可得到個子帶圖像，其中，J表示圖像分解層數(shù)，lj是在尺度j下的方向分解級數(shù).圖3給出了圖像Lena分解成2層的實例.

圖3 圖像Lena的兩層離散的平移不變剪切波分解示例Fig.3 Illustration of decomposing image Lena into 2 levels by discrete shift-invariant shearlet transformation

經(jīng)過平移不變剪切波分解后，每一個剪切波位于近似大小為22j×2j、線方向斜率為2-jl的梯形對上，如圖4(a)所示.根據(jù)小波理論，由一維小波張成的二維小波基具有正方形的支撐區(qū)間，在尺度j，其近似邊長為2j，如圖4(b)所示.當(dāng)尺度j變大時，小波支撐基的面積隨之變大，但是由于其無方向性，非零的小波系數(shù)以指數(shù)形式增長，導(dǎo)致大量不可忽略的系數(shù)出現(xiàn).小波逼近奇異曲線最終表現(xiàn)為圖5(a)所示的“點”逼近過程，它不再是原曲線最優(yōu)的稀疏表示.剪切波以圖5(b)所示的近似長條形即梯形對去逼近曲線，當(dāng)尺度j變大時，其支撐基的線方向斜率2-jl也隨之變化，使得這種近似長條形的支撐基具有各向異性，能夠有效捕獲方向信息，它是圖像中邊緣等特征的真正的稀疏表示.這一性質(zhì)對于醫(yī)學(xué)圖像中病變部位的表示、提取具有重要意義.

圖4 剪切波與小波的支撐基比較Fig.4 Comparison of the support bases for shearlet and wavelet

圖5剪切波與小波逼近曲線Fig.5 Curves for shearlet and wavelet representations

1.2 低頻系數(shù)融合

為避免降低融合圖像的對比度，文中采用如下基于區(qū)域系數(shù)絕對值和權(quán)重的融合規(guī)則.

(1)計算低頻系數(shù)在3×3鄰域的絕對值的和:

式中，Cp(i，j)表示圖像p在(i，j)處的低頻系數(shù).

(2)計算系數(shù)的權(quán)重:

(3)計算融合圖像在(i，j)處的低頻系數(shù):

1.3 高頻系數(shù)融合

SGNN是自組織神經(jīng)網(wǎng)絡(luò)的一個擴(kuò)展，樣本通過非監(jiān)督學(xué)習(xí)形成一棵包括神經(jīng)節(jié)點、節(jié)點間聯(lián)系、節(jié)點權(quán)重的自生成神經(jīng)樹(SGNT)，具有自主性高、學(xué)習(xí)速度快的優(yōu)點.為了充分利用SGNN的優(yōu)點在高頻子帶獲取更多的特征信息，文中將支持向量機(jī)引入SGNN，提出一種基于支持向量值激勵的SGNN高頻子帶融合策略.把每個高頻子帶系數(shù)的支持向量值(T)作為其高頻特征，支持向量值的大小反映了高頻特征的顯著性.令表示圖像分解后第l層、第k方向位置為(i，j)的高頻子帶系數(shù)，基于支持向量值激勵的SGNN高頻子帶融合方法的具體步驟如下.

(1)計算兩幅圖像每個高頻子帶的支持向量值:

式中，V為支持向量值濾波器，其取值參見文獻(xiàn)［13］.

(2)將支持向量值輸入到SGNN中進(jìn)行聚類分析，生成SGNT:

式中，c是以Nnode為根節(jié)點的所有葉子節(jié)點的個數(shù);

(3)對產(chǎn)生的SGNT進(jìn)行水平修剪和垂直修剪.

(4)聚類結(jié)果模糊化:由于圖像的差異，在兩幅圖像高頻子帶的同一位置的支持向量值可能屬于不同的類，因此引入模糊隸屬度函數(shù)，每一個支持向量值不再屬于特定的某一類，而是以一定的模糊隸屬度值屬于某一類，則在位置(m，n)處的支持向量值q(m，n)屬于第j類的隸屬度函數(shù)為

式中，wj是第j類的權(quán)重值，r為分類的總數(shù).

(5)融合后圖像的高頻子帶系數(shù)為

2 實驗結(jié)果及分析

文中通過3組實驗驗證所提方法的融合性能.實驗1比較了文獻(xiàn)［5］中所用的剪切波變換與文中提出的基于平移不變剪切波變換的方法在相同融合規(guī)則下的融合性能，結(jié)果表明文中方法可以有效地克服偽吉布斯現(xiàn)象.實驗2通過對肺部的計算機(jī)斷層掃描(CT)圖像與正電子發(fā)射掃描(PET)圖像融合來驗證文中方法的有效性.實驗3在相同融合規(guī)則下比較了平移不變剪切波變換與NSCT的融合性能，結(jié)果表明平移不變剪切波具有很高的計算效率.實驗中平移不變剪切波分解層數(shù)為4，所有圖像的大小為256×256，并在融合前都已經(jīng)完全配準(zhǔn)對齊.實驗平臺為Matlab7.4，計算機(jī)主機(jī)配置為:Intel(R)處理器，CPU主頻1.0GHz，內(nèi)存512MB.

實驗中對SGNN、小波變換法(Wavelet)、輪廓波變換法(Contourlet)和NSCT 4種方法分別從熵(EN)、互信息(MI)、平均梯度(AvG)、QAB/F4 個方面［14］與文中方法進(jìn)行了比較.熵的大小反映了融合圖像包含的信息量的多少;互信息反映了融合圖像從源圖像獲取信息的多少;平均梯度表示圖像的清晰度，反映了圖像對微小細(xì)節(jié)反差表達(dá)的能力;QAB/F值越大表明融合后圖像的邊緣信息越豐富.在實驗中用Shearlet表示文獻(xiàn)［5］中基于剪切波變換的融合方法.

2.1 剪切波與平移不變剪切波融合結(jié)果比較

圖6(a)是核磁共振(MRI)圖像，其亮度與組織中氫原子數(shù)目有關(guān)，圖中顯示的是軟組織.圖6(b)是CT圖像，它的亮度與組織的密度有關(guān)，圖中清晰顯示的部分是骨架.圖6(c)與6(d)分別是使用剪切波與平移不變剪切波圖像融合結(jié)果.圖6(e)與6(f)是融合圖像圖6(c)和6(d)減去圖6(a)得到的殘差圖像.如圖6(b)所示，大矩形所標(biāo)記處是一個突出的尖點部位，在這個突出點周圍灰度值變化劇烈，可近似看作是一個奇異點.由于剪切波變換不具有平移不變性，在奇異點處會產(chǎn)生偽吉布斯現(xiàn)象，圖像在奇異點處的特征就不能很好地被捕獲，導(dǎo)致融合后的圖像會丟失部分特征信息，如6(e)中的矩形標(biāo)記所標(biāo)記.對比圖6(f)與6(e)可知，平移不變剪切波變換能夠有效抑制偽吉布斯現(xiàn)象，捕獲更多的原始圖像中的特征信息.此外，圖6中的小矩形標(biāo)記部分的差異是平移不變剪切波融合性能優(yōu)于剪切波的又一個有力例證.

圖6 剪切波變換與平移不變剪切波變換的融合結(jié)果比較Fig.6 Comparison of fusion results obtained by shearlet transformation and shift-invariant shearlet transformation

2.2 肺部CT-PET圖像融合

肺部CT-PET圖像融合結(jié)果如圖7和表1所示.圖7(a)是原始CT圖像，它顯示的是人肺部的組織形態(tài)結(jié)構(gòu)，圖7(b)是原始PET圖像，它描述了肺的對應(yīng)位置的功能狀況，圖7(c)-7(h)是各種方法的融合圖像.

通過對比融合結(jié)果可以發(fā)現(xiàn)，圖7(c)丟失了大部分的CT圖像的信息，融合圖像非常不清晰.這主要是因為SGNN融合方法本質(zhì)上是對兩幅待融合圖像的像素加權(quán)求和處理，而PET圖像的像素灰度值取值范圍是0～32767，CT圖像的像素灰度值取值范圍是0～4095.PET圖像的像素取值范圍遠(yuǎn)大于CT圖像像素的取值范圍，簡單的加權(quán)求和處理必然會降低圖像的對比度.圖7(d)邊緣出現(xiàn)模糊現(xiàn)象，這主要是由于小波變換只能獲取圖像的“點”奇異信息，對邊緣、輪廓等不能很好地表示造成的.圖7(e)的邊緣出現(xiàn)了明顯的偽影現(xiàn)象.由表1可知，就所選用的4個評價指標(biāo)而言，文中方法優(yōu)于Contourlet和NSCT.由Shearlet和文中方法的結(jié)果對比可知，文中所用的融合方法能夠獲取更多源圖像中的信息.造成這種差異的主要原因是文中方法所用的平移不變剪切波不僅可以克服偽吉布斯現(xiàn)象的產(chǎn)生，而且不同于文獻(xiàn)［5］中僅僅考慮單個系數(shù)的絕對值的融合規(guī)則，文中提出的融合策略充分考慮了高頻子帶的特征信息，利用支持向量值表征高頻特征的顯著性，用SGNN自動對支持向量值進(jìn)行聚類分析，減少了外在人為因素的干擾，有利于高頻特征信息的獲取，因而更加有效.

圖7 肺部CT-PET圖像融合結(jié)果Fig.7 Fusion results of the CT-PET image of lung

2.3 文中方法與NSCT融合性能的比較

近年來，NSCT由于能夠克服輪廓波變換不具備平移不變性的缺點而被廣泛應(yīng)用于圖像融合中，并且被證明是有效的，如文獻(xiàn)［4］與［15］.因此，將文中提出的基于平移不變剪切波方法的融合性能與NSCT進(jìn)行了比較.為方便起見，融合規(guī)則是低頻取平均，高頻取絕對值最大，這樣避免了過多的融合運算的時間消耗，突出了變換本身的時間消耗.實驗中，分解層數(shù)都設(shè)為4，平移不變剪切波和NSCT的尺度向量分別為(3，3，4，4) 和 (0，1，3，4)，NSCT 的多尺度分解和方向剖分濾波器類型分別為“9-7”濾波器和“pkva”濾波器.實驗結(jié)果如圖8和表2所示.

表1 CT-PET圖像融合結(jié)果量化比較Table 1 Quantitative comparison of fusion results for CT-PET image

圖8 平移不變剪切波變換與NSCT的融合圖像Fig.8 Fusion images of shift-invariant shearlet transformation and NSCT

表2 文中方法與NSCT融合結(jié)果的量化比較Table 2 Quantitative comparison of fusion results obtained by the proposed method and NSCT

由圖8(c)和8(d)可見，CT圖像中的形態(tài)信息和PET圖像中的功能區(qū)域都很好地顯示在融合后的圖像中.從表2也可知兩者具有相似的融合性能.但是，從時間消耗來看，NSCT方法所用時間是平移不變剪切波方法的2倍多，這說明雖然具有類似的平移不變性，但是NSCT方法是以高時間消耗為代價的.造成這種差異的主要原因是:NSCT的逆變換過程需要對方向濾波器組重新進(jìn)行逆合成，而平移不變剪切波則只需對正向變換的剪切濾波器進(jìn)行簡單的加和處理即可完成，具體的細(xì)節(jié)可以參見文獻(xiàn)［11］.雖然降低時間消耗可以通過快速發(fā)展的計算機(jī)硬件技術(shù)的改進(jìn)實現(xiàn)，但是平移不變剪切波方法以其高計算效率無疑為實踐中海量數(shù)據(jù)的融合提供了另一種選擇.

3 結(jié)語

基于平移不變剪切波變換，文中提出一種醫(yī)學(xué)圖像融合方法.相比于傳統(tǒng)的基于剪切波變換的融合方法，文中提出的方法可以有效地克服偽吉布斯現(xiàn)象，并具有很高的計算效率.相比于傳統(tǒng)的融合規(guī)則，文中所用的低頻子帶的融合規(guī)則充分考慮了醫(yī)學(xué)圖像的特點;而對高頻子帶融合而言，支持向量值激勵的自生成神經(jīng)網(wǎng)絡(luò)方法比直接對圖像進(jìn)行像素聚類分析可以提取更多的高頻特征信息，因而融合效果更好.對融合圖像的比較和融合結(jié)果的量化分析表明:平移不變剪切波不僅可以被用于醫(yī)學(xué)圖像融合，而且相比小波等其它常用的多尺度幾何分析工具，具有更大的優(yōu)越性.因此基于平移不變剪切波變換的醫(yī)學(xué)圖像融合方法具有廣闊的應(yīng)用前景.

［1］李彬，田聯(lián)房，陳萍，等.基于模塊化設(shè)計的三維醫(yī)學(xué)圖像體繪制方法［J］.華南理工大學(xué)學(xué)報:自然科學(xué)版，2006，34(1):77-81.Li Bin，Tian Lian-fang，Chen Ping，et al.Volume rendering of 3D medical image based on modularization［J］.Journal of South China University of Technology:Natural Science Edition，2006，34(1):77-81.

［2］Candues D，Donoho D L.New tight frames of curvelets and optimal representations of objects with piecewiseC2singularities［J］.Communications on Pure and Applied Mathematics，2004，57(2):219-266.

［3］Yang L，Guo B L，Ni W.Multimodality medical image fusion based on multiscale geometric analysis of contourlet transform ［J］.Neurocomputing，2008，72(3):203-211.

［4］Li T J，Wang Y Y.Biological image fusion using a NSCT based variable-weight method ［J］.Information Fusion，2011，12(2):85-92.

［5］Miao Q G，Shi C，Xu P，et al.A novel algorithm of image fusion using shearlets［J］.Optics Communications，2010，284(6):1540-1547.

［6］Deng C Z，Wang S Q，Chen X.Remote sensing images fusion algorithm base on shearlet transform［C］∥Proceedings of the 2009 International Conference on Environmen-tal Science and Information Application Technology.Banff:ACM，2009:451-454.

［7］Wang Z B，Ma Y D.Medical image fusion based on m-PCNN ［J］.Information Fusion，2008，9(2):176-185.

［8］Yang S M，Wang M，Lu Y，et al.Fusion of multiparametric SAR images based on SW-nonsubsampled contourlet and PCNN［J］.Signal Processing，2009，89(12):2596-2608.

［9］Qin Z，Bao F M.A novel image fusion method based on SGNN ［J］.Lecture Notes in Computer Science，2005，3497(821):747-752.

［10］Jiang H，Tian Y F.Fuzzy image fusion based on modified self-generating neural network ［J］.Expert Systems with Applications，2011，38(7):8515-8523.

［11］Easley G，Labate D，Lim W-Q.Sparse directional image representations using the discrete shearlet transform［J］.Applied and Computational Harmonic Analysis，2008，25(1):25-46.

［12］Da Cunha A L，Zhou J P，Do M N.The nonsubsampled contourlet transform:theory，design，and applications［J］.IEEE Transactions on Image Processing，2006，15(10):3089-3101.

［13］Zheng S，Shi W Z，Liu J，et al.Multi-source image fusion method using support value transform ［J］.IEEE Transactions on Image Processing，2007，16(7):1831-1839.

［14］Qu G，Zhang D L，Yan P F.Information measure for performance of image fusion ［J］.Electronic Letters，2002，38(7):313-315.

［15］Li S T，Yang B，Hu J.Performance comparison of different multi-resolution transforms for image fusion［J］.Information Fusion，2011，12(2):74-84.