随着医学成像技术的不断发展, 医学影像已成为临床诊断和治疗中不可或缺的一部分. 而这就要求医学影像能够具备较高的空间分辨率, 从而对高密度组织, 软组织有高质量的成像效果, 这一要求仅仅依靠单模态医学图像无法实现. 例如, 计算机断层扫描技术(CT)能够对硬组织的病变位置进行准确定位, 有利于显现高密度的骨骼特征; 磁共振血管成像(MRI)可以对低密度的软组织, 血管等清晰成像; 单光子发射计算机断层成像(SPECT)主要用来显示人体器官, 肿瘤组织血流和代谢功能^[1]. 由于单一模态图像成像机制不同会造成的图像细节信息纹理缺失, 有必要将不同模态的图像进行信息融合, 以实现优势互补, 使图像信息更丰富, 形成全面清晰的医学影像.

融合后的图像应该尽可能的保留源图像的信息, 不造成细节部分的缺失, 同时也要展现纹理特征, 更加有效地诊断医学图像. 近来, 基于变换域的图像融合算法是研究热点, 常用的变换有小波变换^[2], 轮廓波变换, 非下采样轮廓波变换(NSCT)^[3], 剪切波变换等.这些变换或使图像的细节边缘部分表现能力不佳; 或是不具备平移不变性, 会引起图像失真, 无法得到清晰的边缘; 或是计算过于复杂, 需耗费大量时间, 使得实用性大大降低. 为此, NSST被提出, 作为剪切波变换的改进模型, 它突出了图像中的显著特性, 并且成功地应用于医学和非医学图像处理应用中. NSST具有灵敏度高, 方向性强, 运算速度快等特点, 优于小波, 轮廓波等其它变换, 具有处理多方向变化的能力, 是一种非常强大的变换^[1].

国内外学者对基于NSST的图像融合进行了相关研究, 并做出诸多贡献. 陈贞等人^[4]采用基于NSST的医学图像融合算法, 低频采用区域能量加权, 高频采用区域能量和平均梯度加权的融合方法, 取得较好的数据结果, 但在一定程度上忽略了低频子带的边缘信息; 赵丹等人^[5]提出区域特征值的融合方法, 虽然可以综合反映图像的区域特征, 但是没有考虑到高低频系数之间的影响; 吴方等人^[6]提出使用蜘蛛优化算法搜索函数的最优参数集合, 进一步增强融合图像的质量, 但该算法的计算效率比较低; Mohammed等人^[7]提出的一种在改进的NSCT中结合多模态图像的方法, 可以提供高质量的图像, 有助于临床诊断的有效性, 但是SR和PCNN计算比较费时; Padmavathi等人^[8]提出的基于DTCWT-PCA和DTCWT-PSO的多模态医学图像融合方法, DTCWT-PCA与DTCWT-PSO算法具有较好的性能, 特征向量尺寸较小, 计算效率较高, 但是实际融合效果不佳; 文献[9]提出利用空间域主成分分析(PCA)法从大数据集中提取相关信息, 从小波变换中提取图像的方向性细节, 此方法虽然可以增强图像的方向性特征, 但是在一定程度上存在冗余细节, 噪声状伪影, 失真等问题; 郭楚等人^[10]通过多个坐标系的转换关系将超声坐标系配准到CT图像坐标系里, 最终将实时超声图像统一到CT图像中, 并在软件中测量融合误差, 融合误差得到了减小, 但是实时超声图像与CT图像的融合精度还是不够高; 邓志华等人^[11]提出在融合低频成分时, 采用传统的绝对值取大的融合策略, 以保留源图像更多的亮度信息, 但是却没有考虑到图像整体的轮廓清晰程度; Yang等人^[12]提出一种基于NSST和压缩感知的融合算法, 视觉效果较佳, 但是会造成图像信息的部分缺失.

为了提高图像整体的区域强度, 突出图像的边缘信息, 提高图像的对比度, 本文提出一种新的基于NSST的融合算法. 对于低频子带, 合成局域平均能量与局域标准差进行子带选择; 对于高频子带, 利用改进的拉普拉斯能量和的规则进行融合. 主观视觉效果以及客观评价均表明, 融合后的图像包含更丰富的纹理, 轮廓清晰, 对比度较高, 视觉效果更佳.

