肺部肿瘤跨模态图像融合的并行分解自适应融合模型

周涛，刘珊，董雅丽，白静，陆惠玲

1.北方民族大学计算机科学与工程学院，银川 750021；2.北方民族大学图像图形智能处理国家民委重点实验室，银川 750021；3.宁夏医科大学理学院，银川 750004

0 引 言

肺癌是造成全球死亡人数最多的癌症(Sung等，2021)。随着精准医疗技术的发展，正电子发射断层扫描(positron emission tomography，PET)和计算机断层扫描(computed tomography，CT)已成为检测肺部肿瘤的主要工具。CT影像有较高的空间分辨率，对骨组织显像清晰(Li等，2021)，但对病变部位显像效果较差(Khan等，2020)，尤其不能清楚地显示浸润性肿瘤信息(Zhang等，2020)。PET影像对软组织显像清晰，但对骨组织显像较差(Polinati和Dhuli，2020)。CT和PET影像融合技术能够综合病灶和组织器官的解剖学和功能信息，帮助医生对肺部肿瘤进行更准确的定位和定性诊断。因此，PET/CT图像融合技术在临床上具有重要意义。

周涛等人(2021)对经典像素级图像融合算法和框架进行综述。Zhang和Blum(1999)提出经典像素级图像融合框架。在此基础上，Piella(2003)对框架进行完善，将融合过程概括为图像分解和融合规则，使融合过程更加系统和规范。在图像分解方面，常用的有双树复小波变换(魏兴瑜，2015)和NSCT(non-subsampled contourlet transform)(da Cunha等，2006)等。NSCT将原图像分解为不同方向的高频子带，具有平移不变性和多方向性。Zhu等人(2019)在NSCT变换的基础上，使用相位一致高频融合规则，有效地增强图像细节特征。Li和Wu(2018)采用低秩表示(low-rank representation, LRR)方法，有效地提取原图像的全局和局部结构，但对噪声敏感。对此，Li和Wu(2022)提出潜在低秩表示(latent low-rank representation，LatLRR)方法，有效地将噪声部分与显著部分分离。

Piella(2003)将融合规则分为匹配测度、活性测度、决策模块和合成模块。在Piella框架的指导下，Polinati和Dhuli(2020)使用平均融合规则作为活性测度，该方法有效地平滑图像特征，但会降低图像对比度。Liu等人(2018)采用NSML(new sum of modified Laplacian)为活性测度，实验表明该方法比绝对值最大融合规则保留了更多的细节信息。由于图像信号是高度结构化的，人眼视觉系统能够自适应地提取图像背景中的结构信息，为了对图像结构失真程度的度量能够更符合人眼视觉效果，Wang等人(2004)提出结构相似性度量(structural similarity index measure，SSIM)，从亮度、对比度和结构3个方面度量图像相似性。刘卷舒和蒋伟(2018)采用基于SSIM建立的Piella指标作为适应度函数，通过蛙跳算法对一系列的超参数进行寻优，寻优得到的超参数参与融合操作，得到的融合图像效果较好。

因此，为了改善融合图像同时提取图像的细节信息和能量信息不理想问题，本文提出并行分解PET/CT医学图像融合方法，LatLRR具有提取显著特征的细节信息的能力，NSCT变换具有多方向细节信息的能力，所以本文对肺部肿瘤的跨模态医学图像使用NSCT变换分解来提取疾病的多方向细节信息，同时使用LatLRR变换提取原图像的显著特征的细节信息，最后将两个分支的信息进行融合，实现显著特征和多方向细节信息的互补。在融合规则上，使用模糊逻辑融合规则融合低频子带，能够有效解决图像中存在的模糊问题；采用基于Piella框架的高频子带融合规则，引入用SSIM作为活性测度，使融合图像更符合人类视觉系统的感知特性。

1 相关技术

1.1 NSCT变换

图1 NSCT分解示意图

1.2 LatLRR方法

