易利群 盛玉霞 柴 利

正电子发射断层成像(Positron emission to-mography,PET)是一种强大的功能性成像模式,其通过注入特定的放射性示踪剂来观察组织内部的分子水平活动.PET 图像可以显示分子水平的示踪剂分布,从而可以对生理和药理过程进行定量测量.但是,由于各种物理退化因素和检测到的光子数量有限,PET 图像空间分辨率低、噪声程度高,影响了疾病检测和诊断的准确性,因此有必要提高PET图像质量.

对PET 图像后处理可提升图像质量,图像去噪和去模糊是常见的PET 图像后处理方法.PET图像的噪声复杂,没有精确的模型,通常认为具有泊松和高斯混合噪声的特性[1],直接应用传统的去噪方法不能得到好的去噪效果.高斯滤波是一种简单的去噪算法,常用于PET 图像去噪,但其在去噪过程中会平滑重要的图像结构信息.许多学者提出了针对PET 图像的去噪方法,包括双边滤波[2]、各向异性传播和非局部均值(Non-local means,NLM)滤波[1,3-4]、引导滤波[5]等.文献[1]提出了一种改进的基于NLM 的动态PET 图像去噪方法;文献[4] 提出了一种新的结合多尺度小波变换和NLM 的滤波方法;文献[5]研究了一种基于磁共振成像(Magnetic resonance imaging,MRI)引导的脑部PET 图像滤波方法;对PET 图像去模糊问题,Song 等[6]提出了结合MR (Magnetic resonance)图像联合熵先验的PET 图像超分辨率重建方法,相对于传统的全变分和二次罚函数等方法具有更好的去模糊效果;文献[7]提出了一种基于生成对抗网络的自监督PET 超分辨率重建方法;稀疏表示字典学习也常用于提升PET 图像质量[8-9].Li 等[8]通过稀疏表示对低分辨率PET 图像进行复原;Wang等[9]应用半监督字典学习由低剂量PET 图像来预测标准剂量图像;文献[10]通过多层典型相关分析(Canonical correlation analysis,CCA)方法来估计标准剂量PET 图像;全变分方法应用于直接重建模型中,以改善PET 图像重建质量[11-12].Zhang 等[11]提出了一种共边缘加权全变分正则化的PET-MRI联合重建方法,该方法可以同时提高两类图像的重建质量;文献[12]提出了一种基于时空全变分正则化的动态PET 图像重建方法,重建效果优于传统的最大似然期望最大化(Maximum likelihood expectation maximization,MLEM)方法.以上方法都有效地提高了PET 图像空间分辨率.

随着PET/MR 一体化成像设备的应用,高分辨率的MRI 先验信息可用来进一步提高PET 图像去噪效果.Yan 等[5]研究了基于MRI 的引导滤波方法,该方法给出结构MR 图像和功能性PET 图像之间的局部线性模型,并将部分体积效应(Partial volume correction,PVC)融入到模型中,同时对PET 图像去噪和部分体积校正.Cheng 等[13]提出了一种将MRI 与迭代高级约束的局部反投影方法(Iterative highly constrained back-projection local region,IHYPR-LR)相结合的去噪方法,该方法将MR 图像作为IHYPR-LR 方法进行迭代的初始值,在去噪的同时增强了图像的结构和边缘.文献[14]针对脑部成像问题提出了两种不同的PETMR 核方法来修复PET 图像.文献[15]综述了机器学习在PET 图像处理中的应用.

近年来,深度学习广泛应用于医学图像处理中[16-18],也为PET 图像处理提供了新的研究方向[19].文献[16-17]详细综述了深度学习在医学影像中的最新进展;文献[19]综述了基于深度学习的PET 图像重建方法;文献[20]提出了基于预训练VGG (Visual Geometry Group)深度网络的PET图像去噪方法,使用人工生成的数据训练网络,在实际的PET 图像上效果一般;文献[21]通过深度学习方法研究不同级别的低剂量PET 图像的先验信息,用来估计高质量的PET 图像.最近,学者们将生成对抗网络应用于PET 图像处理[7,22-23].Song等[7]提出了基于双重生成对抗网络的自监督PET超分辨率重建算法,无需配对训练数据,具有更广泛的适用性.Gong 等[22]使用混合2D 和3D 网络学习图像先验,复原低剂量PET 图像.文献[23]提出了用于PET 图像去噪的参数传递Wasserstein 生成对抗网络,在保持图像保真度的同时可以有效地抑制噪声.Wang 等[24]提出了一种新颖的基于3D自动上下文的局部自适应多模态融合生成对抗网络模型,用来由低剂量的PET 图像和MR 图像合成高质量的PET 图像,该方法优于深度网络中传统的多模态融合方法和PET 估计方法.但医学图像存在样本少且标注困难等问题,深度学习方法在医学图像上泛化性能差、网络结构缺乏通用性、网络训练成本大、计算复杂度高.非深度学习方法在样本少的情况下有其特有的优势,值得进一步研究.

