基于小波變換的多模態(tài)醫(yī)學(xué)圖像的融合方法

詹玲超　閔芳　郭慧敏

摘 要： 針對(duì)單一模態(tài)的圖像不能提供足夠的信息，可以將互補(bǔ)的圖像進(jìn)行融合。提出一種基于小波變換的醫(yī)學(xué)圖像融合方法。首先對(duì)待融合的醫(yī)學(xué)圖像進(jìn)行小波變換，分別得到低頻和高頻分量。對(duì)于低頻分量，采用基于局部區(qū)域能量的方法進(jìn)行融合；對(duì)于高頻分量，采用基于局部區(qū)域系數(shù)之差的加權(quán)和方法進(jìn)行融合；最后通過(guò)小波逆變換得到融合圖像。實(shí)驗(yàn)結(jié)果表明，相比同類研究方法，該算法在保證圖像質(zhì)量的同時(shí)可增強(qiáng)圖像的空間細(xì)節(jié)表現(xiàn)能力。

關(guān)鍵詞： 圖像融合； 醫(yī)學(xué)圖像； DWT； 區(qū)域能量； 相鄰系數(shù)

中圖分類號(hào)： TN911.73?34； TP391 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 1004?373X（2016）24?0096?04

Multimodal medical image fusion method based on discrete wavelet transform

ZHAN Lingchao， MIN Fang， GUO Huimin

（Department of Information Engineering， Nanhang Jincheng College， Nanjing 211156， China）

Abstract： Since the single modal image can′t provide enough information， it is necessary to fuse the complementary images. A medical image fusion method based on discrete wavelet transform （DWT） is proposed. DWT is performed for the multimodal medical images under fusion to obtain the low?frequency and high?frequency components respectively. A method based on the local regional energy is adopted to fuse low?frequency components. The weighted sum method based on difference of the local regions coefficients is used to fuse the high?frequency components. The fusion image is obtained by means of the inverse DWT. The experimental results prove that， in comparison with the similar research methods， the method can enhance the spatial detail expressive ability of the image while ensuring the image quality.

Keywords： image fusion； medical image； DWT； regional energy； adjacent coefficient

0 引 言

隨著醫(yī)學(xué)影像工程和計(jì)算機(jī)技術(shù)的不斷發(fā)展，出現(xiàn)了很多成像設(shè)備，這樣就有多模態(tài)的醫(yī)學(xué)圖像，這些圖像對(duì)人體臟器和病變組織的分析有不同的作用，比如CT圖像有較強(qiáng)的空間分辨率和幾何特性，對(duì)骨骼成像非常清晰，但對(duì)軟組織的對(duì)比度則較低；MR圖像可清晰反映軟組織等解剖結(jié)構(gòu)，但對(duì)鈣化點(diǎn)不敏感，并且受到磁干擾會(huì)發(fā)生幾何失真。這樣可以看出對(duì)于人體同一解剖結(jié)構(gòu)所得到的不同圖像，在形態(tài)和功能信息上是互補(bǔ)的。單一模態(tài)的圖像往往不能提供足夠的信息，可以將互補(bǔ)的圖像進(jìn)行融合。

圖像融合技術(shù)[1?5]就是將來(lái)自多個(gè)傳感器的多幅源圖像融合成一幅新的圖像，這幅融合圖像具有更多的信息和更高的可信度。圖像融合技術(shù)應(yīng)用在很多領(lǐng)域，醫(yī)學(xué)圖像融合就是其中的一個(gè)。圖像融合可以分為像素級(jí)融合、特征級(jí)融合和決策級(jí)融合。目前的融合方法大多為像素級(jí)融合，常用的有PCA[6]、金字塔變換[7]、小波變換[8?11]等融合方法。其中PCA是一種選取最優(yōu)像素權(quán)值的方法，缺點(diǎn)是無(wú)法突出光譜特性，不適用于相關(guān)性弱的圖像融合。金字塔變換的缺點(diǎn)是層間具有相關(guān)性，導(dǎo)致融合結(jié)果不夠理想。

