基于剪切波變換的醫(yī)學(xué)圖像融合算法

阿都建華 王邦平 王 珂 王 艷

(成都信息工程學(xué)院軟件工程學(xué)院，成都610225)

引言

醫(yī)學(xué)影像技術(shù)已經(jīng)是臨床診斷中不可替代的技術(shù)手段，在臨床診斷中得到了廣泛地應(yīng)用，并取得了巨大的成效。由于不同的成像設(shè)備采用的傳感器原理不同，醫(yī)學(xué)影像技術(shù)有多種成像模式，分別產(chǎn)生不同的醫(yī)學(xué)圖像，這些圖像反映人體臟器信息也有所不同。醫(yī)學(xué)圖像主要包括形態(tài)圖像、結(jié)構(gòu)圖像和功能圖像。醫(yī)學(xué)顯微圖像屬于形態(tài)圖像，X線圖像、超聲、CT、MRI 等屬于結(jié)構(gòu)圖 像，PET、SPECT、FMRT 及EIT 等屬于功能圖像［1］。在臨床上，將不同模態(tài)的醫(yī)學(xué)圖像進(jìn)行適當(dāng)融合，使多種模態(tài)的信息有機(jī)結(jié)合，可獲得信息更為豐富的復(fù)合圖像，為臨床診斷和治療提供更為可靠、準(zhǔn)確的依據(jù)［2-8］。

近年來，由于基于多尺度幾何分析的圖像融合方法具有多分辨率、稀疏描述等特性，被廣泛應(yīng)用于圖像處理領(lǐng)域。小波變換是最典型的多尺度分解工具，然而小波變換不能很好地表示二維圖像中的線奇異性，容易導(dǎo)致偽吉布斯現(xiàn)象。剪切波是一種克服了小波變換缺點(diǎn)的新穎多尺度幾何分析工具，由Guo 及其合作者于2007 年通過特殊形式的具有合成膨脹的仿射系統(tǒng)構(gòu)造的一種接近最優(yōu)的多維函數(shù)稀疏表示方法［9］。剪切波［9-11］具有簡單的數(shù)學(xué)結(jié)構(gòu)，可以通過對一個函數(shù)進(jìn)行伸縮、平移、旋轉(zhuǎn)而生成一組基函數(shù)，使其和多分辨分析關(guān)聯(lián)起來。與圖像融合中常用的其他多尺度分解工具相比，剪切波變換具有自身的優(yōu)勢。剪切波變換不僅具有與曲線波［12］相同的接近最優(yōu)的非線性誤差逼近階，而且在頻率空間中剪切波是逐層細(xì)分的，因而具有更好的表示性能。同時，剪切波在剪切過程中和輪廓波［13-14］不同，它沒有方向數(shù)目和支撐基尺寸大小的限制，而且剪切波逆變換只需對正向變換中的剪切濾波器進(jìn)行加和處理，不需要像輪廓波變換那樣對方向濾波器逆合成，因此剪切波的實現(xiàn)過程比輪廓波具有更高的計算效率。

Miao 等將剪切波變換應(yīng)用于圖像融合，并證明了它的融合結(jié)果比其他多尺度方法捕獲更多的圖像細(xì)節(jié)信息［15］。鄭紅等將剪切波變換應(yīng)用于紅外與可見光圖像融合，并針對紅外與可見光圖像的特性，采用區(qū)域顯著性方法融合分解系數(shù)［16］。筆者在深入研究剪切波變換的基礎(chǔ)上，結(jié)合醫(yī)學(xué)圖像的特性，提出了一種基于剪切波變換的醫(yī)學(xué)圖像融合算法。該算法首先通過剪切波變換，對原始圖像進(jìn)行多尺度、多方向的分解，得到原始圖像的高頻子帶與低頻子帶，然后對低頻子帶的系數(shù)采用非負(fù)矩陣分解［17］(NMF)融合算法進(jìn)行融合，對高頻子帶系數(shù)則設(shè)計了一種基于人類視覺系統(tǒng)特性的視覺能量對比度方法進(jìn)行選擇，最后通過重構(gòu)融合后的高、低頻子帶系數(shù)獲得最終的融合圖像。實驗證明，本研究提出的融合算法所獲取的融合圖像能融合不同模態(tài)醫(yī)學(xué)圖像的主要信息，提高了圖像的臨床診斷和治療價值，取得了良好的效果。