1 NSST模型基本理论

剪切波是在合成小波理论的前提下, 通过经典仿射将几何与多分辨分析结合而成的系统理论. 当维数n=2时, 具有合成膨胀的仿射系统为^[13,14]:

$ {A_{AB}}(\Psi ) \!= \!\left\{ {{\Psi _{j, l, k}}(x)\! =\! |\det A{|^{\frac{j}{2}}}\Psi \left( {{B^l}{A^j}x - k} \right):j, l \in Z, k \in {Z^2}} \right\} $

(1)

式中, $\Psi \in {{{L}}^2}\left( {{{{R}}^2}} \right)$ , A和B均为二维可逆方针, 且 $|\det B| = 1$ . 若对于任意的 $f \in {{{L}}^2}\left( {{{{R}}^2}} \right)$ , ${A_{AB}}(\Psi )$ 都构成一个Parseval框架(也叫紧支撑框架), 即:

$\sum\nolimits_{j,l,k} {{{\left| { < f,{\Psi _{j,l,k}}} \right|}^2} = {{\left\| {{f}} \right\|}^2}} $

(2)

则称 ${A_{AB}}(\Psi )$ 中的元素为合成小波. 其中, 矩阵A^j与尺度变换相关联, 矩阵B^l与区域保持的几何变换相关联, 如旋转与剪切变换. 一般情况下, 这种结构不仅允许像小波一样构造基元素在各个尺度和位置上的Parseval框架, 还可以构造基元素在各个方向上的Parseval框架.

一般令式(1)中的 $A = {A_0} = \left[ {\begin{array}{*{20}{l}} 4&0 \\ 0&2 \end{array}} \right]$ 为各向异性的膨胀矩阵, $B = {B_0} = \left[ {\begin{array}{*{20}{l}} 1&1 \\ 0&1 \end{array}} \right]$ 为剪切矩阵^[13]. 则可以获得的表达式为:

$\left\{ {{\varphi _{a, s, t}}\left( x \right) = {a^{ - \frac{3}{4}}}\phi \left( {A_a^{ - 1}B_s^{ - 1}x - t} \right), a \in {{{R}}^{\rm{ + }}}, s \in {{R}}, t \in {{{R}}^2}} \right\}$

(3)

其中, ${\varphi _{a, s, t}}\left( x \right)$ 为剪切波, 式(3)为剪切波系统. 可以看出, 剪切波是一个涵盖位置, 尺寸和方向的函数集合.

NSST的离散化过程主要分为两步: 多尺度分解和方向局部化. 其中, 多尺度分解是通过非下采样金字塔滤波器组来实现的, 图像经K级非下采样金字塔滤波器分解成1个低频和K个高频子带图像, 而且每个子带图像均与原图像大小相同; 方向局部化是通过剪切滤波器实现的^[13]. 标准的剪切波变换中使用的剪切滤波器是在伪极化坐标系中通过窗函数的平移操作实现的, 此过程中包含下采样操作, 因而不具有平移不变性. 而NSST将标准的剪切滤波器从伪极化网络系统映射回到笛卡尔左边系统, 从而有效地摒弃了下采样操作实现其平移不变性^[15]. 如图1为非下采样剪切波三层分解示意图, 图2为NSST梯形支撑区间.

2 多模态医学图像融合新算法

本文提出的医学图像融合算法流程图如图3所示. 首先, 将从标准库取出的已配准源图像A, B进行NSST分解, 分别得到图像的低频子带系数 $C_L^A(x, y)$ 和 $C_L^B(x, y)$ , 高频子带系数 $D_A^{l, k}(x, y)$ 和 $D_B^{l, k}(x, y)$ (这里, k为方向级数, l为分解层数); 接着, 结合低频子带系数和各层高频方向子带系数的特点, 通过合成局域平均能量与局域标准差的方法进行低频子带的选择, 并利用改进的拉普拉斯能量和的策略进行高频子带的融合; 最后, 对获得的图像融合系数进行逆NSST变换, 得到融合后的图像.