小波變換對(duì)保留圖像信息具有相當(dāng)好的性能，通過(guò)不同的融合規(guī)則可以得到不同的結(jié)果。本文提出了一種基于小波變換的醫(yī)學(xué)圖像的融合方法。先將待融合的兩幅圖像進(jìn)行小波變換，分別得到一個(gè)低頻和若干個(gè)高頻分量，然后分別對(duì)低頻和高頻分量進(jìn)行融合。對(duì)于低頻分量，采用局部區(qū)域能量取大的方法進(jìn)行融合；對(duì)于高頻分量，采用局部區(qū)域系數(shù)之差的加權(quán)和的方法進(jìn)行融合；最后通過(guò)小波逆變換得到融合結(jié)果。

本文算法的融合結(jié)果與幾種算法的融合結(jié)果作對(duì)比，采用熵、平均梯度和互信息三種客觀評(píng)價(jià)標(biāo)準(zhǔn)進(jìn)行評(píng)價(jià)，結(jié)果表明本文算法有不錯(cuò)的效果，融合結(jié)果更加清晰。

1 IVIFDWT算法

本文提出IVIFDWT融合算法，即一種基于小波變換多模態(tài)醫(yī)學(xué)圖像的融合方法（Multimodal Medical Image Fusion Based on Discrete Wavelet Transform）。

1.1 IVIFDWT算法規(guī)則的設(shè)計(jì)

小波變換的層數(shù)越多，融合結(jié)果的細(xì)節(jié)信息就越豐富，但隨著分解層數(shù)的增加，圖像的結(jié)構(gòu)信息量損失變大，這些信息在小波逆變換后都是不能恢復(fù)的，所以小波分解的層數(shù)不能過(guò)高，一般在2～5層之間。通過(guò)實(shí)驗(yàn)證明且考慮到運(yùn)行時(shí)間，本文里選2層小波分解。

定義1：設(shè)待融合的兩幅圖像分別為A和B，將這兩幅圖像分別進(jìn)行兩層小波變換。兩層小波變換后，它們的低頻分量分別表示為[LA（x，y）]和[LB（x，y）]，高頻分量分別表示為[Hh1A（x，y）]，[Hv1A（x，y）]，[Hd1A（x，y）]，[Hh2A（x，y）]，[Hv2A（x，y）]，[Hd2A（x，y）]，[Hh1B（x，y）]，[Hv1B（x，y）]，[Hd1B（x，y）]，[Hh2B（x，y）]，[Hv2B（x，y）]，[Hd2B（x，y）]，分別為兩層小波變換后A圖像第一層和第二層水平、垂直、對(duì)角高頻分量和B圖像第一層和第二層水平、垂直、對(duì)角高頻分量。

這樣對(duì)圖像A和B的低頻和高頻分量，分別以不同的融合規(guī)則進(jìn)行融合。本文算法規(guī)則如下：

（1） 低頻分量融合規(guī)則是將兩幅圖像低頻分量按照系數(shù)的區(qū)域能量來(lái)決定。需要計(jì)算低頻分量的區(qū)域能量，之后根據(jù)區(qū)域能量的大小來(lái)得到低頻系數(shù)的融合結(jié)果。

（2） 高頻分量融合規(guī)則是先計(jì)算高頻系數(shù)局部區(qū)域系數(shù)之差，再求出這些差的加權(quán)和，如果圖像A的加權(quán)和大于圖像B的加權(quán)和，選取圖像A的高頻系數(shù)，否則選取圖像B的高頻系數(shù)。

1.2 IVIFDWT算法的原理圖

已知待融合的兩幅圖像分別為圖像A和圖像B，IVIFDWT算法的原理圖如圖1所示。

1.3 IVIFDWT算法的描述