图信号处理(Graph signal processing,GSP)是近年兴起的高维非规则化数据分析处理新方法,其借助代数图论和谱图理论来处理高维加权图上的信号[25-27].在图像处理领域,可以通过建立非局部图像像素间的图连接关系将图像转化为图信号,图包含的结构信息可以很好地保留图像边缘特性.图信号处理在图像滤波、压缩和复原等应用中取得了较好的效果[28].文献[29]设计了一种低通图滤波器,在图傅里叶域对PET 图像滤波,去噪效果相比于高斯滤波和NLM 都有所提高,但该方法不能清晰地保留病灶点.Hammond 等[30-31]提出了图小波变换(Graph wavelet transform,GWT)方法,并给出了GWT 的切比雪夫多项式逼近快速实现算法,文献[31]详细介绍了图小波变换的基础理论和快速实现方法.相对于图傅里叶变换(Graph Fourier transform,GFT),图小波引入多尺度分析,在一些去噪应用中获得了好的效果[30-32].目前还没有应用图小波来对PET 图像进行去噪的研究.

如上所述,PET 图像空间分辨率低、噪声大.已有的PET 去噪方法不能在去噪的同时很好地保留病灶信息.如何既保留PET 图像的病灶点,又能有效滤除噪声是PET 图像去噪的难点问题.针对此问题,本文提出了一种新的基于图小波的PET图像去噪方法.图滤波的去噪性能依赖于图加权邻接矩阵的构造,直接用低分辨率PET 图像构造邻接矩阵噪声较大,会影响去噪效果.本文给出了结合MRI 结构信息和PET 图像的邻接矩阵构造方法,既能反映MR 图像的结构信息,又能突出PET图像的病灶特点,提高图滤波性能.该方法首先将多帧正弦图相加后重建,得到PET 合成图像,然后将PET 合成图像和MR 图像融合,由融合图像构造邻接矩阵,再对动态PET 图像进行图小波变换去噪.仿真数据实验表明,本文方法与单独的图滤波、图小波变换,以及其他基于模型的PET 图像去噪方法相比具有更好的去噪性能,与VGG 深度神经网络等基于学习的方法去噪效果相当,但不需要大量的训练数据和参数调整.

本文的主要创新性包括: 1)提出了新的基于图小波的PET 图像去噪方法;2) 给出了一种融合MRI 结构信息和PET 图像的构图方法;3)本文方法能更好地去噪和保持病灶信息,优于图滤波方法和结合MRI 的传统PET 图像去噪方法;4)与深度学习方法相比,本文方法不受标注样本限制,无需训练和迭代,计算复杂度低.

1 理论基础

1.1 图谱理论

在图信号处理中,信号x=[x1,x2,···,xN]T通常被转换为一个在无向加权图G={V,ε,A}上的信号;其中V={v1,v2,···,vN}是图上顶点的集合,ε⊆V ×VT是一组边,A∈RN×N为加权邻接矩阵.如果有一条连接了顶点i和顶点j的边ε=(i,j),则Ai,j ＞0 表示这条边的权重;否则Ai,j=0.得到邻接矩阵A后,可以定义图拉普拉斯矩阵L

其中,D是对角度矩阵,D的对角元素dii是矩阵A第i行元素之和.

图傅里叶变换将传统的傅里叶变换推广到图谱域,给出了图信号顶点域与图谱域之间的傅里叶变换关系.将图拉普拉斯矩阵分解成特征值和特征向量,令λℓ和χℓ=[χℓ(0),χℓ(1),···,χℓ(N-1)]T,ℓ=0,1,···,N -1分别为图拉普拉斯矩阵L对应的特征值和特征向量.特征值表示频率,特征向量形成图傅里叶变换基,任何图信号都可以通过图拉普拉斯矩阵的特征向量表示.令f=[f(1),f(2),···,f(N)]T为一个图信号,则它的第i个元素f(i) 是图G的第i个顶点上的信号,信号f的图傅里叶变换定义为[25]