Liu和Yan(2011)提出LatLRR方法，针对LRR得到的图像包含过多噪声数据，无法有效提取图像特征值的问题，LatLRR方法将输入的图像分解为低秩部分、显著部分和噪声部分。噪声部分在低秩的约束条件下被去除，其表达式为

X=ZX+LX+E

(1)

式中，X是输入图像，Z是低秩系数，L是显著性系数；ZX和LX分别表示低秩部分和显著部分；E表示稀疏噪声。LatLRR问题再次进行最优化求解，其最终表达式为

s.t.X=ZX+LX+E

(2)

2 并行分解自适应融合模型

Piella框架下的并行分解图像自适应融合模型的主要步骤如下：

1)[La,Ha]=NSCTdec(A)；[Lb,Hb]=NSCTdec(B)。图像A和B是已配准的CT和PET图像，NSCTdec()是NSCT变换。La和Lb分别是图像A和B经过NSCT变换后得到的低频子带，Ha={Hai(x,y)|i=1,2,…,8}和Hb={Hbi(x,y)|i=1,2,…,8}分别是图像A和B经过NSCT变换后得到的8个高频子带。

2)[Lrra,Sa,Ea]=LatLRR(A)；[Lrrb,Sb,Eb]=LatLRR(B)。LatLRR()是LatLRR分解。Lrra，Sa，Ea和Lrrb，Sb和Eb分别是图像A和B经过LatLRR分解后得到的低秩部分、显著部分和噪声部分。

3)FL=FLA(La,Lb)。FLA()函数是基于模糊逻辑自适应的低频子带融合规则(fuzzy logic adaptive fusion rules，FLA)，FL是低频子带融合图像。

4)FH=Piella(Ha,Hb)；Mab=Match(Ha,Hb)；αa=Active(Ha)；αb=Active(Hb)。

计算高频子带融合图像

FH=d×Ha+(1-d)×Hb

Piella()函数是基于Piella框架的高频子带融合规则，Match()函数是计算匹配测度的方法，Mab是两幅高频子带之间的相似值，Active()是计算活性测度的函数，αa和αb是两幅图像A和B经过Active函数计算得到的图像相异度，d是决策因子，FH是Piella函数的返回值，表示高频子带的融合图像。

算法融合框架如图2所示。

图2 并行分解图像自适应融合模型

2.1 低频融合规则

低频子带融合规则的设计从4个方面考虑。1)在图像融合过程中，图像融合过程是灰度值多对一的映射过程，在映射过程中存在不确定性；2)在成像过程中，由于PET和CT影像成像过程的复杂性，以及在成像过程中由于人体呼吸、血液流动和器官组织间相互重叠等无法避免的因素，导致图像存在混淆轮廓特征的噪声，使图像的模糊性被放大；3)在图像分解方法方面，原图像经过NSCT分解后，低频子带保留原图像的轮廓、背景等主要能量信息，映射的不确定关系也被保留，需要设计合理的融合规则对映射关系进行处理(蔡怀宇 等，2018)；4)基于模糊集理论的融合规则将整幅图像用模糊矩阵表示，通过一定的算法对模糊矩阵进行求解，能够有效解决图像融合过程中的模糊问题。高斯隶属度函数能够有效地描述低频子带包含的模糊信息(王艳 等，2019)。因此，使用高斯隶属度函数作为低频子带的自适应加权系数，采用基于模糊逻辑的自适应加权融合规则，计算式为

FL(i,j)=FLA(La,Lb)=

(3)

式中，FL(i,j)为低频子带融合图像；FLA(La,Lb)为基于模糊逻辑的自适应加权融合规则；La(i,j)和Lb(i,j)分别为图像A和B的低频图像；ωa(i,j),ωb(i,j)分别为La(i,j)和Lb(i,j)在(i,j)的高斯隶属度，ωb(i,j)的计算与ωa(i,j)相同，ωa(i,j)的计算式为

(4)