本文設(shè)計(jì)的IVIFDWT算法先對(duì)待融合的兩幅圖像分別進(jìn)行小波變換，然后使用不同的融合規(guī)則分別對(duì)低頻系數(shù)和高頻系數(shù)進(jìn)行融合。具體算法描述如下：

1.3.1 低頻系數(shù)融合方法

將[A]，[B]兩幅圖像的低頻系數(shù)的分別劃分成[N×N]的區(qū)域窗口（嘗試將圖像劃分為[3×3]和[5×5]區(qū)域窗口分別計(jì)算，實(shí)驗(yàn)證明[3×3]的效果相對(duì)較好，以下都劃分為[3×3]大小的區(qū)域窗口），然后計(jì)算每個(gè)窗口的能量。假設(shè)用[EA]和[EB]分別表示圖像[A]和圖像[B]低通子帶系數(shù)每個(gè)窗口的能量，能量公式為[8]：[EA（x，y）=i，j=-N2N2（LA（x+i，y+j））2log（（LA（x+i，y+j））2）] （1）

[EB（x，y）=i，j=-N2N2（LB（x+i，y+j））2log（（LB（x+i，y+j））2）] （2）

據(jù)以上的能量公式就可以得到一個(gè)初步的融合圖（Initial Fusion Decision Map），其公式為：

