肺癌是造成全球死亡人数最多的癌症(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作为活性测度，使融合图像更符合人类视觉系统的感知特性。

2006 年,da Cunha 等人(2006) 在 contourlet 变换的基础上进行改进,提出NSCT变换。NSCT使用 两个非下采样滤波器组,即非下采样金字塔滤波器组( non-subsampled pyramid filter banks, NSPFB) 和非下采样方向滤波器组( non-subsampled directional filter banks, NSDFB)。NSCT变换取消 contourlet变换中的下采样操作,避免 contourlet 变换过程中容易出现的伪吉普斯现象。NSCT变换分解得到的高低 频子带尺寸与原图像相同。单层NSCT分解过程如下:首先,原图像经过 NSPFB 分解为一个低频子带和一个高频子带。然后,对高频子带使用 NSDFB 进行$K$级方向分解, 得到$2^k$个方向子带。$J$级NSCT分解可以得到1个低频子带和$\sum\limits_{j=1}^J 2^{k_j}$个高频子带。其中, $k_j$为尺度$j$下的方向分解级数。NSCT分解示意图如图 1所示。

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

$ \boldsymbol{X}=\boldsymbol{Z} \boldsymbol{X}+\boldsymbol{L} \boldsymbol{X}+\boldsymbol{E} $

(1)

式中, $\boldsymbol{X}$是输人图像, $\boldsymbol{Z}$是低秩系数, $\boldsymbol{L}$是显著性系数; $\boldsymbol{Z X}$和$\boldsymbol{L X}$分别表示低秩部分和显著部分; $\boldsymbol{E}$表示稀疏噪声。LatLRR问题再次进行最优化求解，其最终表达式为

$ \begin{gathered} \min _{\boldsymbol{Z}, \boldsymbol{L}, \boldsymbol{E}}\|\boldsymbol{Z}\|_*+\|\boldsymbol{L}\|_*+\lambda\|\boldsymbol{E}\|_1 \\ \text { s. t. } \quad \boldsymbol{X}=\boldsymbol{Z} \boldsymbol{X}+\boldsymbol{L} \boldsymbol{X}+\boldsymbol{E} \end{gathered} $

(2)

式中, $\|\cdot\|_*=\sum \sigma_i$表示矩阵的核范数, $\sigma_i$为矩阵的奇异值, $\|\cdot\|_1$表示L1范数; $\lambda$表示平衡系数且$\lambda>0$。

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

1) $\left[\boldsymbol{L}_a, \boldsymbol{H}_a\right]={NSCTdec}(\boldsymbol{A}); \left[\boldsymbol{L}_b, \boldsymbol{H}_b\right]=NSCTdec (\boldsymbol{B})$。图像$\boldsymbol{A}$和$\boldsymbol{B}$ 是已配准的CT和PET图像, $NSCTdec$ () 是NSCT变换。$\boldsymbol{L}_a$和$\boldsymbol{L}_b$ 分别是图像$\boldsymbol{A}$和$\boldsymbol{B}$ 经过NSCT变换后得到的低频子带, $\boldsymbol{H}_a=$ $\left\{H_{a i}(x, y) \mid i=1, 2, \cdots, 8\right\}$和$\boldsymbol{H}_b=\left\{H_{b i}(x, y) \mid i=\right.$ $1, 2, \cdots, 8\}$ 分别是图像$\boldsymbol{A}$和$\boldsymbol{B}$ 经过NSCT变换后得到的8个高频子带。

2) $\left[\boldsymbol{L r r}_a, \boldsymbol{S}_a, \boldsymbol{E}_a\right]={LatLRR}(\boldsymbol{A})$; $\left[\boldsymbol{L r r}_b, \boldsymbol{S}_b\right.$, $\left.\boldsymbol{E}_b\right]={LatLRR}(\boldsymbol{B}) 。{LatLRR}(\;)$是LatLRR分解。$\boldsymbol{L r r}_a, \boldsymbol{S}_a, \boldsymbol{E}_a$和$\boldsymbol{L r r}_b, \boldsymbol{S}_b$和$\boldsymbol{E}_b$分别是图像$\boldsymbol{A}$和$\boldsymbol{B}$经过LatLRR分解后得到的低秩部分、显著部分和噪声部分。

3) $\boldsymbol{F}_{\mathrm{L}}=F L A\left(\boldsymbol{L}_a, \boldsymbol{L}_b\right) 。F L A($) 函数是基于模糊逻辑自适应的低频子带融合规则(fuzzy logic adaptive fusion rules, FLA), $\boldsymbol{F}_{\mathrm{L}}$是低频子带融合图像。