其中,f代表顶点域信号,其对应的图傅里叶逆变换为

1.2 图小波变换

图小波是定义在加权图顶点域数据上的小波变换.图小波变换需要定义一个图小波核函数,通过在图谱域的伸缩操作构造不同的图谱域小波簇,并在图谱域调制信号.在文献[30]中,作者定义图小波函数算子为Tg=g(L),信号f的图小波变换在图谱域的表达式为

其中,g是图小波核函数,满足g(0)=0 和limx→∞g(x)=0的条件,它是图谱域的带通滤波器,类似于传统小波变换母函数的傅里叶变换,根据图傅里叶逆变换(3)得到图小波逆变换

给定信号f,其图小波系数和尺度系数定义如下[30]:

1.3 图小波变换逆变换

图小波变换通过切比雪夫多项式逼近得到图小波近似系数,其可以看作是将一个大小为N的输入信号f映射为大小为N(J+1) 的系数c,即c=W f,其中f∈RN×1为原始信号,J为小波尺度数目.W为矩阵算子,且W存在很多左逆M,使得MW f=f.已知系数c,通过求解方程 (W*W)f=W*c可得到W的伪逆 (W*W)-1W*[30],并将求解出的伪逆作为GWT 的逆[31].为了加快计算速度,采用共轭梯度算法,并应用切比雪夫多项式逼近方法得到重建信号.关于图小波更详细的介绍,可参考文献[30-31].

2 方法

2.1 融合MRI 信息的动态PET 图像去噪方法

将多个PET 正弦图短时间帧相加后重建得到合成图像,可以改善计数统计信息并减少噪声[33].合成帧图像尽管丢失了时间信息,但不同时间帧的空间信息可以很好地保留在合成帧中,用于提高短时间帧的重建质量.在本文中,对于有病灶点的PET 图像,将多帧正弦图相加得到合成的正弦图,再通过最大似然期望最大化(Maximum likelihood expectation maximization,MLEM)方法[34]重建得到PET 合成图像Ic,以保留病灶点信息.为了保留图像更多的结构信息,本文提出了一种新的结合MRI 信息的融合图像构造方法,在此基础上对PET 图像进行图小波变换,用来去除PET 图像中的噪声.通过硬阈值方法将MR 图像与PET 合成图像进行融合,融合方法为

其中,Ic(m,n),IMRI(m,n),IF(m,n) 分别表示PET合成图像、MR 图像和融合图像在点 (m,n) 处的像素值,Threshold为阈值参数.对合成图像中大于阈值的点,融合图像取PET 合成图像的像素值,这部分保留了PET 图像中的病灶信息;对小于阈值的点,融合图像取MR 图像的像素值,以保留结构信息.PET 图像病灶点的像素值相对于对应位置的MR 图像像素值偏大,通过设置适当的阈值,可以保留病灶信息.新的融合图像既保留了PET图像的病灶信息,也包含有MR 图像的部分结构信息.

在融合图像上构造邻接矩阵A,将融合图像中每个像素点当作是图G上的一个顶点,顶点i和顶点j之间的边缘权重由式(14)的阈值高斯核权重函数定义,即

其中,fin表示PET 图像Iin的向量形式,切比雪夫多项式阶数M=51,j表示图小波第j个尺度,j=1,2,···,J.对PET 图像Iin进行第1 次图小波多尺度分解后,得到表示低频的图谱尺度系数1,0和表示高频的图小波系数1,j,通过保留图小波域的低频段尺度系数1,0,并丢弃高频段的图小波系数1,j来达到去噪的效果.同时,为了尽可能保留在高频段部分丢弃的图像细节,保留图小波系数的第1 个尺度系数1,1,进行第2 次图小波多尺度分解

通过两次图小波多尺度分解得到两组图小波变换系数,保留两次分解得到的低频图谱尺度系数1,0和2,0.将系数组合=[1,0,2,0,0] 进行图小波逆变换,通过共轭梯度法求解式(19),得到去噪后的PET 图像[30]

本文方法的整个实现过程如下:

2.2 参数选择