式中，μa为图像A的均值，σa为图像A的方差，k是高斯函数参数，参考范围为[1, 2.5]，一般取值为k=1.2；均值μb与μa的计算相同，方差σa与σb的计算相同，μa和σa的计算式为

(5)

(6)

式中，低频子带尺寸经过NSCT分解后和原图像相同，为M×N。

2.2 高频融合规则

高频子带的融合规则从3个方面考虑。1)考虑到高频子带自身的特点，原图像经过NSCT分解后，得到的8个高频子带包含原图像的组织器官的轮廓、边缘细节特征，具有多方向性和结构相似性，系数间存在较强的相关性；2)SSIM是衡量两幅图像相似性的指标，能够较好地反映高频子带系数之间的相关性，因此，采用平均结构相似性指标衡量两幅高频子带之间的系数相关性；3)肺部肿瘤的病变区域范围一般不超过100个像素点，所以采用基于区域的融合规则可以使病变区域的特性表现的更完整，区域方差能够表示局部区域的灰度值变化程度，方差越大，反映图像的细节信息就越丰富。因此，选取区域方差作为计算图像活性测度的依据。

高频子带融合规则设计具体为

FH=Piella(Ha,Hb)

(7)

式中，Piella()函数指基于Piella框架的融合规则。融合过程包含4个模块，分别为：Match模块(即计算两幅图像匹配测度的函数)、Active模块(即计算图像活性测度的函数)、决策模块和合成模块。

1)Match模块。本文使用平均结构相似性指标度量图像A和B高频系数之间的相似度Mab。首先利用滑动窗口将图像分块，考虑到窗口形状对分块的影响，简单的加窗会使映射矩阵出现不良的“分块”效应，因此使用高斯加权计算每一个窗口的均值、方差和协方差。考虑到PET图像的病灶约在100个像素以内，将滑动窗口大小设置为11 × 11，并逐像素遍历整幅图像，得到由局部SSIM组成的映射矩阵SSIMj(Ia,Ib)，最后对映射矩阵的所有局部SSIM指数取平均值。

设A和B图像分成M个图像块Ia,Ib，对M个图像块Ia,Ib的所有SSIMj(Ia,Ib)(j=1,2,…,M)求平均值，记为Mab。计算式为

(8)

式中，Hai和Hbi分别为图像A和B的第i个高频子带，Ia和Ib分别为图像A和B的图像块。SSIMj(Ia,Ib)(j=1,2,…,M)是局部SSIM组成的映射矩阵，图像块的SSIM值分别从亮度、对比度和结构相似性3个部分计算得到，计算式为

SSIM(Ia,Ib)=

[l(Ia,Ib)]α×[c(Ia,Ib)]β×[s(Ia,Ib)]γ

(9)

式中，α>0，β>0，γ>0是调整亮度、对比度和结构度所占权重的重要参数。设α=β=γ=1，Ia,Ib为Hai和Hbi相同位置提取的窗口大小为11 × 11的图像块，l(Ia,Ib)为Ia,Ib图像块的亮度相似性函数，c(Ia,Ib)为对比度相似性函数，s(Ia,Ib)为结构相似性函数，分别定义为

(10)

(11)

(12)

式中，C1，C2，C3为维持公式稳定的变量。设定C1=(K1×L)2，C2=(K2×L)2，C3=C2/2，K1=0.01，K2=0.03，L= 255。μIa，μIb，σIa，σIb，cov(Ia,Ib)分别为图像块Ia,Ib的均值、标准差和协方差。其中，μIa与μIb计算相同，σIa与σIb计算相同。μIa，σIa，cov(Ia,Ib)分别定义为

(13)

σIa=

(14)

cov(Ia,Ib)=

(Ib(x,y)-μIb)

(15)

式中，ω(x,y)(x,y=1,2,…,N)为11×11的高斯加权函数，标准差为1.5，N=11。

2)Active模块。图像相异度计算式为

αa=Active(Ha)

(16)