[iL_map（x，y）=1，EA（x，y）≥EB（x，y）0，EA（x，y）

融合規(guī)則是，如果圖像[A]的區(qū)域窗口能量大于圖像[B]的區(qū)域窗口能量，則對(duì)應(yīng)的融合圖里的值就是1，否則就是0。

最終的融合圖[fL_map（x，y）]（Final Fusion Decision Map）是通過(guò)濾波運(yùn)算的一致性驗(yàn)證得到的。對(duì)于之前得到的初步融合圖，如果在一個(gè)[3×3]大小的窗口區(qū)域內(nèi)，中間的是0，其他大部分是0，這時(shí)中間的值不變；如果中間的為0，其他大部分是1，這時(shí)將中間的值改為1。

可以根據(jù)最終的融合圖來(lái)得到低頻系數(shù)的融合結(jié)果，其結(jié)果為：

[L（x，y）=LA（x，y）， fL_map（x，y）=1LB（x，y）， fL_map（x，y）=0] （4）

如果在最終的融合圖中值為1，表示這個(gè)位置取圖像[A]的低頻系數(shù)，如果值為0，則取圖像[B]的低頻系數(shù)。

1.3.2 高頻系數(shù)融合方法

高頻系數(shù)反映圖像的邊緣和細(xì)節(jié)信息，這里充分考慮了相鄰系數(shù)，采用相鄰系數(shù)之差加權(quán)和的融合方法。首先計(jì)算每個(gè)系數(shù)和其相鄰系數(shù)之差，以兩幅圖像的第一層水平系數(shù)為例，其公式為[4]：

[Mh1A（x，y）=2Hh1A（x，y）-Hh1A（x-1，y）-Hh1A（x+1，y）+ 2Hh1A（x，y）-Hh1A（x，y-1）-Hh1A（x，y+1）+ 2Hh1A（x，y）-Hh1A（x-1，y-1）-Hh1A（x+1，y+1）+ 2Hh1A（x，y）-Hh1A（x-1，y+1）-Hh1A（x+1，y-1）] （5）

[Mh1B（x，y）=2Hh1B（x，y）-Hh1B（x-1，y）-Hh1B（x+1，y）+2Hh1B（x，y）-Hh1B（x，y-1）-Hh1B（x，y+1）+2Hh1B（x，y）-Hh1B（x-1，y-1）-Hh1B（x+1，y+1）+2Hh1B（x，y）-Hh1B（x-1，y+1）-Hh1B（x+1，y-1）] （6）

則其相鄰系數(shù)之差加權(quán)和的公式[4]為：

[WMh1A（x，y）=i，j=-N2N2ω（i+1，j+1）Mh1A（x+i，y+j）] （7）

[WMh1B（x，y）=i，j=-N2N2ω（i+1，j+1）Mh1B（x+i，y+j）] （8）

式中，[N=3]，[ω]為[4]：

[ω=116121242121] （9）

根據(jù)式（7）和式（8）計(jì)算出的局部區(qū)域系數(shù)之差加權(quán)和，得到一個(gè)初步的高頻系數(shù)融合圖，其公式為：

[ih1H_map（x，y）=1，WMh1A（x，y）≥WMh1B（x，y）0，WMh1A（x，y）

融合規(guī)則是，如果圖像[A]的局部區(qū)域系數(shù)之差加權(quán)和大于圖像[B]的局部區(qū)域系數(shù)之差加權(quán)和，則對(duì)應(yīng)的融合圖里的值就是1，否則就是0。

最終的融合圖[fh1H_map（x，y）]同樣是通過(guò)濾波運(yùn)算的一致性驗(yàn)證得到的。對(duì)于之前得到的初步融合圖，如果在一個(gè)[3×3]大小的窗口區(qū)域內(nèi)，中間的是0，其他大部分是0，這時(shí)中間的值不變；如果中間的為0，其他大部分是1，這時(shí)將中間的值改為1。

得到最終的融合圖之后，可以根據(jù)最終的融合圖來(lái)得到高頻系數(shù)的融合結(jié)果，其結(jié)果為：

[Hh1（x，y）=Hh1A（x，y），fh1H_map（x，y）=1Hh1B（x，y），fh1H_map（x，y）=0] （11）

如果在最終的融合圖中值為1，表示這個(gè)位置取圖像[A]的高頻系數(shù)，如果值為0，則取圖像[B]的高頻系數(shù)。高頻的其他分量依次類推。

對(duì)融合后的低頻系數(shù)和高頻系數(shù)進(jìn)行小波逆變換，得到最終的融合結(jié)果。

2 實(shí)驗(yàn)及結(jié)果分析

2.1 融合結(jié)果

本文算法在Matlab環(huán)境中進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)中的圖像來(lái)自http：//www.med.harvard.edu/aanlib/home.html。本文算法應(yīng)用于CT和MR的融合，MR?T1和MR?T2的融合。

2.1.1 CT和MR的融合

圖2列出了其中兩組CT和MR圖像，其中圖2（a），圖2（c）是CT圖像，圖2（b），圖2（d）是MR圖像。

<實(shí)驗(yàn)結(jié)果如圖3所示。實(shí)驗(yàn)里有三組對(duì)比實(shí)驗(yàn)分別是像素取平均值法（P_avg），如圖3（a），圖3（e）所示；小波變換低頻分量取平均，高頻分量取大值法（DWT_avg_max），如圖3（b），圖3（f）所示；NSCT變換低通子帶系數(shù)取平均，帶通子帶系數(shù)取大值法（NSCT_avg_max），如圖3（c），圖3（g）所示；本文算法如圖3（d），圖3（h）所示。

2.1.2 MR?T1和MR?T2的融合

圖4列出了其中兩組MR?T1和MR?T2圖像，其中圖4（a），圖4（c）是MR?T1 images，圖4（b），圖4（d）是MR?T2 images。

實(shí)驗(yàn)結(jié)果如圖5所示。對(duì)比實(shí)驗(yàn)同樣分別是像素取平均值法（P_avg），如圖5（a），圖5（e）所示；小波變換低頻分量取平均，高頻分量取大值法（DWT_avg_max），如圖5（b），圖5（f）所示；NSCT變換低通子帶系數(shù)取平均，帶通子帶系數(shù)取大值法（NSCT_avg_max），如圖5（c），圖5（g）所示；本文算法如圖5（d），圖5（h）所示。

2.2 結(jié)果分析

將實(shí)驗(yàn)結(jié)果用信息熵、平均梯度和互信息三種評(píng)價(jià)標(biāo)準(zhǔn)進(jìn)行評(píng)價(jià)。信息熵是衡量圖像信息豐富程度的一個(gè)重要指標(biāo)，熵值越大說(shuō)明信息量越大。平均梯度反映了圖像中的微小細(xì)節(jié)方差、紋理變化特征和圖像的清晰度，平均梯度越大表示圖像清晰度越好?；バ畔⑹菍?duì)兩幅圖像共享信息的測(cè)量，互信息的值越大，圖像的融合效果越好。表1是對(duì)CT和MR的融合結(jié)果的評(píng)價(jià)，表2為MR?T1和MR?T2的融合結(jié)果的評(píng)價(jià)。