1 基于剪切波變換的圖像融合方法

1.1 剪切波變換

在論文［10］中，剪切波可定義為

令ψ ∈L2(R2)，需滿足以下條件:

1)對于ξ = (ξ1，ξ2)∈(ξ2≠0)，有(ε)=(ε1，ε2)=(ε1(ε2/ε1)，其中(ε)為ψ(ε)的傅里葉變換;

2)ψ1為連續(xù)小波，∈C∞(R)，suppψ1?［-2，- 1/2］∪［1/2，2］;

假定

對每一個j ≥0，存在

圖1 剪切波頻域剖分圖以及剖分子區(qū)域的幾何特征。(a)剪切波的頻域支撐;(b)剪切波頻域剖分圖Fig． 1 The tiling of the frequency by the shearlets and the geometric characteristics of the split region． (a)The tiling of the frequency by the shearlets;(b)The size of the frequency support of a shearlet

集合:

是L2(D0)∨= {f ∈L2(R2):supp ^f ?D0}的一個Parseval 框架。

1.2 基于Shearlet 變換的融合規(guī)則

首先對經(jīng)過嚴(yán)格配準(zhǔn)的原始圖像進(jìn)行剪切波變換，將原圖分解成由剪切波系數(shù)表示的低頻子帶Sc=0(i，j)與一系列高頻子帶(i，j)，其中c = 1，2，…，N，表示圖像的第c 個高頻子帶，k = 1，2，…，M，(i，j)為分解方向參數(shù)，表示剪切波系數(shù)所在方向的位置，從而實現(xiàn)圖像的多尺度，多方向分解;然后將低頻子帶系數(shù)和高頻子帶系數(shù)融合，再進(jìn)行逆運(yùn)算，得到最后的融合圖像。剪切波的具體分解步驟如下［16，18］:

1)通過非降采樣塔式分解，將原始圖像分解為各個尺度圖像T0，T1，…，TN，其中T0表示分解后的粗尺度圖像，T1～TN表示分解后的第1 細(xì)分尺度圖像至第N 細(xì)分尺度圖像;

2)通過快速二維傅里葉變換，將各細(xì)分尺度圖像T1，…，TN變換至頻域FT1，…，F(xiàn)TN;

3)將各細(xì)分尺度頻域圖像FT1，…，F(xiàn)TN作為剪切波濾波器組的輸入，通過的多方向分解得到各高頻子帶頻域系數(shù)，通過對的快速二維傅里葉反變換得到一系列高頻子帶系數(shù)(i，j);

4)基于步驟2 的結(jié)果T0，直接使S0(i，j)= T0，綜合步驟4 的結(jié)果(i，j)，即得到剪切波變換的全部系數(shù);

5)將分解后源圖像對應(yīng)的低頻子帶系數(shù)和高頻子帶系數(shù)按特定的融合算法進(jìn)行融合處理，分別得到融合圖像的低頻子帶系數(shù)和高頻子帶系數(shù);

6)通過剪切波逆變換，得到融合圖像。

圖2 融合算法的流程Fig．2 The flow diagram of fusion algorithm

2 基于剪切波變換的醫(yī)學(xué)圖像融合算法

2.1 融合算法流程

首先對經(jīng)過嚴(yán)格配準(zhǔn)的原始圖像進(jìn)行剪切波變換，將原圖分解成由剪切波系數(shù)表示的低頻子帶Sl0(x，y)與一系列高頻方向子帶Si，l(x，y)，其中i= 1，2，…，N 為圖像的第i 個高頻方向子帶，l = 1，2，…，M 為分解方向參數(shù)，(x，y)表示剪切波系數(shù)所在方向的位置，從而實現(xiàn)了原始圖像的多尺度、多方向分解。本研究采用NMF 算法完成低頻子帶系數(shù)的融合，得到融合后的低頻子帶系數(shù);對于高頻方向子帶系數(shù)融合則采用視覺能量對比度方法完成，得到融合后的高頻方向子帶系數(shù);最后通過Shearlet 逆運(yùn)算，得到融合圖像。本研究提出的融合算法流程如圖2 所示。