2.1 低频子带融合策略

低频子带反映图像的轮廓与基本信息. 融合后的图像要想要最大程度地保留源图像的基础信息, 重点在于恰当地处理低频系数. 常用的低频区域融合方法有: 加权平均法, 绝对值取大法, 标准差选择法等^[16]. 加权平均法容易造成图像整体亮度的下降, 视觉效果上, 对比度略显不足, 导致图像边缘纹理和内部细节信息的模糊; 绝对值取大法的优势在于, 可以较大程度地保留源图像的边缘纹理细节信息, 但是低亮度区域的纹理细节弱化较多; 标准差选择法可以优先选择边缘细节丰富的区域信息, 但是对噪声现象抵御力较弱^[5,14]. 为了解决这一问题, 本文采取局域平均能量与局域标准差相结合的方法进行低频部分的融合. 具体步骤如下:

1) 计算低频系数的均值 ${C_L}$ , 区域能量 ${E_L}$ 以及局域标准差 ${D_L}$ . 其中:

${C_L}(x, y) = \frac{1}{{M \times N}}\sum\limits_{(r = - (M - 1)/2)}^{(M - 1)/2} {\sum\limits_{c = - (N - 1)/2}^{(N - 1)/2} {{C_L}} } (x + r, y + c)$

(4)

${E_L}(x, y) = \frac{1}{{{{M}} \times {{N}}}}\sum\limits_{(r = - (M - 1)/2)}^{(M - 1)/2} {\sum\limits_{c = - (N - 1)/2}^{(N - 1)/2} {C_L^2} } (x + r, y + c){{w}}(x + r, y + c)$

(5)

$ {D_L}\left( {x, y} \right) = \sqrt {\frac{1}{{(M - 1)(N - 1)}}\sum\limits_{(r = - (M - 1)/2)}^{(M - 1)/2} {\sum\limits_{c = - (N - 1)/2}^{(N - 1)/2} {{{\left( {{C_L}(x + r, y + c) + {{\bar C}_L}(x, y)} \right)}^2}} } } $

(6)

窗口大小为M×N, 窗口函数为 $W = \left[ {\begin{array}{*{20}{l}} 1&2&1 \\ 2&3&2 \\ 1&2&1 \end{array}} \right]$ .

2) 根据区域能量 ${E_L}$ 及局域标准差 ${D_L}$ 构造图像区域特征 ${H_1}_I, {H_{\rm{2}}}_I\left( {I = A, B} \right)$ , 及 ${R_P}(P = A, B)$ . 其中,

${H_1}(x, y) = \frac{{{E_L}(x, y)}}{{\left( {1 + {D_L}(x, y)} \right)}}$

(7)

${H_2}(x, y) = {E_L}(x, y) + {D_L}(x, y)$

(8)

${R_P} = H_1^2 + H_2^2, \left( {P = A, B} \right)$

(9)

式(7)中的“1”用于防止 ${D_L}$ 为零时融合系数失效.R的值越大, 对应图像包含的信息就越丰富.

3) 设计低频系数融合方案如下:

$ {F_L}(x, y)\! =\! \left\{\!\!\!\! {\begin{array}{*{20}{c}} {C_L^A(x, y), \;\;{H_{1a}} > {H_{1b}}, {H_{2a}} > {H_{2b}}} \\ {C_L^B(x, y), \;\;{H_{1b}} > {H_{1a}}, {H_{2b}} > {H_{2a}}} \\ {{q_1} \times C_L^A(x, y) \!+ \!{q_2} \times C_L^B(x, y), \;{\rm{ otherwise }}} \end{array}} \right. $

(10)

其中, ${q_1} = {R_A}{{/}}{{{R}}_A} + {{{R}}_B} + 0.01$ , ${q_2} = 1 - {q_1}$ 为权值.

2.2 高频子带融合策略

高频子带图像含有边缘, 纹理等细节信息, 对准确地诊断病情十分重要^[17]. 改进的拉普拉斯能量和是图像清晰度的有效度量, 与方差, 梯度能量, 拉普拉斯能量等清晰度评价指标相比, 性能更为优越^[16,18], 融合后图像的边缘和细节保留更加丰富. 改进的拉普拉斯(Modified Laplacian, ML)定义为:

$\begin{aligned}[b] ML\left( {i, j} \right) =& |2f\left( {i, j} \right) - f\left( {i - 1, j} \right) - f\left( {i + 1, j} \right)| \\ & + |2f\left( {i, j} \right) - f\left( {i, j - 1} \right) - f\left( {i, j + 1} \right)| \\ &+ \frac{1}{{\sqrt 2 }}|2f\left( {i, j} \right) - f\left( {i - 1, j + 1} \right) - f\left( {i + 1, j - 1} \right)| \\ &+ \frac{1}{{\sqrt 2 }}|2f\left( {i, j} \right) - f\left( {i - 1, j - 1} \right) - f\left( {i + 1, j + 1} \right)| \end{aligned} $

(11)

这里将两个对角项的权重设为 $1/\sqrt 2 $ , 考虑到拉普拉斯对角系数的影响, 充分应用了相邻位置的信息^[19]. 若取大小为3×3的窗口, 那么区域内改进的拉普拉斯能量和定义为:

$NSM{L_{\left( {i, j} \right)}} = \sum\limits_{m = - 1}^{m = 1} {\sum\limits_{n = - 1}^{n = 1} {(ML(} } i + m, j + n){)^2}$

(12)

$f\left( {i, j} \right)$ 为带通子带图像在 $\left( {i, j} \right)$ 处的系数, “1”为步长.

最终, 设计高频系数融合方案如下:

${F_H} \!=\! \left\{ {\begin{aligned} & D_A^{l, k}(x, y),\left( {NSM{L_1} \!> \!NSM{L_2}} \right), \\ &\quad\quad\left| {NSM{L_1}\! -\! NSM{L_2}} \right| \!> \!0.5 \\ & D_B^{l, k}(x, y),\left( {NSM{L_1} < NSM{L_2}} \right),\\ & \quad\quad \left| {NSM{L_1} - NSM{L_2}} \right| > 0.5 \\ & {o_1} \times D_A^{l, k}(x, y) + {o_2} \times D_B^{l, k}(x, y), \\ & \quad\quad \left| {NSM{L_2}} \right| < 0.5 \end{aligned}} \right.$

(13)

其中, $NSM{L_1}$ 是源图像A的改进拉普拉斯能量和值, $NSM{L_2}$ 是源图像B的改进拉普拉斯能量和值, 另外 ${o_1} \!= \!NSM{L_1}/\left( {NSM{L_1} \!+\! NSM{L_2} \!+ \!0.01} \right)$ , ${o_2} \!=\! 1 \!-\! {o_1}$ 为权值.

3 实验结果与分析 3.1 灰度图像实验结果与分析

为了验证本文算法的有效性, 选用3种已有算法的数据作为实验对比数据, 并选择信息熵, 空间频率, 标准差和平均梯度对融合后的图像进行评价. 图4(c)和图5(c)是文献[4]提出的基于NSST的方法分解图像, 低频子带根据区域能量加权融合系数进行融合, 高频子带利用平均梯度和区域能量加权的融合策略: 图4(d)和图5(d)是文献[20]所提出的低频子带利用区域能量加权, 高频子带利用区域能量取大的融合方法^[20]; 图4(e)和图5(e)是文献[21]提出的低频区域采用信息熵加权, 高频子带图像梯度能量加权的融合方法^[21]; 图4(f)和图5(f)是应用本文算法所融合得到的融合图像. 表1和表2是应用不同融合算法产生的医学融合图像的质量评价数据.

为了实现本文算法的可比性, 所有实验均采用Window 10操作系统, 使用Matlab R2014a平台进行仿真. 实验选取3组经过配准的源图像, 大小均为256×256, 颜色深度为8 bit灰度^[5]. 如图所示, 图4(a)和4(b)是正常脑部CT/MRI图像, 图5(a)和图5(b)是多发性脑梗塞MR-T1/MR-T2图像.

从主观角度分析, 图4(d)和图4(e)的组织信息亮度稍显暗淡, 图4(d)中轮廓较模糊, 而图4(c)在轮廓和对比度上要优于4(d)和4(f), 但是纹理边缘不够清晰., 图5(c)细节结构不够清晰, 轮廓信息不完整, 对比度稍优于图5(d)和图5(e), 图5(d)对比度明显不够高, 边缘轮廓也不够清晰, 图5(e)对比度优于5(d), 但是细节信息有所缺失. 而从本文算法融合的图4(f)和图5(f)可以看出融合质量较高, 不仅很大程度上弥补了其他图像的不足, 而且细节信息全面, 轮廓细致, 纹理清晰, 视觉效果良好.

从客观角度分析, 观察表1和表2, 图4(d)和图5(d)以及图4(e)和图5(e)的各项评价指标数值普遍较低, 尤其是图5(d)的标准差和平均梯度指数过低, 说明对应算法得到的融合图像纹理特征不够清晰, 对比度不高; 图5(e)的信息熵稍稍高于前面两种, 但是标准差指数较低, 细节信息丢失过多. 相比其他3种算法, 本文算法在客观评价指数上大部分都优于其余几种算法, 使融合后的图像信息丰富, 轮廓完整, 更具实用性.

3.2 彩色图像实验结果与分析

灰度图像过后, 为了进一步验证本文算法的通用性, 将本文算法运用于彩色图像与灰度图像之间的融合. 由于彩色图像是三维的, 所以需要先将24位彩色图像从RGB空间转换到HIS空间, 然后利用NSST分解出I分量图像与灰度图像, 并用本文算法融合, 将融合后得到的灰度图像代替原来的I分量, 最后通过HIS逆变换得到最终的RGB彩色医学图像. 融合步骤如图6所示, 本文选取两组SPECT/MRI彩色图像进行仿真实验, 如图7(a)和图7(b)是转移性肺MRI/SPECT图像, 而图8(a)和图8(b)是正常脑组织MRI/SPECT图像.

从主观角度分析, 图7(c), 图8(c)与图7(e), 图8(e)整体轮廓不够清晰, 内部信息丢失过多, 纹理边缘不够理想, 图7(d), 图8(d)整体性效果稍稍好于图7(c), 图8(c)与图7(e), 图8(e), 但是该图像的对比度不佳, 降低了临床诊断的实用性. 而本文算法融合后的图像, 轮廓清晰, 纹理清楚, 信息全面.

从客观数据分析, 观察表3和表4, 本文算法融合的图像空间频率, 平均梯度等评价指数都是较高的, 表明图像像素活跃程度高, 图像清晰. 表3中(d), (e)的标准差较低, 说明图像细节有所缺失, 表4的(d)中信息熵与空间频率较低, 说明图像中包含的信息量较少, 空间活跃度较低. 本文算法得出的评价指标大部分都优于其他算法, 所以本文算法也适用于彩色图像融合, 且融合质量较高.

综上分析, 本文提出的基于NSST域的医学图像融合算法, 在有效保留源图像概貌的同时, 融合图像的灰度分布均匀, 细节信息丰富, 与源图像之间相似度高, 融合质量优, 效果好. 是一种效果良好, 具有一定有效性与通用性地医学影像融合方法.

4 结论与展望

医学图像融合要求融合后的图像尽可能保真, 准确丰富, 全面的体现源图像含有的信息. 基于此, 本文提出一种基于NSST的医学图像融合算法, 低频子带系数采用局域平均能量与局域标准差相结合的方法进行融合, 这样图像的基本信息可以较为完善的保留; 高频子带系数采用改进的拉普拉斯能量和的策略进行融合, 有利于提高图像抓取细节的能力. 实验结果表明, 本文融合算法得到的图像信息全面, 轮廓完整, 纹理清晰. 在人眼视觉体验和客观评价数据上都优于其他3种算法.但是也存在一定的局限性, 本文算法是针对医学图像展开研究的, 将此算法拓宽到其他领域如遥感, 红外与可见光等是下一阶段研究的重点.