本文方法有9 个关键参数: 图像块的大小n,最近邻居数目k,平滑参数θ,图像块相似度影响参数η,小波尺度数目J,尺度核函数h,小波核函数g,小波尺度参数s和阈值参数Threshold.经过反复实验,选取相对较优的参数设置,其中尺度核函数h,小波核函数g和小波尺度参数s设置与文献[30]中相同,其他参数具体设置见表1.

表1 本文方法参数设置Table 1 Parameter setting in this paper

3 实验与分析

我们在两个PET 图像仿真数据集上进行实验,将本文方法和MRI 引导滤波方法[5]、MRI 引导的IHYPR-LR 方法[13]、直接由PET 合成图像构图的图滤波[29]和图小波四种方法进行比较.其中图小波方法是直接由PET 合成图像构图,没有融合MRI信息,相当于阈值为0.其他参数设置与本文融合MRI 信息的方法相同.为公平起见,将对比方法[5,13]的实验参数都调整到最优.

3.1 评价指标

为了评估和比较各种去噪方法,使用信噪比(Signal-to-noise ratio,SNR)和均方根误差(Root mean squared error,RMSE) 来评估去噪后的PET 图像质量.

其中,xi表示原始PET 图像体素i的强度值,yi表示去噪后的PET 图像体素i的强度值.

3.2 模拟数据生成

本文使用3D 大脑模型通过计算机模拟Brain-Web 的幻像[35],该模型由脊椎液、白质和灰质组成,并选择两个基于18F-FDG 的葡萄糖代谢的组织隔室模型作为动力学模型.本文从正常MR 图像中选择了第90 个切片,从异常MR 图像中选择了第90个切片,并采用文献[36]的方法产生区域时间活动曲线(Time activity curve,TAC),在注射示踪剂时开始动态PET 扫描.动态PET 数据由60 min内的24 个时间帧组成: 4×20 s,4×40 s,4×60 s,4×180 s 和 8×300 s.为了生成带噪声的动态PET图像,本文采用PET 成像领域常用的方法[1],首先通过正向投影生成原始的动态PET 正弦图,然后将60 min 内产生 7×107个光子总数的泊松过程应用于原始的无噪声正弦图中,最后通过60 次MLEM 方法迭代重建得到无病灶动态PET 图像和单病灶动态PET 图像.实验中PET 图像大小为175×175,MR 图像大小为 1 81×217,通过Brain-Web 中的配准参数得到已配准的MR 图像.

3.3 实验结果与分析

3.3.1 无病灶数据实验结果与分析

第1 个实验从光子数为 7×107时生成的无病灶动态PET 大脑图像中,选择第6 帧、第12 帧、第18 帧和第24 帧进行去噪效果演示.MR 图像如图1(a),从动态PET 图像中得到的合成图像如图1(b)所示,该PET 图像没有病灶点,直接使用MR 图像构图.观察图1(a)可知,MR 图像含有更多的结构信息,可以帮助我们对动态PET 图像进行去噪.无病灶情况下PET 图像去噪实验结果如图2 所示,无病灶情况下结合MRI 的PET 图像去噪方法比较结果见表2,直接由PET 合成图像构图的方法与本文方法比较结果见表3.由图2 可以看出,与其他方法相比,本文方法得到的PET 去噪图像具有更好的视觉效果.同时在客观评估指标方面,与结合MRI 的PET 图像去噪方法、直接由 PET合成图像构图的去噪方法相比,本文方法的SNR值最大,RMSE 值最小,说明有更好的去噪效果.

图1 两种图像示例Fig.1 Example of two kinds of images

图2 无病灶PET 图像去噪结果Fig.2 Denoising results of normal PET images

表2 无病灶情况下结合MRI 的PET 图像去噪方法比较Table 2 Comparison of PET image denoising methods incorporated with MRI on the normal dataset

表3 无病灶情况下由PET 合成图像构图的方法与本文方法比较Table 3 Comparison of the methods of constructing graph by PET composite image with the proposed method on the normal dataset

3.3.2 单病灶数据实验结果与分析