αb=Active(Hb)

(17)

经过NSCT分解得到的高频子带包含图像的边缘等细节信息，各个方向的高频子带系数间具有很强的区域相关性。因此，选取区域方差作为计算图像活性测度的依据，以更好地保留图像的清晰度。区域方差VARai与VARbi计算相同。VARai计算式为

(18)

式中，区域窗口尺寸设为3×3，μbi和μai为图像A和B高频子带Hai和Hbi的均值，μbi与μai计算相同。μai的计算式为

(19)

3)决策模块。使用匹配测度和活性测度共同构造决策模块，通过计算得到决策因子d对高频子带进行融合。决策因子d的计算式为

d=

(20)

设置阈值T= 0.5。当Mab

4)合成模块。

FHi(x,y)=d×Hai(x,y)+

(1-d)×Hbi(x,y)

(21)

式中，d为决策因子，对各个高频子带按照决策模块进行融合，即合成模块。

3 实验结果及分析

实验硬件环境为Intel(R)Core(TM)i7-9700CPU@3.00 GHz，64位Windows10操作系统；软件环境为MATLAB R2020a。

实验数据为某三甲医院的10例已配准的肺部肿瘤PET和CT影像，图像大小为365 × 365像素。

为验证本文方法的有效性进行了两组实验。实验1为5例CT肺窗和PET影像；实验2为5例CT纵膈窗和PET影像。对比实验选择4种像素级图像融合方法。方法1为NSCT-SSIM，用于证明LatLRR并行分解的有效性，分解方法为NSCT变换，高低频子带融合规则为本文方法一致；方法2为NSCT-CS(王惠群 等，2016)，使用NSCT变换，低频子带系数使用高斯加权平均融合规则，高频子带系数使用基于压缩感知(compression perception，CS)的融合规则；方法3为NSCT-CS-PCNN(陆惠玲 等，2017)，使用NSCT变换，低频子带系数使用PCNN融合规则，高频子带系数使用基于CS融合规则；方法4为LatLRR(Li和Wu，2022)，使用LatLRR变换，低频子带系数使用平均值融合规则，高频子带系数使用区域能量融合规则。

选用6种客观评价指标对融合图像质量进行评价，分别是：

1)平均梯度(average gradient，AG)(敬忠良 等，2007)，反映融合图像细节信息变化，平均梯度越大，图像细节越多，图像清晰度越高；

2)边缘强度(edge intensity，EI)(Zhu等，2016)，边缘强度越大，融合图像质量越高；

3)信息熵(information entropy，IE)(敬忠良 等，2007)，反映图像信息量，信息熵越大，融合图像的细节信息越多；

4)空间频率(spatial frequency，SF)(敬忠良 等，2007))，反映图像的整体清晰度；

5)标准差(standard deviation，STD)(敬忠良 等，2007)，衡量信息的丰富程度，标准差越大，图像的灰度级分布越分散，图像的信息量越多；

6种评价指标均是值越大，融合质量越好。

3.1 实验1：PET/CT肺窗

实验1为5组PET影像和CT肺窗影像的融合，结果如图3所示。

从图3可以看出，在肺部支气管清晰度显示方面，与4组对比实验相比，本文算法的融合结果更清晰地区分了支气管与周围软组织之间的位置关系，支气管区域的边缘效果表现更好。在病变区域显示方面，方法4的融合图像中的病变区域对比度较低，边缘模糊，对代谢信息显示效果较差，方法1、方法2、方法3和本文方法生成的融合图像能更清楚地显示病变区域的生理代谢信息，但方法3生成的融合图像伪影变多，病灶区域的阴影区域变大，对医生误诊带来一定隐患。实验结果表明，本文算法的融合性能较好。实验1融合图像指标评价结果如表1所示。

表1 PET/CT肺窗融合图像指标评价结果

图3 PET/CT肺窗融合结果