2.2 低頻子帶系數(shù)融合規(guī)則

2.2.1 非負(fù)矩陣分解理論

非負(fù)矩陣分解(non-negative matrix factorization，NMF)是一種新的矩陣分解方法［17］，由Lee 和Seung 于1999 年在Nature 上的一篇論文中提出。在矩陣分解的過程中，該方法始終約束所有的元素，將基作為非負(fù)數(shù)存在，即要求所有分量始終為純加性的描述，同時使分解后的所有分量也均為非負(fù)值，并且降低了矩陣的維度。非負(fù)矩陣分解(NMF)問題可描述為:已知一個非負(fù)矩陣V，求非負(fù)的n × r 矩陣W 和非負(fù)的r × m 矩陣H，使

對于 NMF 問題的求解，常用歐氏距離和Kullback-liebler 散度作為目標(biāo)函數(shù)，有

因為對于一個給定矩陣V，矩陣W 和H 的最佳選擇是要使得V 和WH 之間的重構(gòu)誤差最小。因此，NMF 的求解問題實際上是個優(yōu)化問題，可以描述為

上述公式可以應(yīng)用梯度下降法來求解得到結(jié)果，求解過程也是收斂的，Lee 和Seung 已經(jīng)在理論上對此做出了證明［19］，并得到求解矩陣W 和矩陣H的迭代規(guī)則，即

式中，“.* ”和“. /”分別表示矩陣中各元素的相乘和相除。

2.2.2 基于非負(fù)矩陣分解的融合算法

在傳感器成像的過程中，由于傳感器自身的因素或外界的各種影響，常常會引入各種噪聲，因此待融合的原始觀測圖像實際是就是客觀真實世界在成像過程中引入了這些噪聲而形成的。在非負(fù)矩陣分解算法中，可以假設(shè)V = WH + ε(ε 表示噪聲)，此時噪聲ε 會在式(12)所示的迭代算法中趨于收斂，這個過程恰好符合圖像融合的過程［20-21］。因此，聯(lián)系圖像融合過程，如果假設(shè)觀測圖像為V，真實圖像為W，噪聲為ε，那么V 可以理解為W 與ε 之和，這樣NMF 可以有效地應(yīng)用于圖像融合。

由NMF 算法理論可知，通過迭代運(yùn)算方法能夠針對原始數(shù)據(jù)矩陣V 得到一個基于部分的近似表示形式WH。其中，W 的列數(shù)(即特征基的數(shù)量r)是一個待定量，將直接決定特征子空間的維數(shù)。對于特定的數(shù)據(jù)集，隱藏在數(shù)據(jù)集內(nèi)部特征空間的維數(shù)是確定的。也就是說，當(dāng)選取的r 與實際數(shù)據(jù)集特征空間的維數(shù)一致時，所得到的特征空間以及特征空間的基最有意義。當(dāng)r =1 時，通過迭代算法將得到唯一一個含有源數(shù)據(jù)全部特征的特征基。

由上述內(nèi)容可知，NMF 與圖像融合能很好地結(jié)合在一起應(yīng)用。假設(shè)有k 幅來自于多傳感器的大小為m × n 的觀測圖像f1，f2，…，fk，那么可以將觀測圖像逐個元素地按照行優(yōu)先的方式存儲到一個列向量中，這樣就可以得到一個mn × k 的矩陣V。V中包含k 個列向量v1，v2，…，vk，每個列向量代表k幅觀測圖像中的一幅圖像的信息，有

對這個觀測矩陣V 進(jìn)行非負(fù)矩陣分解，分解時取r =1，則可得到一個唯一的特征基W。顯然，此時的W 包含了參與融合的k 幅圖像的完整特征，將特征基W 還原到源圖像的像素級上，即可得到比源圖像效果都好的圖像。

2.3 高頻子帶系數(shù)融合規(guī)則

2.3.1 人類視覺系統(tǒng)