第2 个实验从光子数为 7×107时生成的单病灶动态PET 大脑图像中,选择第6 帧、第12 帧、第18 帧和第24 帧进行去噪效果演示.BrainWeb中相应的MR 图像如图3(a)所示,动态PET 合成图像如图3(b)所示,由MR 图像和PET 合成图像得到的融合图像如图3(c)所示.观察图3(c)可以看出,融合图像既保留了MRI 图像的结构信息,也保留了合成图像中的病灶点信息.单病灶情况下5种方法去噪后的PET 图像如图4 所示,由图4 可以看出,本文方法在去噪的同时还保留了更多的结构信息和病灶信息,有更好的视觉效果.各种方法的SNR 和RMSE 值如表4 和表5 所示.可以看出,本文方法的SNR 值最大,RMSE 值最小.为进一步评估本文方法的去噪性能,我们分别画出了噪声图像和5 种方法去噪后的PET 图像残差图,如图5所示.观察图5 可知,本文方法得到的PET 图像与原图更相近,并且病灶点和边缘细节保留得更完整.

图3 单病灶图像Fig.3 Single-lesion images

图4 单病灶PET 图像去噪结果Fig.4 Denoising results of single-lesion PET images

图5 单病灶PET 图像去噪残差图Fig.5 Denoising residual map of single-lesion PET images

表4 单病灶情况下结合MRI 的PET 图像去噪方法比较Table 4 Comparison of PET image denoising methods incorporated with MRI on the single-lesion dataset

表5 单病灶情况下由PET 合成图像构图的方法与本文方法比较Table 5 Comparison of the methods of constructing graph by PET composite image with the proposed method on the single-lesion dataset

为了验证不同噪声级别下本文方法的有效性,在单病灶情况下,本文设置了光子数分别为7×108、7×109两种情况来模拟生成不同噪声级别的PET 图像,各种方法的SNR 和RMSE 值如表6 所示.可以看出,本文方法在不同噪声程度下,SNR值最大,RMSE 值最小,有更好的去噪效果.

表6 单病灶PET 图像在不同光子数时各种去噪方法比较Table 6 Comparison of different denoising methods for single-lesion PET images with different photon numbers

为了定量比较不同方法在病灶点区域的重建效果,用病灶点区域对比度恢复系数(Contrast recovery coefficient,CRC)与白质区域(背景区域)标准差(Standard deviation,STD)作为衡量指标.CRC 越接近1,说明病灶点区域的重建去噪效果越好;STD 值越小,说明背景区域的噪声越少.CRC与STD 的计算式为[20]

本文对比高斯(Gaussian)滤波、MRI 引导滤波(GuideMRI)[5]、局部反投影方法(HYPRMRI)[13]、图小波(GWT)和本文方法(GWTMRI)对病灶点的恢复效果.指标参数的设置与文献[20]相同.用以上5 种方法对MLEM 重建图像去噪,MLEM 共迭代120 次,每迭代24 次计算一个CRC-STD 值,改变迭代次数得到的CRC-STD 曲线如图6 所示.可以看出,本文方法相对于高斯滤波、GuideMRI、HYPRMRI、GWT 有更好的CRC-STD 折中值.在同样水平的STD 下,本文方法有更大的CRC 值.因没有算法实现代码,不能直接与基于深度学习的方法对比,我们间接比较各种方法的CRC-STD 图.从文献[20]中的深度学习方法与高斯滤波的CRCSTD 曲线对比图可以看出,文献[20]中的方法和本文方法相对于高斯滤波性能都有很大程度的提升.值得指出的是,该仿真示例仅表明了本文所提方法具有与基于深度学习方法相当的去噪性能.但本文方法的优势在于可解释性好,不需要大规模数据训练.

图6 不同方法病灶点CRC 与背景区域STD 曲线图Fig.6 CRC-STD curves of the different denoising methods for single-lesion PET images

4 结束语

提出了一种融合MRI 信息的图小波变换方法用于动态PET 图像去噪.本文将MR 图像与PET合成图像通过硬阈值方法进行融合,在融合图像上构造图拉普拉斯矩阵,最后通过图小波变换对动态PET 图像去噪.在两种不同的BrainWeb 仿真数据上进行实验,实验结果表明,与结合MRI 的PET图像去噪方法以及基于图滤波的去噪方法相比,本文方法在图像视觉效果和评估指标方面的表现都是最优的.这说明本文提出的融合MRI 信息的图小波方法可以更大程度地降低噪声,同时保留更多的PET 图像病灶信息,具有更好的去噪性能.本文使用的图小波函数是文献中给定的常用函数,根据PET 图像的特点,能否设计性能更好的专门针对医学图像的图小波函数,值得进一步研究.