4) $\boldsymbol{F}_{\mathrm{H}}={Piella}\left(\boldsymbol{H}_a, \boldsymbol{H}_b\right) ; M_{a b}={Match}\left(\boldsymbol{H}_a, \right. $$ \left.\boldsymbol{H}_b\right) ; \alpha_a={Active}\left(\boldsymbol{H}_a\right) ; \alpha_b={Active}\left(\boldsymbol{H}_b\right)$。

计算高频子带融合图像

$ \boldsymbol{F}_{\mathrm{H}}=d \times \boldsymbol{H}_a+(1-d) \times \boldsymbol{H}_b $

$Piella$ () 函数是基于Piella框架的高频子带融合规则, $Match$ () 函数是计算匹配测度的方法, $M_{a b}$是两幅高频子带之间的相似值, $Active$() 是计算活性测度的函数, $\alpha_a$和$\alpha_b$是两幅图像$\boldsymbol{A}$和$\boldsymbol{B}$经过$Active$函数计算得到的图像相异度, $d$是决策因子, $\boldsymbol{F}_{\mathrm{H}}$是$Piella$函数的返回值, 表示高频子带的融合图像。

5) $\overline{\boldsymbol{F}}={InverseNSCT}\left(\boldsymbol{F}_{\mathrm{L}}, \boldsymbol{F}_{\mathrm{H}}\right) 。$ $InversenSCT$ () 是NSCT逆变换, $\overline{\boldsymbol{F}}$是$InverseNSCT$函数的返回值, 表示融合图像。

6) $\boldsymbol{F}=\overline{\boldsymbol{F}}+\boldsymbol{S}_a+\boldsymbol{S}_b \circ \boldsymbol{F}$是图像$\boldsymbol{A}$和$\boldsymbol{B}$的最终融合图像。

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

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

$ \begin{gathered} \boldsymbol{F}_{\mathrm{L}}(i, j)=F L A\left(\boldsymbol{L}_a, \boldsymbol{L}_b\right)= \\ \frac{\omega_a}{\omega_a+\omega_b} \boldsymbol{L}_a(i, j)+\frac{\omega_b}{\omega_a+\omega_b} \boldsymbol{L}_b(i, j) \end{gathered} $

(3)

式中, $\boldsymbol{F}_{\mathrm{L}}(i, j)$为低频子带融合图像; $F L A\left(\boldsymbol{L}_a, \boldsymbol{L}_b\right)$为基于模糊逻辑的自适应加权融合规则; $\boldsymbol{L}_a(i, j)$和$\boldsymbol{L}_b(i, j)$分别为图像$\boldsymbol{A}$和$\boldsymbol{B}$的低频图像; $\omega_a(i$, $j), \omega_b(i, j)$分别为$\boldsymbol{L}_a(i, j)$和$\boldsymbol{L}_b(i, j)$在$(i, j)$的高斯隶属度, $\omega_b(i, j)$的计算与$\omega_a(i, j)$相同, $\omega_a(i, j)$的计算式为

$ \omega_a(i, j)=\exp \left\{-\frac{\left(\boldsymbol{L}_a(i, j)-\boldsymbol{\mu}_a\right)^2}{2\left(k \sigma_a\right)^2}\right\} $

(4)

式中, $\boldsymbol{\mu}_a$为图像$\boldsymbol{A}$的均值, $\sigma_a$为图像$\boldsymbol{A}$的方差, $k$是高斯函数参数, 参考范围为$[1, 2.5]$, 一般取值为$k=1.2$; 均值$\boldsymbol{\mu}_b$与$\boldsymbol{\mu}_a$的计算相同, 方差$\sigma_a$与$\sigma_b$的计算相同, $\boldsymbol{\mu}_a$和$\sigma_a$的计算式为

$ \boldsymbol{\mu}_a=\frac{1}{M \times N} \sum\limits_{i=1}^M \sum\limits_{j=1}^N \boldsymbol{L}_a(i, j) $

(5)

$ \sigma_a=\frac{\sum\limits_{i=1}^M \sum\limits_{j=1}^N\left(\boldsymbol{L}_a(i, j)-\boldsymbol{\mu}_a\right)^2}{M \times N} $

(6)

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

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

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

$ \boldsymbol{F}_{\mathrm{H}}={Piella}\left(\boldsymbol{H}_a, \boldsymbol{H}_b\right) $

