图像去噪是图像有关领域内的基本处理, 目标是去除图像中的噪声或最大程度地逼近原始图像, 作为进一步处理和分析图像的前导作业^[1]. 图像去噪结果的优劣除可用肉眼判定也可客观计算, 图像质量评价是判定处理后图像的一个基本方法, 在不同领域均有广阔的应用前景^[2]. 在医学诊断中, 根据医学图像判断病情, 质量较高的诊断图像对制定合理有效的治疗方案有至关重要的作用.

医学图像评价同样包括主观和客观两类. 主观评价基本以肉眼判定图像, 但各种因素影响, 评价结果往往不稳定, 且不便于集成和实现; 客观质量评价是目前研究的重要方向. 客观质量评价通常依据数学模型给出的量化指标来衡量图像质量. 峰值信噪比(Peak Signal Noise Ratio, PSNR)^[3]和结构相似度(Structure SIMilarity, SSIM)^[4]是常用的客观质量评价方法. 段影影等提出了用于医学图像质量评估的基于梯度加权的SSIM (gradient-Weighted SSIM, WSSIM)^[5]和基于梯度方向信息的SSIM (Gradient Direction SSIM, GDSSIM)^[6], 均是对SSIM算法的改进; 张剑华等^[7]提出基于结构显著性的医学图像质量评价, 在医学图像评价上取得了一定的成效; 笔者^[8]早期提出基于灰关联分析的PET重建图像评价方法, 取得了与PSNR和SSIM相当的评价效果. 陈宗哲等^[9]基于空间自相关性提出的自适应性非线性滤器在图像局部进行有选择的去噪计算, 有效地改进了传统滤器的缺陷, 提高了图像去噪效果. 但医学图像质量评价在临床应用上还需要进一步研究和优化.

PET是将生物生命代谢中必须的物质标记上短寿命的放射性核素标记物注入人体, 通过对该物质在代谢中的聚集, 来反映生命代谢活动的情况, 达到诊断的目的, 是无创的、动态的、定量观测活体生理变化的医学成像工具之一^[10,11]. 但目前重建出高质量图像的方法一直是PET的研究热点. 进来发展的基于期望最大化EM(Expectation Maximization)迭代算法的最大似然期望最大化MLEM(Maximum-Likelihood EM)算法或有序子集最大期望值OSEM(Ordered Subsets EM)算法^[12], 在重建过程中考虑了测量资料的统计性质, 是目前最常用的PET图像重建算法, 相对于滤波反投影FBP(Filtered Back Projection)重建法有很多优点^[13]. 但迭代重建图像质量并不是随着迭代的次数越多越好. 重建图像初始随着迭代次数增加, 图像细节逐渐恢复, 随着迭代次数的增加, 图像噪声逐渐上升, 某次迭代后, 重建图像的质量随着迭代的进行反而开始下降^[14]. 目前临床一般根据PET仪器制造商设定的迭代次数或依实际经验作为停止依据, 所获得图像质量还有进一步提升的空间. 本研究探讨PET重建图像噪声与质量间的关系, 分析PET重建图像噪声随迭代次数的变化情况, 通过评估噪声对图像质量的影响, 探索可能最适当的迭代次数或停止依据.

1 PET图像与噪声

注入体内的放射性核素在衰变过程放出正电子, 正电子与人体组织中的负电子相结合发生湮灭, 产生一对反向传播的511 keV伽马射线, 被环绕在人体周围的探测器组接收, 生成正弦图^[13]. 使用重建算法处理正弦图而获取诊断图像, 了解核素在体内的分布及代谢情况, 获得患者身体组织的病理信息^[13]. PET扫描仪一般使用OSEM迭代重建算法, 其存在的主要问题是无法得到收敛解且随着迭代次数增加引入了大量的噪声. 图1显示使用OSEM重建的不同迭代次数Huffman图像(图像大小128×128, OSEM 01表示第一次迭代结果). 可见随着迭代次数增加, 图像质量越来越好, 纹理更加清晰, 但噪声也逐渐增加.

2 方法 2.1 MLEM最大似然最大期望值和OSEM有序子集最大期望值

MLEM是基于像素的迭代算法, 每个像素值被看作待估计参数, 通过不断迭代更新, 使似然函数逼近最大, 得到极大似然参数估计值. 由于MLEM算法计算量大, 收敛缓慢, 在1994年, Hudson和Larkin提出了OSEM^[15], 使用有序子集来加速MLEM算法. 其思想是将投影数据分成多个小组, 称为有序子集, 每次迭代处理其中一个小组的数据, 多个小组可以同时进行. 相比MLEM加速效果明显, OSEM算法被广泛应用在各种PET机器上.

PET重建使用MLEM算法的迭代公式为:

$f_j^{(k + 1)} = \frac{{f_j^{(k)}}}{{\sum\nolimits_{i = 1}^I {{a_{ij}}} }}\sum\nolimits_j^I {\frac{{{a_{ij}}{y_j}}}{{\sum\nolimits_j^J {{a_{ij}}f_j^k} }}} $

(1)

其中, y表示观测到的投影数据. $A = \left\{ {{a_{ij}}} \right\}$ 是系统矩阵, ${a_{ij}}$ 表示从像素j发出的光子被第i条投影线所在探测单元接收到的概率, f为重建图像.

假设重建图像大小为n×n, 360个投影角度, 每角度下n条投影射线. 将投影数据按投影角度划分成T个有序子集, S_I表示第I个子集, I =1, 2, …, T. OSEM算法迭代公式为:

$f_j^{(l + 1)} = f_j^{(l)}\frac{{\sum\nolimits_{j \in {s_l}} {\frac{{{P_i}{\omega _{ij}}}}{{\sum\nolimits_{j = 1}^N {{\omega _{ij}}x_j^{(l)}} }}} }}{{\sum\nolimits_{j \in {s_l}} {{\omega _{ij}}} }}$

(2)

2.2 蒙特卡罗模拟和假体图像

本研究使用操作系统Ubuntu14.04LST, 计算机硬件为配置Intel I5 CPU, 16 G内存和独显的图形工作站, 使用GATE蒙特卡罗模拟软件(geant4版本4.10) 模拟Siemens ECAT PET扫描仪. 该扫描仪内径43 cm, 576个LSO材质的晶体围成一圈, 晶体大小为4×4×20 mm, 晶体与晶体间的距离为0.5 mm, 每个探测器模块(Detector Block)为12×12个晶体, 48个模块围成一圈, 总共有4个圈. 为缩短时间模拟时只模拟射源及探测, 不记录散射及衰减. Hoffman和Utah 图像分别模拟至900 Mega次计数(counts)后停止.

研究中使用的MLEM、OSEM算法利用IRT (Image Reconstruction Toolbox, http://www.eecs.umich.edu/~fessler/code/)软件在x64 PC利用Matlab R2009a运行.

2.3 图像质量评价指标 2.3.1 PSNR峰值信噪比图像评价指标

PSNR是基于像素点的图像质量评价指标. 过去研究表明^[3,4], PSNR图像质量评价方法是可有效的评价图像质量且已经广泛的使用. PSNR的计算公式如下所示:

$PSNR = 10\lg \left( {\frac{{{a_{\max }}^2}}{{MSE}}} \right)$

(3)

$MSE = \frac{1}{{MN}}\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^N {{{\left[ {a(i,j) - \mathop a\limits^ \wedge (i,j)} \right]}^2}} } $

(4)

MSE(Mean Square Error)称为均方误差, 是衡量“平均误差”的一种较方便的方法, 可以评价数据的变化程度. $a{{(i,j)}}$ 和 $\mathop a\limits^ \wedge {{(i,j)}}$ 分别代表原始图像与重建图像中对应像素点的灰度值, MN为图像的总像素, ${a_{\max }} = {2^\ell } - 1$ 通常就是图像的灰度级, 常取值为255.

2.3.2 SSIM结构相似度图像评价指标

SSIM与PSNR不同之处是考虑了图像的结构信息, 即像素与其周边像素间的关联信息. SSIM将像素的强度指针分割为标准的亮度和对比度, 从图像组成的角度将图像质量评价建模为三个分量的比较: 亮度、对比度和结构^[4,5].

$SSIM({x},{y}) = l{({x},{y})^\alpha } \cdot c{({x},{y})^\beta } \cdot s{({x},{y})^\gamma }$

(5)

其中, x和y分别代表两幅比较图像, l、c、s分别代表图像的亮度比较函数、对比度比较函数和结构比较函数, α、β、γ>0, 是亮度、对比度和结构信息的权重系数, 本研究中取值为1.

2.4 噪声评估方法: STD标准偏差

Olsen等人曾经对比了6种噪声评估方法, 发现通过估算均值滤器过滤的图像与原图像相减所得差值图像的标准偏差(Standard Deviation, STD)最能代表图像的高斯噪声量^[16]. 通过比较不同迭代次数的重建图像以均值滤器处理后再与原图像运算所得差值图像的标准偏差平均值, 可用以说明图像噪声的变化情况. 我们将不同迭代次数的PET图像分别以3×3均值滤波器过滤后, 与原图像相减得到相应的差值图像, 再以3×3的滑动窗口计算差值图像的平均STD.

3 结果

经由OSEM算法重建, 获取不同迭代次数的Huffman重建图像, 计算重建图像与标准Huffman图像的PSNR、SSIM和STD, 分析其质量和噪声的变化情况.