從表1和表2中可以看出，將本文提出的融合方法和另外三種融合方法相比，本文方法效果比較好。表1是CT和MR融合結(jié)果的評(píng)價(jià)，除了第一組數(shù)據(jù)DWT_avg_max中的信息熵，本文方法都是最優(yōu)的。表2是MR?T1和MR?T2融合結(jié)果的評(píng)價(jià)，除了第二組數(shù)據(jù)DWT_avg_max中的信息熵，本文方法都是最優(yōu)的。

3 結(jié) 語(yǔ)

本文提出了一種基于小波變換的紅外和可見光圖像的融合方法。對(duì)于低頻系數(shù)，采用基于局部區(qū)域能量的方法進(jìn)行融合；對(duì)于高頻系數(shù)，采用基于局部區(qū)域系數(shù)之差加權(quán)和的方法進(jìn)行融合。結(jié)果表明不管是從視覺效果方面來(lái)看，還是從融合結(jié)果的客觀評(píng)價(jià)來(lái)看，相對(duì)同類研究方法，本文方法獲得的圖像都更好，更為清晰，達(dá)到較為滿意的效果。

參考文獻(xiàn)

[1] 段先華，操佳佳，劉佳.改進(jìn)的PCNN模型在多光譜與全色圖像融合中的應(yīng)用研究[J].現(xiàn)代電子技術(shù)，2014，37（3）：55?60.

[2] 許良鳳，林輝，胡敏.基于差分進(jìn)化算法的多模態(tài)醫(yī)學(xué)圖像融合[J].電子測(cè)量與儀器學(xué)報(bào)，2013，27（2）：110?114.

[3] 康晨，曾丹，沈潔，等.圖像融合中的彩色圖像顏色校正[J].電子測(cè)量技術(shù)，2014，37（3）：54?57.

[4] GANASALA P， KUMAR V. CT and MR image fusion scheme in nonsubsampled contourlet transform domain [J]. Journal of digit imaging， 2014， 27（3）： 407?418.

[5] BHATNAGAR G， WU Q M J， LIU Z. Directive contrast based multimodal medical image fusion in NSCT domain [J]. IEEE transactions on multimedia， 2013， 15（5）： 1014?1024.

[6] KUMAR M， DASS S. A total variation?based algorithm for pixel?level image fusion [J]. IEEE transactions on image processing， 2009， 18（9）： 2137?2143.

[7] BURT P J. A gradient pyramid basis for pattern selective image fusion [C]// Proceedings of 1992 Society for Information Display Conference. San Jose： SID Press， 1992： 467?470.

[8] 王健，鄭少鋒.基于YUV與小波變換的可見光與紅外圖像融合[J].西北工業(yè)大學(xué)學(xué)報(bào)，2013，33（3）：208?211.

[9] 徐小軍，王友仁，陳帥.基于下采樣分?jǐn)?shù)階小波變換的圖像融合新方法[J].儀器儀表學(xué)報(bào)，2014，35（9）：2061?2069.

[10] SINGH R， KHARE A. Fusion of multimodal medical images using Daubechies complex wavelet transform： a multi resolution approach [J]. Information fusion， 2012， 19（3）： 49?60.

[11] ROY S， HOWLADER T， RAHMAN S M M. Image fusion technique using multivariate statistical model for wavelet coefficients [J]. Signal， image and video processing， 2013， 7（2）： 355?365.