越来越多的医疗设备能采集人体不同特征数据并形成三维图像。在临床应用中, 通常需要将多幅含有不同重要信息的源图像融合到一幅图像中以辅助治疗。而传统的图像融合方法主要针对二维图像, 这些方法直接应用于三维数据会导致第三维信息的损失。本文结合近期出现的超小波变换——三维带限剪切波变换和四组传统融合准则, 提出新的三维磁共振图像融合方法。最后通过四组人脑T2*和磁量图(QSM)源数据, 将本文的方法分别与基于二维和三维的小波、双树复数小波的融合方法进行对比。实验显示基于三维变换的融合方法性能普遍优于基于二维变换的方法; 三维方法中, 剪切波克服了传统小波变换缺乏表达方向的缺陷, 能有效地提高传统图像融合算法的性能, 因此, 本文的方法质量指标最高。
引用本文： 段昶, 汪学刚, 王洪, 王帅. 基于三维带限剪切波变换的磁共振图像融合. 生物医学工程学杂志, 2015, 32(1): 181-186, 196. doi: 10.7507/1001-5515.20150033 复制
引言
医学图像融合是图像融合的一种，已经研究了数十年，许多方法也被广泛地应用于临床诊断中[1-2]。融合是指将不同设备[如X线计算机断层摄影(computer tomography，CT)、磁共振成像(magnetic resonance imaging，MRI)等]采集的重要信息提取并合并成为一幅图像。不同设备或同类设备不同配置所生成的图像中包含的重要信息是不同的，有些信息具有相似性，但大部分信息是互补的。CT图像能提供人体稠密、坚硬组织的信息，而MRI图像则主要提供软组织的信息，如T2*能提供组织弛豫时间的对比信息，磁量图(quantitative susceptibility mapping，QSM)可提供由多种内部磁性生物标记物或铁、钙、钆等对比剂引起的磁敏感对比信息[3]。通常而言，图像融合需要先将源图像配准，而少数情况如T2*和QSM图像产生于同一次扫描数据，即二者已经完全配准了。
当前医学图像融合研究主要考虑的是二维图像的处理，但是现在许多医疗设备都能直接采集和生成三维图像。三维图像中每个点的灰度值不仅与同层邻近点互相关，也与相邻层中的邻近点互相关。传统的二维融合方法直接应用于三维数据会导致第三维相关信息的丢失，因此有必要研究能直接处理三维图像的融合方法。
融合算法可以在空域或变换域进行处理。空域中，融合图像通常是源数据的加权平均，这样的方法简单易于实现，但融合图像质量不高，边缘等信息都有损失。变换域方法遵循以下步骤：源图像变换到变换域，按融合准则对图像系数进行处理，得到融合后的系数，最后将系数变回到空域，输出即为融合图像。这类算法中研究重点主要集中在两点：变换的选取，融合准则的设计。许多多尺度变换都能应用于融合算法中，如离散小波变换(discrete wavelet transform，DWT)、双树复数小波变换(dual-tree complex wavelet transform，DTCWT)、曲波(curvelet)[4]、剪切波(shearlet)[5]等等。
剪切波变换是近几年提出并逐步成熟的对多维数据高效表示的多尺度变换。实际上，针对小波变换缺乏对边缘等方向性特征稀疏表示的缺点，已提出和研究了许多其他的多尺度变换。但剪切波变换是所有方法中唯一同时拥有以下优点的工具：只有一个或有限个产生函数集合，能几乎最优地表示高维数据，对连续和离散数据采用相同方式处理，拥有紧支实现等。剪切波变换已广泛应用于图像处理中，如去噪[5]、缘检测[6]、增强[7]等。剪切波也同样适用于图像融合。
本文提出基于三维带限剪切波(band limited shearlet transform，BLST)的医学图像融合方法。融合准则方面，本文采用了最大点模值(max point’s modulus，MPM)、最大区域能量(maximum regional energy，MRE)[8]、最大绝对差值(maximum absolute difference，MAD)和最大改进拉普拉斯能量(maximum sum of modified Laplacian energy，MSMLE)[4] 4个准则。这4个准则在二维情况下是比较高效的融合准则，本文将其扩展到三维。
1 剪切波基础理论
三维剪切波是二维剪切波的扩展，能几乎最有效地表示三维数据中面状(二维)和线状(一维)奇异性特征。设：
| $ \begin{align} & {{A}_{{{2}^{j}}}}=\left(\begin{matrix} {{2}^{j}} & 0 & 0 \\ 0 & {{2}^{j/2}} & 0 \\ 0 & 0 & {{2}^{j}} \\ \end{matrix} \right), {{{\tilde{A}}}_{{{2}^{j}}}}=\left(\begin{matrix} {{2}^{j/2}} & 0 & 0 \\ 0 & {{2}^{j}} & 0 \\ 0 & 0 & {{2}^{j}} \\ \end{matrix} \right), \\ & {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{A}}}_{{{2}^{j}}}}=\left(\begin{matrix} {{2}^{j}} & 0 & 0 \\ 0 & {{2}^{j}} & 0 \\ 0 & 0 & {{2}^{j/2}} \\ \end{matrix} \right), {{S}_{R}}=\left(\begin{matrix} 1 & {{k}_{1}} & {{k}_{2}} \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{matrix} \right), \\ & {{{\tilde{S}}}_{k}}=\left(\begin{matrix} 1 & 0 & 0 \\ {{k}_{1}} & 1 & {{k}_{2}} \\ 0 & 0 & 1 \\ \end{matrix} \right), {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{S}}}_{k}}=\left(\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ {{k}_{1}} & {{k}_{2}} & 1 \\ \end{matrix} \right)\\ \end{align} $ |
其中j∈ Z。Mc=diag(c1, c2, c2), =diag(c2, c1, c2), c=diag(c2, c2, c1), c1＞0, c2＞0。三维剪切波可被定义为
| $ \begin{matrix} SH\left(\Phi, \psi, \tilde{\psi }, \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }, c \right)=\Phi \left(\Phi; {{c}_{1}} \right)\cup \Psi \left(\psi; c \right)\cup \tilde{\Psi }\left(\tilde{\psi }; c \right)\cup \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\Psi }\left(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }; c \right), 214\Phi \left(\Phi; {{c}_{1}} \right)=\left\{ {{\Phi }_{m}}=\Phi \left(\cdot-m \right):m\in {{c}_{1}}{\boldsymbol{{\text{Z}}}^{3}} \right\}, \\ \Psi \left(\psi; c \right)=\left\{ {{\psi }_{j, k, m}}={{2}^{j}}\psi \left({{S}_{k}}{{A}_{{{2}^{j}}}}\cdot-m \right)\ge 0, \left| k \right|\le \left\lceil {{2}^{j/2}} \right\rceil, m\in {{M}_{c}}{\boldsymbol{{\text{Z}}}^{3}} \right\}, \\ \tilde{\Psi }\left(\tilde{\psi }; c \right)=\left\{ {{{\tilde{\psi }}}_{j, k, m}}={{2}^{j}}\tilde{\psi }\left({{{\tilde{S}}}_{k}}{{{\tilde{A}}}_{{{2}^{j}}}}\cdot-m \right)\ge 0, \left| k \right|\le \left\lceil {{2}^{j/2}} \right\rceil, m\in {{{\tilde{M}}}_{c}}{{\text{Z}}^{3}} \right\} \\ \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\Psi }\left(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }; c \right)=\left\{ {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }}}_{j, k, m}}={{2}^{j}}\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }\left({{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{S}}}_{k}}{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{A}}}_{{{2}^{j}}}}\cdot-m \right)\ge 0, \left| k \right|\le \left\lceil {{2}^{j/2}} \right\rceil, m\in {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{M}}}_{c}}{\boldsymbol{{\text{Z}}}^{3}} \right\}, j\in {{\text{N}}_{0}}, k\in {\boldsymbol{{\text{Z}}}^{2}} \\ \end{matrix}$ |
对于频域实现的带限剪切波，其生成函数(ξ1)= 1(ξ1)2(ξ2/ξ1) 2(ξ3/ξ1)，其中Ψ1是一个带限小波，而Ψ2是尺度函数，表示Ψ的傅里叶变换。其频域划分情况由图 1表示，关于带限剪切波的更多知识可以参考文献[9-10]。
图1
金字塔三维剪切波频域划分
Figure1.
Partition of frequency domain for 3D pyramid adapted shearlet transform
2 融合方法
首先介绍融合准则。本文算法采用的融合准则如下：融合后的低频系数为源数据低频系数的均值，而高频部分采用4个准则，分别是MPM、MRE[8]、MAD和MSMLE准则[4]，其计算公式为
| $ {{C}_{f}}=\left\{ \begin{align} & {{C}_{a}}, {{S}_{a}}\ge {{S}_{b}} \\ & {{C}_{b}}, {{S}_{a}} < {{S}_{b}} \\ \end{align} \right., $ |
其中MPM准则中St=abs(Ct)，t∈{a, b}；MRE准则中St=σt=，t∈{a, b}，Ω为一个区域，
融合算法流程图如图 2所示，两幅源图像分别用Va和Vb表示，经过三维带限剪切波正变换后得到变换域系数Ca和Cb，通过融合准则，得到融合系数Cf，最后再通过三维带限剪切波逆变换，得到最终的融合图像Vf。
图2
本文提出的三维图像融合算法的流程图
Figure2.
Steps of proposed 3D image fusion method
3 实验分析
本节通过4组人脑磁共振T2*和QSM图像数据对本文提出算法进行性能分析和比较。所有人脑的研究数据经过了学院审查委员会认可。MR扫描使用GE Signa HDxt 3.0T系统，采用8通道头线圈。扫描序列使用3D T2*加权的多回波梯度回波序列，扫描参数为FA=20°；TR=57 ms；TE数=8；第一个TE=5.7 ms；TE间隔(ΔTE)=6.7 ms；BW=±41.67 kHz；FOV=24 cm；扫描分辨率0.57 mm×0.75 mm×2 mm。三维T2*幅度和QSM使用非线性形状限定偶极子反演方法来重建，并且插值到128×128×128。因为在QSM生成过程中，脑组织之外的磁场被噪声破坏，信息没有意义，所以QSM的区域使用一个掩膜(Mask)来裁剪。该掩膜使用了一个大脑裁剪工具[11]来生成。
由于源数据为三维数据，传统的二维算法只在x、y两个方向上分别对每层数据进行融合，而三维融合方法则能直接对整个数据源进行融合。两类算法的主要差别在于z轴信息的保持。对于z轴保持可通过帧间差分互信息量(inter frame difference mutual information，IFD_MI)来衡量[12-13]。本文采用的数据由于掩膜的存在，不能像文献[12-13]那样直接通过图像中所有点进行计算，计算指标时只使用掩膜内的点：设上述文献中互信息量用MIi(Dai, Dbi, Dfi)表示，Dai, Dbi, Dfi表示源图像和融合图像沿z轴方向的差分图像，Dti=V ti+1 -Vti，t∈{a, b, f}，i表示序号。则本文的帧间差分互信息量用式(2)进行计算，即
| $ \text{IDF }\!\!\_\!\!\text{ MI}=\frac{1}{N-1}\sum\limits_{i=1}^{N-1}{\left(M{{I}_{i}}\left(D_{a}^{i}, D_{b}^{i}, D_{f}^{i} \right)\cdot \left(Mas{{k}^{i+1}}\cup Mas{{k}^{i}} \right)\right), } $ |
其中N为数据第三维层数量(128)，·表示点乘，Mask表示掩膜，i表示序号。
从图 3可以看出，二维变换的帧间差分图像存在许多与源数据不一致的区域，不同变换和准则得到的结果也有很大区别。而三维变换的结果则与源图像非常相似，不同变换和准则的结果也几乎没有区别。这样的差异也可以从表 1看出，三维变换的IFD_MI值明显高于二维变换。而基于三维剪切波加MPM准则的方法拥有最高的IFD_MI值。由此可以得出结论，基于三维变换的方法相对基于二维变换的方法而言，具有更好的第三维信息保持能力。
图3
二维和三维算法的帧间差分图像
(a)源图像帧间差分图；(b)图例；(c)2D DWT；(d) 2D DTCWT；(e) 3D BLST；(f) 3D DWT；(g) 3D DTCWT
Figure3. Inter frame difference image for 2D and 3D methods(a) inter frame difference image for source; (b) legend; (c) 2D DWT; (d) 2D DTCWT; (e) 3D BLST; (f) 3D DWT; (g) 3D DTCWT
表1
4组数据的IFD_MI
Table1.
IFD_MI for the four subjects
图 4从左到右显示了冠状面、轴状面、矢状面各一层图像的源图像和不同方法的处理结果图像。从图中可以看出融合图像中都保持了两幅源图像中的重要特征：相对T2*图像，融合图像的亮度更一致，而纹理信息更多来自QSM。不同方法的融合图像的视觉效果几乎一致，单纯通过目测很难判定算法的优劣，还需采用性能评价指标。本文采用了互信息量(mutual information，MI)和QAB｜F[14]来评价。在文献[14]中，QAB｜F只给出二维的情况，本文将其扩展到三维。同样，所有的值均是基于掩膜内的数据计算得到。表 2列出4组数据的MI和QAB｜F。通过观察可以发现，除了第二组一个指标外，指标最高的方法是基于三维带限剪切波变换和MRE融合准则的方法。
图4
轴状面、冠状面、矢状面源图像和不同方法的融合图像
(a)轴状面、冠状面、矢状面源图像(左：T2*图；右：QSM图)；(b)3D BLST波变换融合结果(左上：MPM；右上：MRE；左下：MAD；右下：MSMLE)；(c)其他变换与最大区域能量准则的融合结果(左上：2D DWT；右上：2D DTCWT；左下：3D DWT；右下：3D DTCWT)
Figure4. Sagittal, coronal, axial source images and fused slices of different methods(a) sagittal, coronal, axial source images (Left: T2* images; Right: QSM images); (b) the fused images based on 3D BLST (Up-left: MPM; Up-right: MRE; Bottom-left: MAD; Bottom-right: MSMLE); (c) the fused images based on other transform with MRE rule (Up-left: 2D DWT; Up-right: 2D DTCWT; Bottom-left: 3D DWT; Bottom-right: 3D DTCWT)
表2
4组数据的质量指标
Table2.
Performance indices for four subjects
4 结论
现代医疗设备能提供人体结构三维图像，而传统的二维图像融合算法如果直接应用于三维图像，会导致第三维信息的损失。本文提出了基于带限剪切波变换的三维图像融合方法而言，并通过4组人脑MRI图像对不同的融合方法性能进行了评估。实验结果表明：三维图像融合方法相对二维图像融合方法而言，在第三维信息的保持能力上更具优势；相对三维小波变换、三维双树复数小波变换而言，本文提出的融合方法具有更高的性能指标且同样适用于其他设备(如CT、PET等)采集并已配准的医学图像。
引言
医学图像融合是图像融合的一种，已经研究了数十年，许多方法也被广泛地应用于临床诊断中[1-2]。融合是指将不同设备[如X线计算机断层摄影(computer tomography，CT)、磁共振成像(magnetic resonance imaging，MRI)等]采集的重要信息提取并合并成为一幅图像。不同设备或同类设备不同配置所生成的图像中包含的重要信息是不同的，有些信息具有相似性，但大部分信息是互补的。CT图像能提供人体稠密、坚硬组织的信息，而MRI图像则主要提供软组织的信息，如T2*能提供组织弛豫时间的对比信息，磁量图(quantitative susceptibility mapping，QSM)可提供由多种内部磁性生物标记物或铁、钙、钆等对比剂引起的磁敏感对比信息[3]。通常而言，图像融合需要先将源图像配准，而少数情况如T2*和QSM图像产生于同一次扫描数据，即二者已经完全配准了。
当前医学图像融合研究主要考虑的是二维图像的处理，但是现在许多医疗设备都能直接采集和生成三维图像。三维图像中每个点的灰度值不仅与同层邻近点互相关，也与相邻层中的邻近点互相关。传统的二维融合方法直接应用于三维数据会导致第三维相关信息的丢失，因此有必要研究能直接处理三维图像的融合方法。
融合算法可以在空域或变换域进行处理。空域中，融合图像通常是源数据的加权平均，这样的方法简单易于实现，但融合图像质量不高，边缘等信息都有损失。变换域方法遵循以下步骤：源图像变换到变换域，按融合准则对图像系数进行处理，得到融合后的系数，最后将系数变回到空域，输出即为融合图像。这类算法中研究重点主要集中在两点：变换的选取，融合准则的设计。许多多尺度变换都能应用于融合算法中，如离散小波变换(discrete wavelet transform，DWT)、双树复数小波变换(dual-tree complex wavelet transform，DTCWT)、曲波(curvelet)[4]、剪切波(shearlet)[5]等等。
剪切波变换是近几年提出并逐步成熟的对多维数据高效表示的多尺度变换。实际上，针对小波变换缺乏对边缘等方向性特征稀疏表示的缺点，已提出和研究了许多其他的多尺度变换。但剪切波变换是所有方法中唯一同时拥有以下优点的工具：只有一个或有限个产生函数集合，能几乎最优地表示高维数据，对连续和离散数据采用相同方式处理，拥有紧支实现等。剪切波变换已广泛应用于图像处理中，如去噪[5]、缘检测[6]、增强[7]等。剪切波也同样适用于图像融合。
本文提出基于三维带限剪切波(band limited shearlet transform，BLST)的医学图像融合方法。融合准则方面，本文采用了最大点模值(max point’s modulus，MPM)、最大区域能量(maximum regional energy，MRE)[8]、最大绝对差值(maximum absolute difference，MAD)和最大改进拉普拉斯能量(maximum sum of modified Laplacian energy，MSMLE)[4] 4个准则。这4个准则在二维情况下是比较高效的融合准则，本文将其扩展到三维。
1 剪切波基础理论
三维剪切波是二维剪切波的扩展，能几乎最有效地表示三维数据中面状(二维)和线状(一维)奇异性特征。设：
| $ \begin{align} & {{A}_{{{2}^{j}}}}=\left(\begin{matrix} {{2}^{j}} & 0 & 0 \\ 0 & {{2}^{j/2}} & 0 \\ 0 & 0 & {{2}^{j}} \\ \end{matrix} \right), {{{\tilde{A}}}_{{{2}^{j}}}}=\left(\begin{matrix} {{2}^{j/2}} & 0 & 0 \\ 0 & {{2}^{j}} & 0 \\ 0 & 0 & {{2}^{j}} \\ \end{matrix} \right), \\ & {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{A}}}_{{{2}^{j}}}}=\left(\begin{matrix} {{2}^{j}} & 0 & 0 \\ 0 & {{2}^{j}} & 0 \\ 0 & 0 & {{2}^{j/2}} \\ \end{matrix} \right), {{S}_{R}}=\left(\begin{matrix} 1 & {{k}_{1}} & {{k}_{2}} \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{matrix} \right), \\ & {{{\tilde{S}}}_{k}}=\left(\begin{matrix} 1 & 0 & 0 \\ {{k}_{1}} & 1 & {{k}_{2}} \\ 0 & 0 & 1 \\ \end{matrix} \right), {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{S}}}_{k}}=\left(\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ {{k}_{1}} & {{k}_{2}} & 1 \\ \end{matrix} \right)\\ \end{align} $ |
其中j∈ Z。Mc=diag(c1, c2, c2), =diag(c2, c1, c2), c=diag(c2, c2, c1), c1＞0, c2＞0。三维剪切波可被定义为
| $ \begin{matrix} SH\left(\Phi, \psi, \tilde{\psi }, \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }, c \right)=\Phi \left(\Phi; {{c}_{1}} \right)\cup \Psi \left(\psi; c \right)\cup \tilde{\Psi }\left(\tilde{\psi }; c \right)\cup \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\Psi }\left(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }; c \right), 214\Phi \left(\Phi; {{c}_{1}} \right)=\left\{ {{\Phi }_{m}}=\Phi \left(\cdot-m \right):m\in {{c}_{1}}{\boldsymbol{{\text{Z}}}^{3}} \right\}, \\ \Psi \left(\psi; c \right)=\left\{ {{\psi }_{j, k, m}}={{2}^{j}}\psi \left({{S}_{k}}{{A}_{{{2}^{j}}}}\cdot-m \right)\ge 0, \left| k \right|\le \left\lceil {{2}^{j/2}} \right\rceil, m\in {{M}_{c}}{\boldsymbol{{\text{Z}}}^{3}} \right\}, \\ \tilde{\Psi }\left(\tilde{\psi }; c \right)=\left\{ {{{\tilde{\psi }}}_{j, k, m}}={{2}^{j}}\tilde{\psi }\left({{{\tilde{S}}}_{k}}{{{\tilde{A}}}_{{{2}^{j}}}}\cdot-m \right)\ge 0, \left| k \right|\le \left\lceil {{2}^{j/2}} \right\rceil, m\in {{{\tilde{M}}}_{c}}{{\text{Z}}^{3}} \right\} \\ \overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\Psi }\left(\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }; c \right)=\left\{ {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }}}_{j, k, m}}={{2}^{j}}\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{\psi }\left({{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{S}}}_{k}}{{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{A}}}_{{{2}^{j}}}}\cdot-m \right)\ge 0, \left| k \right|\le \left\lceil {{2}^{j/2}} \right\rceil, m\in {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{M}}}_{c}}{\boldsymbol{{\text{Z}}}^{3}} \right\}, j\in {{\text{N}}_{0}}, k\in {\boldsymbol{{\text{Z}}}^{2}} \\ \end{matrix}$ |
对于频域实现的带限剪切波，其生成函数(ξ1)= 1(ξ1)2(ξ2/ξ1) 2(ξ3/ξ1)，其中Ψ1是一个带限小波，而Ψ2是尺度函数，表示Ψ的傅里叶变换。其频域划分情况由图 1表示，关于带限剪切波的更多知识可以参考文献[9-10]。
图1
金字塔三维剪切波频域划分
Figure1.
Partition of frequency domain for 3D pyramid adapted shearlet transform
2 融合方法
首先介绍融合准则。本文算法采用的融合准则如下：融合后的低频系数为源数据低频系数的均值，而高频部分采用4个准则，分别是MPM、MRE[8]、MAD和MSMLE准则[4]，其计算公式为
| $ {{C}_{f}}=\left\{ \begin{align} & {{C}_{a}}, {{S}_{a}}\ge {{S}_{b}} \\ & {{C}_{b}}, {{S}_{a}} < {{S}_{b}} \\ \end{align} \right., $ |
其中MPM准则中St=abs(Ct)，t∈{a, b}；MRE准则中St=σt=，t∈{a, b}，Ω为一个区域，
融合算法流程图如图 2所示，两幅源图像分别用Va和Vb表示，经过三维带限剪切波正变换后得到变换域系数Ca和Cb，通过融合准则，得到融合系数Cf，最后再通过三维带限剪切波逆变换，得到最终的融合图像Vf。
图2
本文提出的三维图像融合算法的流程图
Figure2.
Steps of proposed 3D image fusion method
3 实验分析
本节通过4组人脑磁共振T2*和QSM图像数据对本文提出算法进行性能分析和比较。所有人脑的研究数据经过了学院审查委员会认可。MR扫描使用GE Signa HDxt 3.0T系统，采用8通道头线圈。扫描序列使用3D T2*加权的多回波梯度回波序列，扫描参数为FA=20°；TR=57 ms；TE数=8；第一个TE=5.7 ms；TE间隔(ΔTE)=6.7 ms；BW=±41.67 kHz；FOV=24 cm；扫描分辨率0.57 mm×0.75 mm×2 mm。三维T2*幅度和QSM使用非线性形状限定偶极子反演方法来重建，并且插值到128×128×128。因为在QSM生成过程中，脑组织之外的磁场被噪声破坏，信息没有意义，所以QSM的区域使用一个掩膜(Mask)来裁剪。该掩膜使用了一个大脑裁剪工具[11]来生成。
由于源数据为三维数据，传统的二维算法只在x、y两个方向上分别对每层数据进行融合，而三维融合方法则能直接对整个数据源进行融合。两类算法的主要差别在于z轴信息的保持。对于z轴保持可通过帧间差分互信息量(inter frame difference mutual information，IFD_MI)来衡量[12-13]。本文采用的数据由于掩膜的存在，不能像文献[12-13]那样直接通过图像中所有点进行计算，计算指标时只使用掩膜内的点：设上述文献中互信息量用MIi(Dai, Dbi, Dfi)表示，Dai, Dbi, Dfi表示源图像和融合图像沿z轴方向的差分图像，Dti=V ti+1 -Vti，t∈{a, b, f}，i表示序号。则本文的帧间差分互信息量用式(2)进行计算，即
| $ \text{IDF }\!\!\_\!\!\text{ MI}=\frac{1}{N-1}\sum\limits_{i=1}^{N-1}{\left(M{{I}_{i}}\left(D_{a}^{i}, D_{b}^{i}, D_{f}^{i} \right)\cdot \left(Mas{{k}^{i+1}}\cup Mas{{k}^{i}} \right)\right), } $ |
其中N为数据第三维层数量(128)，·表示点乘，Mask表示掩膜，i表示序号。
从图 3可以看出，二维变换的帧间差分图像存在许多与源数据不一致的区域，不同变换和准则得到的结果也有很大区别。而三维变换的结果则与源图像非常相似，不同变换和准则的结果也几乎没有区别。这样的差异也可以从表 1看出，三维变换的IFD_MI值明显高于二维变换。而基于三维剪切波加MPM准则的方法拥有最高的IFD_MI值。由此可以得出结论，基于三维变换的方法相对基于二维变换的方法而言，具有更好的第三维信息保持能力。
图3
二维和三维算法的帧间差分图像
(a)源图像帧间差分图；(b)图例；(c)2D DWT；(d) 2D DTCWT；(e) 3D BLST；(f) 3D DWT；(g) 3D DTCWT
Figure3. Inter frame difference image for 2D and 3D methods(a) inter frame difference image for source; (b) legend; (c) 2D DWT; (d) 2D DTCWT; (e) 3D BLST; (f) 3D DWT; (g) 3D DTCWT
表1
4组数据的IFD_MI
Table1.
IFD_MI for the four subjects
图 4从左到右显示了冠状面、轴状面、矢状面各一层图像的源图像和不同方法的处理结果图像。从图中可以看出融合图像中都保持了两幅源图像中的重要特征：相对T2*图像，融合图像的亮度更一致，而纹理信息更多来自QSM。不同方法的融合图像的视觉效果几乎一致，单纯通过目测很难判定算法的优劣，还需采用性能评价指标。本文采用了互信息量(mutual information，MI)和QAB｜F[14]来评价。在文献[14]中，QAB｜F只给出二维的情况，本文将其扩展到三维。同样，所有的值均是基于掩膜内的数据计算得到。表 2列出4组数据的MI和QAB｜F。通过观察可以发现，除了第二组一个指标外，指标最高的方法是基于三维带限剪切波变换和MRE融合准则的方法。
图4
轴状面、冠状面、矢状面源图像和不同方法的融合图像
(a)轴状面、冠状面、矢状面源图像(左：T2*图；右：QSM图)；(b)3D BLST波变换融合结果(左上：MPM；右上：MRE；左下：MAD；右下：MSMLE)；(c)其他变换与最大区域能量准则的融合结果(左上：2D DWT；右上：2D DTCWT；左下：3D DWT；右下：3D DTCWT)
Figure4. Sagittal, coronal, axial source images and fused slices of different methods(a) sagittal, coronal, axial source images (Left: T2* images; Right: QSM images); (b) the fused images based on 3D BLST (Up-left: MPM; Up-right: MRE; Bottom-left: MAD; Bottom-right: MSMLE); (c) the fused images based on other transform with MRE rule (Up-left: 2D DWT; Up-right: 2D DTCWT; Bottom-left: 3D DWT; Bottom-right: 3D DTCWT)
表2
4组数据的质量指标
Table2.
Performance indices for four subjects
4 结论
现代医疗设备能提供人体结构三维图像，而传统的二维图像融合算法如果直接应用于三维图像，会导致第三维信息的损失。本文提出了基于带限剪切波变换的三维图像融合方法而言，并通过4组人脑MRI图像对不同的融合方法性能进行了评估。实验结果表明：三维图像融合方法相对二维图像融合方法而言，在第三维信息的保持能力上更具优势；相对三维小波变换、三维双树复数小波变换而言，本文提出的融合方法具有更高的性能指标且同样适用于其他设备(如CT、PET等)采集并已配准的医学图像。