3.1 Huffman重建图像的PSNR

图2展示了重建图像的PSNR随着迭代次数的变化情况. 迭代一开始, PSNR显著增加, 重建图像质量先增加后, 质量再开始下降, 当迭代至第36次之后(峰值), PSNR开始减小.

3.2 Huffman重建图像的SSIM

图3表明, OSEM迭代一开始, 根据SSIM的计算结果, 重建图像质量不断增加, 第22次以后, 质量开始缓慢下降, 结果与PSNR相同.

3.3 Huffman重建图像的STD

随着OSEM迭代次数的增加, 重建图像噪声增加. 如图4所示, 迭代一开始, STD快速增加, 第22次迭代后达到中等水平, 之后噪声缓慢上升.

由PSNR和SSIM的计算结果可以看出, OSEM在迭代至某一次时产生的图像质量较好. 此后, PSNR、SSIM随着迭代次数增加而降低. 由STD计算结果可以看出, PSNR与SSIM取得最大值时, 噪声也处于一个中等合适水平, 随后, 噪声继续增加, 而图像质量开始退化.

4 讨论 4.1 PET图像质量与噪声的关系

均值滤波主要过滤了图像中的高频部分, 也就是噪声, 所以滤波后的图像主要包含的是有用的质量信息, 理论上滤波后的图像质量更好. 图5、图6展示了经过均值滤波后的图像质量指标PSNR及SSIM, 滤波后图像的计算值均明显低于未经滤波的图像(图中实线所示), 且最优值分别出现在第97次(PSNR)及48次(SSIM), 晚于未滤波的图像第36次(PSNR)及22次(SSIM). 这说明滤波不仅仅过滤了噪声, 也影响了图像质量. 噪声滤除前与后的图像质量计算的变化, 说明了噪声也影响到PSNR、SSIM对于图像质量的计算. 最优图像出现晚于未滤波的图像, 这说明了噪声缓慢上升将造成最优图像错误的计算. 如果我们能够精准的估算噪声并且有效的滤除, 将可能获得真正最优迭代图像, 作为停止迭代依据. 均值滤波是一种最简易的滤波方法, 这种滤波方法是否能有效的将PET图像噪声滤除还需要进一步改进研究. STD计算结果也只是粗略的噪声估算. 未来改进PET图像噪声滤除以及精准的估算噪声是一个重要的研究方向.

4.2 MLEM重建Utah图像的质量与噪声

为了进一步验证PET图像质量、噪声的变换规律, 我们采用另一核医常用Utah图像来进行模拟验证.

图7展示了Utah标准图像和MLEM迭代图像, 图8计算了Utah的标准偏差STD, 图9比较了Utah迭代图像与均值滤波图像(虚线)的SSIM. 可以看出, Utah重建图像质量随着迭代次数不断提高, 第22次之后开始退化, 噪声持续增加; 均值滤波图像质量指标低于原重建图像, 最优图像较迟出现, 这些现象与Huffman图像估算结果一致.

这表明PET重建图像有以下规律: 第一, 重建图像质量在迭代初始阶段迅速提高, 某次迭代之后开始退化. 第二, 重建图像噪声持续增加, 某次迭代之后增加平缓. 第三, 噪声对重建图像质量影响较大. 第四, 均值滤波最优图像出现较晚, 迭代需要晚于PSNR或SSIM最优迭代次数后停止.

5 结语

PET不同迭代次数图像的图像质量和噪声与迭代次数都有较大关联. 迭代一开始, 质量不断提高, 某次迭代之后质量开始退化, 噪声不断增加. SSIM评价标明, Huffaman和Utah重建图像质量在第22次取得最大值, 均值滤波后, 最大SSIM在Huffaman和Utah中并不一致, 出现在第48次(Utah出现在第36次), 但都在 一个相对合理的范围内. 噪声滤除前后图像质量的变化, 说明噪声影响了图像质量的计算, 噪声缓慢上升将造成最优图像计算错误. 如果我们能够精准的估算噪声并且有效的滤除, 将可能获得真正最优迭代图像, 作为停止迭代依据. 均值滤波是一种最简易的滤波方法, 是否能有效的将PET图像噪声滤除还需要进一步改进研究, STD计算结果也只能是粗略的噪声估算. PET图像噪声的精确估算和有效滤除是未来研究的重要方向.

实际诊断中, 迭代应在PSNR或SSIM最优迭代次数之后某次停止, 需要医生根据具体情况做出判断. 我们需要结合医生的判断, 研究适合PET噪声滤除的方法, 探索基于噪声考虑的迭代停止规则和图像质量评价途径, 以找到合适的PET迭代停止规则.