(7)

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

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

设$\boldsymbol{A}$和$\boldsymbol{B}$图像分成$M$个图像块$\boldsymbol{I}_a, \boldsymbol{I}_b$, 对$M$个图像块$\boldsymbol{I}_a, \boldsymbol{I}_b$的所有$\boldsymbol{S S I M}_j\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)(j=1, 2$, $\cdots, M)$求平均值, 记为$M_{a b}$。计算式为

$ M_{a b}={Match}\left(\boldsymbol{H}_a, \boldsymbol{H}_b\right)=\frac{1}{M} \sum\limits_{j=1}^M {\boldsymbol{SSIM}}_j\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right) $

(8)

式中, $\boldsymbol{H}_{a i}$和$\boldsymbol{H}_{b i}$分别为图像$\boldsymbol{A}$和$\boldsymbol{B}$的第$i$个高频子带, $\boldsymbol{I}_a$和$\boldsymbol{I}_b$分别为图像$\boldsymbol{A}$和$\boldsymbol{B}$的图像块。$\boldsymbol{S S I M}_j\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)(j=1, 2, \cdots, M)$是局部$\boldsymbol{S S I M}$组成的映射矩阵, 图像块的$\boldsymbol{S S I M}$值分别从亮度、对比度和结构相似性3个部分计算得到, 计算式为

$ \begin{gathered} {\boldsymbol{SSIM}}\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)= \\ {\left[l\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)\right]^\alpha \times\left[c\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)\right]^\beta \times\left[s\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)\right]^\gamma} \end{gathered} $

(9)

式中, $\alpha>0, \beta>0, \gamma>0$是调整亮度、对比度和结构度所占权重的重要参数。设$\alpha=\beta=\gamma=1$, $\boldsymbol{I}_a, \boldsymbol{I}_b$为$\boldsymbol{H}_{a i}$和$\boldsymbol{H}_{b i}$相同位置提取的窗口大小为$11 \times$ 11的图像块, $l\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)$为$\boldsymbol{I}_a, \boldsymbol{I}_b$图像块的亮度相似性函数, $c\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)$为对比度相似性函数, $s\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)$为结构相似性函数, 分别定义为

$ l\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)=\frac{2 \boldsymbol{\mu}_{I a} \boldsymbol{\mu}_{I b}+C_1}{\boldsymbol{\mu}_{I a}^2+\boldsymbol{\mu}_{I b}^2+C_1} $

(10)

$ c\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)=\frac{2 \sigma_{I a} \sigma_{I b}+C_2}{\sigma_{I a}^2+\sigma_{I b}^2+C_2} $

(11)

$ s\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)=\frac{{cov}\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)+C_3}{\sigma_{l a} \sigma_{l b}+C_3} $

(12)

式中, $C_1, C_2, C_3$为维持公式稳定的变量。设定$C_1=\left(K_1 \times L\right)^2, C_2=\left(K_2 \times L\right)^2, C_3=C_2 / 2, K_1=$ $0.01, K_2=0.03, L=255_{\circ} \boldsymbol{\mu}_{I a}, \boldsymbol{\mu}_{I b}, \sigma_{I a}, \sigma_{l b}$, ${cov}\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)$分别为图像块$\boldsymbol{I}_a, \boldsymbol{I}_b$的均值、标准差和协方差。其中, $\boldsymbol{\mu}_{I a}$与$\boldsymbol{\mu}_{I b}$计算相同, $\sigma_{I a}$与$\sigma_{I b}$计算相同。$\boldsymbol{\mu}_{I a}, \sigma_{I a}, {cov}\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)$分别定义为

$ \boldsymbol{\mu}_{I a}=\frac{1}{M \times N} \sum\limits_{x=1}^M \sum\limits_{y=1}^N \omega(x, y) \boldsymbol{I}_a(x, y) $

(13)

$ \begin{gathered} \sigma_{I a}= \\ \sqrt{\frac{1}{M \times N-1} \sum\limits_{x=1}^M \sum\limits_{y=1}^N \omega(x, y)\left(\boldsymbol{I}_a(x, y)-\boldsymbol{\mu}_{I a}\right)^2} \end{gathered} $

(14)