圖像對比度是指對一幅圖像中灰度值反差大小的測度，反差越大意味著對比度越大，反差越小意味著對比度越小。人眼能夠分辨的亮度差異所要求的最小光亮度差值ΔL 稱為亮度辨別閾值，對比靈敏度是指人眼具有的分辨亮度差異的能力。亮度辨別閾值ΔL 也不是固定不變的，其大小在不同的背景亮度L 下并不相同。也就是說，即使客觀的亮度是增強(qiáng)的，但在亮度L 的變化沒有增加到L+ ΔL 之前，人眼是感受不到這些亮度變化的，只有當(dāng)亮度L 變化達(dá)到或超過L + ΔL 時，人眼才能感受到這些亮度變化。因此，在研究中，采用ΔL/L 這一相對靈敏度閾值才更為合理。實驗表明，在無背景光時，ΔL/L 的值近似為一個常數(shù)，其數(shù)值為0.02。

在對視覺系統(tǒng)的進(jìn)一步研究中發(fā)現(xiàn)，對比度敏感門限與背景亮度的關(guān)系更接近指數(shù)規(guī)律，同時得到了一些更準(zhǔn)確的對比度敏感函數(shù)(contrast sensitivity function)［22-25］。其中，根據(jù)文獻(xiàn)［24］，對比度敏感門限可表示為

式中:I0為當(dāng)I = 0 (背景亮度為零)時的對比度敏感門限;α 為冪函數(shù)的指數(shù)，是常數(shù)，可根據(jù)視覺生理實驗獲得，其值為0.6 ～0.7。

2.3.2 高頻子帶系數(shù)融合規(guī)則

人類視覺系統(tǒng)對圖像局部對比度的變化十分敏感，因為對比度反映了圖像內(nèi)紋理、邊緣等信息的變化特征，同時包含了圖像的高頻信息以及相對于局部圖像背景的強(qiáng)度。鑒于此，融合方法可有選擇地突出被融合圖像的對比度信息，以求達(dá)到良好的視覺效果。在高頻子帶系數(shù)的融合算法中，筆者利用人類視覺特性來提取圖像的細(xì)節(jié)信息，并根據(jù)高頻子帶系數(shù)的特點(diǎn)設(shè)計了一種視覺能量對比度方法，進(jìn)行融合圖像系數(shù)的選擇，有

式中:j 表示第j 個高頻子帶，l 表示第l 個方向;Rj，l(x，y)表示像素點(diǎn)(x，y)的灰度與局部背景平均灰度的局部對比度;Ej，l(x，y)是以像素點(diǎn)(x，y)為中心的局部區(qū)域能量和;ˉmj，l(x，y)為低頻子帶系數(shù)以(x，y)為中心的某一局部區(qū)域均值;α 為視覺常數(shù)，通常取0.6 ～0.7 間的值。

局部對比度Rj，l(x，y)表示高頻子帶圖像中的系數(shù)與通過上節(jié)方法融合后獲得的低頻子帶圖像的局部區(qū)域灰度平均值的比值，即局部對比度。由前述人類視覺系統(tǒng)對比度敏感度的內(nèi)容可知，其值的大小在一定程度上能反映出待融合高頻子帶圖像的像素點(diǎn)融合于圖像后的效果。Rj，l(x，y)的計算方式為

式中:Sj，l(x，y)為j 尺度、l 方向的 高頻子帶系數(shù);為以像素點(diǎn)(x，y)為中心、尺寸大小為m × n 的低頻子帶圖像的局部灰度均值，其中SL(x +i，y +j)為融合后的低頻子帶系數(shù)。

式中，Sj，l(x + i，y + j)為j 尺度、l 方向的高頻子帶系數(shù)。

局部區(qū)域能量和Ej，l(x，y)表示高頻子帶圖像的局部區(qū)域能量和，分別計算圖像f1和f2的局部區(qū)域能量值(m，n)和(m，n)，有

對于經(jīng)過嚴(yán)格配準(zhǔn)的源圖像f1和f2，則在經(jīng)過shearlet 分解后，先分別計算得到它們的視覺能量對比度，然后再選擇高頻子帶系數(shù)，有

在高頻子帶系數(shù)的選擇中，將選取視覺能量對比度較大者作為最后的融合圖像的低頻子帶系數(shù)。

3 結(jié)果