图4是PET/CT肺窗融合图像评价指标柱状图。可以看出，在AG、EI、SF、Qabf中，本文方法的融合图像优于其他4组对比实验，但在IE和STD数值中低于方法3，因为方法3提高了融合图像的对比度，从而保证IE和STD数值大，但方法3提高对比度时，图像中伪影变多，带来一定的误诊隐患。本文方法在增大对比度的情况下，有效避免了图像中伪影变多和病灶区域的阴影区域变大问题。虽然本文方法在这两个指标中的优势较弱，但在其他4种评价指标上取得了较好的结果。

图4 PET/CT肺窗融合图像评价指标柱状图

图5是可视化显示结果，是第3组图像的评价指标雷达图。可以看出，在AG、EI、SF、Qabf等客观评价指标中，本文方法生成的融合图像明显好于4组对比实验方法。

图5 第3组图像评价指标雷达图

3.2 实验2: PET/CT纵膈窗

实验2为5组CT纵膈窗和PET影像融合，融合结果如图6所示。

从图6可以看出，在骨骼、器官软组织清晰度显示方面，本文方法和方法3生成的融合图像中的骨骼、器官的轮廓信息的对比度更高。在病变区域显示方面，本文方法、方法1和方法2生成的融合图像更好地保留了PET影像中肿瘤病变区域的生理代谢信息。与方法1和方法2相比，本文方法对病变区域的定位线的显示有更好的效果。实验结果表明，本文方法能够较好地融合PET影像中的病灶信息和CT影像中的纵膈等组织、骨骼的轮廓细节信息。实验2融合图像指标评价结果如表2所示。

图6 PET/CT纵膈窗图像融合结果

图7是 PET/CT纵膈窗融合图像评价指标柱状图。从表2和图7可以看出，本文方法与4组对比实验生成的融合图像相比，在平均梯度、边缘强度、信息熵、空间频率和边缘保持度方面均有较大提升。方法3在标准差指标上(图7(e))更占优势，但在视觉效果上，方法3不能反映病变区域的位置。实验结果表明，本文方法生成的融合图像具有优势。

图7 PET/CT纵膈窗融合图像评价指标柱状图

表2 PET/CT纵膈窗融合图像指标评价结果

图8是第1组图像的评价指标雷达图。可以看出，相比4组对比实验，本文方法生成的融合图像在AG、EI、IE、SF和Qabf等客观评价指标中更占优势。

图8 第1组图像评价指标雷达图

4 结 论

针对CT图像对病灶的成像效果较差，尤其是浸润性肿瘤信息无法清晰显示，以及PET的骨组织的成像效果较差，传统像素级图像融合技术对PET/CT的融合图像不能较好地突出病灶区域信息，存在对比度低、边缘细节信息不能较好保留等问题，本文在Piella框架基础上，提出并行分解图像自适应融合模型，用于肺部肿瘤多模态图像融合。对已配准的CT和PET影像进行并行NSCT变换和LatLRR分解，经过NSCT变换后得到低频和高频子带，考虑到不同频率子带的特性以及图像自身的特点，低频子带使用自适应模糊逻辑融合规则，突出显示病灶的代谢信息；高频子带使用基于Piella框架的融合规则，保留原图像的边缘细节信息；对经过NSCT逆变换得到的融合图像与经过LatLRR分解得到的显著图相加，同时保留图像边缘细节信息和整体对比度，得到最终的融合图像。通过5组CT肺窗/PET和5组CT纵膈窗/PET的融合实验表明，与对比实验相比，本文方法生成的融合图像细节更清晰，视觉效果更好，在客观评价指标中取得较好的结果，有助于医生更快速、更精准地进行诊疗。

影像数据获取较难，本文方法仅在一个数据集上进行测试，未对其他尺寸的影像进行检测，因此，本文方法的普适性具有很大的提升潜力。此外，本文方法使用了并行分解方法，虽然提高了图像融合质量，但在运行中占用了较多的时间和内存，未来将研究时间复杂度和图像融合质量都更优的算法。