$ \begin{gathered} {cov}\left(\boldsymbol{I}_a, \boldsymbol{I}_b\right)= \\ \frac{1}{M \times N-1} \sum\limits_{x=1}^M \sum\limits_{y=1}^N \omega(x, y)\left(\boldsymbol{I}_a(x, y)-\boldsymbol{\mu}_{I a}\right) \times \\ \left(\boldsymbol{I}_b(x, y)-\boldsymbol{\mu}_{I b}\right) \end{gathered} $

(15)

式中, $\omega(x, y)(x, y=1, 2, \cdots, N)$为$11 \times 11$的高斯加权函数, 标准差为$1.5, N=11$。

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

$ \alpha_a={Active}\left(\boldsymbol{H}_a\right) $

(16)

$ \alpha_b={Active}\left(\boldsymbol{H}_b\right) $

(17)

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

$ {VAR}_{a i}(x, y)=\frac{\sum\limits_{x=1}^M \sum\limits_{y=1}^N\left(\boldsymbol{H}_{a i}(x, y)-\boldsymbol{\mu}_{a i}\right)^2}{M \times N} $

(18)

式中, 区域窗口尺寸设为$3 \times 3, \boldsymbol{\mu}_{b i}$和$\boldsymbol{\mu}_{a i}$为图像$\boldsymbol{A}$和$\boldsymbol{B}$高频子带$\boldsymbol{H}_{a i}$和$\boldsymbol{H}_{b i}$的均值, $\boldsymbol{\mu}_{b i}$与$\boldsymbol{\mu}_{a i}$计算相同${ }_{\circ} \boldsymbol{\mu}_{a i}$的计算式为

$ \boldsymbol{\mu}_{a i}=\frac{1}{M \times N} \sum\limits_{x=1}^M \sum\limits_{y=1}^N \boldsymbol{H}_{a i}(x, y) $

(19)

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

$ \begin{align} & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \mathit{d}\text{=} \\ & \left\{ \begin{array}{*{35}{l}} 1\quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{M}_{ab}}<T\text{ 且 } \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \mathit{VA}{{\mathit{R}}_{\mathit{ai}}}(\mathit{x, y})>\mathit{VA}{{\mathit{R}}_{\mathit{bi}}}(\mathit{x, y} ) \\ 0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{M}_{ab}}<T\ \text{且 } \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \mathit{VA}{{\mathit{R}}_{\mathit{ai}}}(\mathit{x, y})<\mathit{VA}{{\mathit{R}}_{\mathit{bi}}}(\mathit{x, y}) \\ \frac{\mathit{VA}{{\mathit{R}}_{\mathit{ai}}}(\mathit{x, y})}{\mathit{VA}{{\mathit{R}}_{\mathit{ai}}}(\mathit{x, y})+\mathit{VA}{{\mathit{R}}_{\mathit{bi}}}(\mathit{x, y})}\quad {{\mathit{M}}_{\mathit{ab}}}\ge \mathit{T} \\ \end{array} \right. \\ \end{align} $

(20)

设置阈值$T $= 0.5。当$M_{ab}＜T $时，两个高频子带的灰度值差异明显，选择区域方差较大的高频子带系数作为融合系数；$ M_{ab}≥T $时，两个对应高频子带的相似度较高，使用区域方差自适应加权方法计算高频子带的融合系数。

4) 合成模块。

$ \begin{gathered} \boldsymbol{F}_{H i}(x, y)=d \times \boldsymbol{H}_{a i}(x, y)+ \\ (1-d) \times \boldsymbol{H}_{b i}(x, y) \end{gathered} $

(21)

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

5) $\overline{\boldsymbol{F}}= InversenNSCT \left(\boldsymbol{F}_{\mathrm{L}}, \boldsymbol{F}_{\mathrm{H}}\right)$。其中, $InverseNSC T$() 表示NSCT逆变换, $\boldsymbol{F}_{\mathrm{L}}$表示低频部分融合图像, $\boldsymbol{F}_{\mathrm{H}}$表示高频部分融合图像, $\overline{\boldsymbol{F}}$为融合图像。

6) $\boldsymbol{F}=\overline{\boldsymbol{F}}+\boldsymbol{S}_a+\boldsymbol{S}_b$。其中, $\boldsymbol{S}_a$和$\boldsymbol{S}_b$为$\boldsymbol{A}$和$\boldsymbol{B}$的显著图, $\boldsymbol{F}$为最终的融合图像。

实验硬件环境为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) 边缘保持度(Qabf)，值越高表示融合图像保留原图像的边缘信息越丰富(Bondžulić和Petrović, 2008)。

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

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

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

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

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

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

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

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

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

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

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