為了驗證上述圖像融合方法的有效性，選取兩組經(jīng)過嚴(yán)格配準(zhǔn)的醫(yī)學(xué)圖像進(jìn)行融合實驗。為了能充分比較本研究提出的融合算法，實驗1 和實驗2 中都采用了3 種圖像融合方法作為對比方法。，包括拉普拉斯金字塔(LP)圖像融合方法、基于小波變換(DWT)的圖像融合方法和基于簡單非下采樣Contourlet 變換(NSCT)的圖像融合方法。在實驗中，拉普拉斯金字塔分解的層次為4 層;基于小波變換的圖像融合方法和基于NSCT 的簡單融合方法都采用簡單的平均值方法選取低頻子帶系數(shù)，采用最大絕對值方法選取高頻子帶系數(shù)，NSCT 的多尺度分解層次為3，方向分解級數(shù)分別為4、3、2。對于所提出的算法，剪切波變換分解層數(shù)為3，2 個細(xì)分尺度分解方向數(shù)分別為6，區(qū)域大小尺寸取3 × 3 。在實驗中，還采用了信息熵(information entropy，IE)、空間頻率(spatial frequency，SF)、標(biāo)準(zhǔn)差(standard deviation，SD)、平均梯度(average grads，AG)和均方根交叉熵(root mean square cross entropy，RCE)作為融合圖像的客觀評價指標(biāo)。

實驗1 采用嚴(yán)格配準(zhǔn)后的MRI 圖像和CT 圖像作為實驗的源圖像，如圖3 所示。其中，圖3(a)是CT 圖像，圖3(b)是MRI 圖像，圖像的尺寸大小均為256 像素×256 像素，且已經(jīng)過嚴(yán)格配準(zhǔn)。圖3(c)～圖3(f)分別為基于DWT 方法獲取的融合圖像、基于LP 方法獲取的融合圖像、基于NSCT 簡單方法獲取的融合圖像和基于本研究提出方法獲取的融合圖像。從主觀視覺角度可以明顯看出，圖3(c)的邊緣位置明顯出現(xiàn)了偽吉布斯現(xiàn)象，圖3(e)則非常暗淡，MRI 圖像信息明顯融入較少。圖3(e)相對圖3(f)亮度較暗，可視性效果略差。表1 是實驗1 的客觀評價指標(biāo)，可以看出本研究提出的方法在IE、SF、AG 三項指標(biāo)中均取得了最大值，在RCE指標(biāo)中取得了最小值。IE 最大，說明融合圖像包含的信息量最多，融合效果最好;SF 最大，說明圖像在空間域的變化最活躍，融合圖像的質(zhì)量最好;AG 最大，說明圖像中的細(xì)節(jié)最豐富，圖像清晰程度最高;RCE 最小，說明融合圖像和原始圖像之間的差異水平最小，即融合效果最好。SD 指標(biāo)的值小于DWT方法，由于SD 表示像素點(diǎn)的灰度值與融合圖像平均值的偏離程度，其值越大說明融合圖像中細(xì)節(jié)變化越豐富，融合質(zhì)量就越好，因此僅該項客觀指標(biāo)不是最優(yōu)值。由上述分析可知，圖3(f)在這4 項客觀指標(biāo)中取得了最好的融合結(jié)果，僅SD 一項指標(biāo)的計算結(jié)果不是最優(yōu)。因此，從客觀評價指標(biāo)同樣可以看出，本研究提出方法的獲取的融合圖像效果最好。

圖3 MRI 圖像與CT 圖像的融合結(jié)果比較。(a)CT 圖 像;(b)MRI 圖像;(c)DWT 融合圖像;(d)LP 融合圖像;(e)NSCT 融合圖像;(f)本方法融合圖像Fig． 3 Fusion results on MRI image and CT image． (a)CT image;(b)MRI image;(c)Fusion result of the DWT-based method;(d)Fusion result of the LP-based method;(e)Fusion result of the NSCT-based method;(f)Fusion result of the proposed method

表1 實驗1 融合圖像客觀評價指標(biāo)比較Tab． 1 The objective evaluation on the experiment 1

實驗2 采用Brain Web［26］中提供的不同形態(tài)大腦圖像作為實驗圖像，融合結(jié)果如圖4 所示。其中，圖4(a)是T1 圖像，圖4(b)是PD 圖像，圖像的尺寸大小均為216 像素×181 像素，且已經(jīng)過嚴(yán)格配準(zhǔn)。圖4(c)～圖4(f)分別為基于DWT 方法獲取的融合圖像、基于LP 方法獲取的融合圖像、基于NSCT簡單方法獲取的融合圖像和基于本研究提出方法獲取的融合圖像。從主觀視覺角度可以明顯看出，圖4(c)～圖4(e)的圖像效果明顯較差，圖像中邊緣信息比較模糊，整體對比度差。圖4(f)的融合質(zhì)量相對較好，融合信息豐富，整體對比度高，顯得更為清晰。因此，在實驗2 中所提出的融合方法取得了更好的主觀視覺融合效果。表2 是實驗2 的客觀評價指標(biāo)，從中可以看出，本方法同樣在IE、SD、SF

及RCE 指標(biāo)中取得了最優(yōu)值，僅AG 一項指標(biāo)的計算結(jié)果不是最優(yōu)。由實驗1 中對各項客觀評價指標(biāo)的分析可知，圖4(f)在大部分客觀指標(biāo)中取得了最好的融合結(jié)果。因此，從客觀評價指標(biāo)同樣可以看出，本方法獲取的融合圖像效果最好。

圖4 T1 圖像與PD 圖像的融合結(jié)果比較。(a)T1 圖像;(b)PD 圖像;(c)DWT 融合圖像;(d)LP 融合圖像;(e)NSCT 融合圖像;(f)本方法融合圖像Fig．4 Fusion results on T1 image and PD image．(a)T1 image;(b)PD image;(c)Fusion result of the DWT-based method;(d)Fusion result of the LPbased method;(e)Fusion result of the NSCT-based method;(f)Fusion result of the proposed method

表2 實驗2 融合圖像客觀評價指標(biāo)比較Tab． 2 the objective evaluation on the experiment 2

4 討論和結(jié)論

由于圖像傳感器種類繁多，應(yīng)用環(huán)境和目的各不相同，所以圖像融合算法也具有一定的針對性。筆者主要針對不同模態(tài)醫(yī)學(xué)圖像進(jìn)行融合研究，使多種模態(tài)的信息能更好地有機(jī)結(jié)合，形成一幅新的信息豐富、易于判別的融合圖像，以提高臨床診斷和治療的準(zhǔn)確性。多尺度幾何分析工具在圖像融合中已經(jīng)廣泛應(yīng)用，其有效性也已經(jīng)在實踐中得到證明。不同的多尺度幾何分析工具具有不同的特點(diǎn)，相應(yīng)的圖像融合效果也各不相同。本研究采用最新的多尺度幾何分析工具——剪切波變換作為圖像分解變換的數(shù)學(xué)工具，有效克服了小波變換的缺陷。此外，在分解后得到的高頻和低頻子帶圖像融合中，采用的融合規(guī)則對融合圖像的質(zhì)量至關(guān)重要，也是基于多尺度幾何分析圖像融合中的核心研究內(nèi)容。

本研究以一種新穎的多尺度幾何分析工具——剪切波變換為基礎(chǔ)，提出一種新的醫(yī)學(xué)圖像融合算法。本算法利用剪切波變換的多分辨率和多方向性特點(diǎn)，將待融合的原始圖像進(jìn)行分解，得到低頻子帶圖像和高頻方帶圖像。在低頻子帶系數(shù)的融合中，采用非負(fù)矩陣分解方法，再利用人類視覺系統(tǒng)對局部對比度較為敏感的特點(diǎn)，提出最大視覺能量對比度方法，進(jìn)行高頻方向子帶系數(shù)的融合，最后通過剪切波逆變換得到融合圖像。實驗證明，本研究提出的融合方法比傳統(tǒng)的小波方法、拉普拉斯金字塔方法、NSCT 方法，在主觀融合質(zhì)量和客觀評價指標(biāo)方面都有較大提高。

［1］ 陳曉艷，李健楠，王化祥. 一種電阻抗圖像與CT 圖像融合方法研究［J］.中國生物醫(yī)學(xué)工程學(xué)報，2012，30(6):892 -896.

［2］ Constantinos S，Pattichis MS，Micheli-Tzanakou E. Medical imaging fusion applications:An overview［C］ //Conference Record of the Thirty-Fifth Asilomar Conference on Signals，Systems and Computers． Pacific Grove:IEEE，2001:1263 -1267.

［3］ Das S，Kundu MK. NSCT-based multimodal medical image fusion using pulse-coupled neural network and modified spatial frequency ［J］. Medical and Biological Engineering and Computing，2012，50(10):1105 -1114.

［4］ Li Shutao， Yin Haitao， Fang Leyuan. Group-sparse representation with dictionary learning for medical image denoising and fusion［J］. IEEE Transactions on Bio-medical Engineering，2012，59(12):3450 - 3459.

［5］ Singh R，Srivastava R，Prakash O，et al． Multimodal medical image fusion in dual tree complex wavelet transform domain using maximum and average fusion rules［J］. Journal of Medical Imaging and Health Informatics，2012，2(2):168 -173.

［6］ Townsend DW，Beyer T. A combined PET/CT scanner:the path to true image fusion［J］. British journal of Radiology，2002，75(9):S24 -S30.

［7］ Yang Liu，Guo Baolong，Ni Wei. Multimodality medical image fusion based on multiscale geometric analysis of contourlet transform［J］.Neurocomputing，2008，72(1):203 -211.

［8］ 楊立才，劉延梅，劉欣，等. 基于小波包變換的醫(yī)學(xué)圖像融合方法［J］. 中國生物醫(yī)學(xué)工程學(xué)報，2009，28(1):12 -16.

［9］ Guo Kanghui，Labate D. Optimally sparse multidimensional representation using shearlets［J］. SIAM journal on mathematical analysis，2008，39(1):298 -318.

［10］ Kutyniok G，Labate D. Resolution of the wavefront set using continuous shearlets［J］. Trans Amer Math Soc，2009，361(5):2719 -2754.

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

［12］ Candes EJ，Donoho DL. Curvelets:A surprisingly effective nonadaptive representation for objects with edges［C］//in Curve and Surface Fitting:Saint-Malo，Nashville，USA:Vanderbilt University Press，2000:105 -120.

［13］ Do MN，Vetterli M. The contourlet transform:an efficient directional multiresolution image representation ［J］. IEEE Transactions on Image Processing，2005 ，14 (12):2091 -2106.

［14］ Cunha AL，Zhou J，Do MN. The nonsubsampled contourlet transform: theory， design， and applications ［J］. IEEE Transactions on Image Processing，2006 ，15 (10):3089 -3101.

［15］ Miao Qiguang，Shi Cheng，Xua Pengfei，et al． A novel algorithm of imagefusion using shearlets［J］. Optics Communications，2011，284(6):1540 -1547.

［16］ 鄭紅，鄭晨，閆秀生，陳海霞. 基于剪切波變換的可見光與紅外圖像融合算法［J］. 儀器儀表學(xué)報，2012，33(7):1613 -1619.

［17］ Lee DD，Seung HS. Learning the parts of objects by non-negative matrix factorization［J］. Nature，1999，401(6755):788 -791.

［18］ Lim WQ. The discrete shearlet transform:A new directional transform and compactly supported shearlet frames［J］. IEEE Transactions on Image Processing，2010 ，19(5):1166 -1180.

［19］ Seung D，Lee L. Algorithms for non-negative matrix factorization［J］. Advances in Neural Information Processing Systems，2001，13(2):556 -562.

［20］ 苗啟廣，王寶樹. 基于非負(fù)矩陣分解的多聚焦圖像融合研究［J］. 光學(xué)學(xué)報，2005，25(6):755 -759.

［21］ 苗啟廣，王寶樹. 圖像融合的非負(fù)矩陣分解算法［J］. 計算機(jī)輔助設(shè)計與圖形學(xué)學(xué)報，2005，17(7):2029 -2032.

［22］ Legge GE，F(xiàn)oley JM. Contrast masking in human vision［J］.JOSA，1980，70(12):1458 -1471.

［23］ Watson AB. DCT quantization matrices visually optimized for individual images［C］// Human Vision，Visual Processing and Digital Display IV. San Jose:SPIE，1993:202 -216.

［24］ Watson AB. Efficiency of a model human image code［J］. JOSA A，1987，4(12):2401 -2417.

［25］ Legge GE. A power law for contrast discrimination［J］. Vision Research，1981，21(4):457 -467.

［26］ http://www.bic.mni.mcgill.ca/brainweb/.