编码板成像系统MLEM算法优化

李汉平 王锋 艾宪芸



编码板成像系统MLEM算法优化

李汉平 王锋 艾宪芸

（国民核生化灾害防护国家重点实验室防化研究院 北京 102205）

为增强γ辐射编码成像系统中最大似然期望最大化(Maximum Likelihood Expectation Maximization, MLEM)算法对噪声的抑制，提高算法重建质量，基于互补编码板成像相减去除噪声的思想，对原有MLEM算法进行改进，提出了MLEM互补算法，并引入修正因子对MLEM互补算法的收敛程度进行控制。通过模拟及实验方式验证了MLEM互补算法的有效性，给出了不同情况下MLEM互补算法达到最优值时修正因子的拟合曲线。结果表明，MLEM互补算法可有效抑制噪声，提高重建图像质量。修正因子的拟合曲线可以有效确定MLEM算法的最佳值。

编码板，图像重建，噪声抑制，最大似然期望最大化法，蒙特卡罗模拟

在γ射线编码板成像系统中，重建算法的好坏直接影响重建图像的质量。对于修正均匀冗余阵列(Modified Uniformly Redundant Arrays, MURA)编码模式下的编码板成像系统，研究人员对一些重建算法进行过对比研究[1]，而最大似然期望最大化(Maximum Likelihood Expectation Maximization, MLEM)算法由于能对噪声进行很好的抑制，且在投影数据不完全时可以对图像进行很好的重建而成为主流算法。MLEM算法的缺点是：在编码板投影图像所含噪声较大时，随着迭代次数的增加，噪声对重建图像的影响会相应放大。目前对于MLEM算法的改进主要集中在核医学领域，如有序子集期望最大化(Ordered Subset Expectation Maximization, OSEM)算法、子集序列期望最大化(Subset Sequence Expectation Maximization, SSEM)算法以及计算调控有序子集期望最大化(Count Regulated Ordered Subset Expectation Maximization, CROSEM)算法等[2‒3]。由于编码板成像与医学成像在数据采集方式上的不同，医学成像中的相关算法不能直接应用于编码板成像。在编码板成像中，有关MLEM算法改进的论述较少。

对于编码板成像系统，噪声大小直接影响其图像质量。为了抑制噪声对重建图像的影响，通过理论证明[4]，可以用两个互补的编码孔经重建所得图像相减来减小噪声对重建图像的影响。本文基于该思想对编码板成像中的MLEM算法进行改进，提出MLEM互补算法。该算法在传统MLEM算法公式中加入编码板修正因子，在迭代算法中进行互补编码板相减计算。MLEM互补算法中最佳修正因子值的选取与编码板线性衰减系数与厚度的乘积有关。通过蒙特卡罗模拟给出算法最佳修正因子的t拟合公式，并进行实验验证。实验结果证明，与传统MLEM算法相比，MLEM互补算法在成像时可以更有效地抑制噪声。采用模拟所得-拟合公式，可以给出不同情况下MLEM互补算法最佳值。

1 算法简介

1.1 MLEM算法介绍

MLEM算法是基于泊松统计模型的最大贝叶斯后验概率的图像复原方法，该算法由Shepp等[5]于1982年提出。MLEM算法主要分为两步：1) 计算完全投影数据的似然函数在测定数据及当前参数估计值下的期望；2) 求出使期望最大化的参数估计值。

在γ辐射编码成像中MLEM算法[6‒8]表示为：

1.2 基于互补码板的MLEM算法

基于互补编码板相减去除噪声的思想，在进行MLEM迭代时，对迭代公式(1)中的编码函数(,) 进行修正。令(,)=′(,)+h″(,)，其中：′(,)为正常编码函数；″(,)为′(,)逆时针旋转90°所得编码函数；为算法修正因子，取值范围为(−0.8,0)。将(,)带入迭代公式，可得：

在计算第次放射源强度分布的迭代估计值和编码函数所计算得出的投影估计值时使用互补编码板相减去噪的方式对噪声进行抑制。由式(2)可知，由于在计算时对进行了噪声抑制，使得实际测定值与修正理论投影值的比值对重建图像的修正效果更好。将该值与修正后的(,)做相关运算，使得迭代算法再次对噪声进行抑制，进一步减少噪声对重建图像的影响。

2 模拟及实验结果分析

2.1 模拟条件及模拟结果分析

使用蒙特卡罗模拟软件MCNP5 (Monte Carlo N Particle Transport Code 5)对γ辐射编码板成像过程进行模拟，得到实验所需编码板投影数据。建立的γ辐射编码板成像模型示意图如图1所示。编码板以11×11扩展模式的MURA编码板为基础，几何尺寸为50mm×50mm×10mm，采用方孔模式，孔径2mm，材料选用铅；点源位于视野中心轴线，距编码板500mm处，能量为662keV；探测器材料选用锗酸铋(Bi2O3-GeO2, BGO)闪烁晶体，几何尺寸为50mm×50mm×15mm，距编码板20mm，计数时，将BGO晶体模型按100×100×1进行网格划分，并采用*F8计数卡对每个网格内沉积能量进行记录。通过蒙特卡罗模拟所得理想点源图像及其投影图像如图2所示。

图1 编码板模型示意图

图2 理想点源(a)和点源投影(b)图像

采用蒙特卡罗模拟所得编码板投影数据，分别通过MLEM算法和MLEM互补算法对图像进行重建。图3为MLEM算法与MLEM互补算法(=−0.6)在迭代次数=5时的重建图像。由图3可知，两种算法都可对模拟编码图像进行有效重建。采用MLEM互补算法对模拟图像进行重建，其对噪声的抑制效果更好，重建图像更接近模拟的点源图像。通过直接观察法可知，MLEM互补算法的成像效果要优于传统的MLEM算法。

对重建图像进行评估时，主观评价往往很难判断算法的优劣。为进一步验证算法的有效性，在对重建图像的质量进行评价时，采用归一化均方误差(Normalized Root Mean Square Error, NRMSE)以及相关系数(Correlation Coefficient, CC)这两种定量方法对重建图像的质量进行评价。用归一化均方误差衡量算法的接近程度，相关系数衡量重建图像与原图像相似程度[9]。NRMSE与CC评价标准的值和定义分别为：

表1 MLEM算法和MLEM互补算法NRMSE与CC值

2.2 实验条件及实验结果分析

为验证模拟结果的正确性，通过实验对算法的有效性进行验证。实验中编码板以11×11扩展的MURA编码模式为基础，编码板几何尺寸分别为50mm×50mm×10mm、50mm×50mm×20mm、50mm×50mm×30mm，孔型为方孔，孔径2mm，材料选用铅；探测器由50mm×50mm×10mm的BGO阵列晶体，耦合滨松H8500位置灵敏光电倍增管组成，探测器距编码板20mm。点源采用强度3.7×104Bq的137Cs源，位于视野中心轴线，距编码板500mm处。

采用实验所得编码板投影数据，分别通过MLEM算法和MLEM互补算法对图像进行重建。图4为MLEM算法和MLEM互补算法在迭代次数=5时对三种厚度编码板实测投影重建图像。由图4(a)可知，随着编码板厚度的增加，MLEM算法重建图像的背景噪声越来越小，图像质量逐步提高。由图4(b)可知，采用MLEM互补算法进行图像重建时，在编码板厚度相同的情况下，MLEM互补算法的重建结果比MLEM算法的重建结果对噪声的抑制更明显。

图4 MLEM算法(a)与MLEM互补算法(b)实验数据重建图像

表2为MLEM算法和MLEM互补算法在迭代次数=5时对三种厚度编码板实测投影重建图像的NRMSE和CC值。由表2可知，在编码板厚度相同的情况下，MLEM互补算法的NRMSE和CC值皆优于传统的MLEM算法，实验结果与模拟结果相一致。与MLEM算法相比，采用MLEM互补算法对编码图像进行重建，可以更有效地抑制噪声，提高重建图像的质量。

表2 MLEM算法和MLEM互补算法的NRMSE、CC值

2.3值对MLEM互补算法的影响

通过MLEM互补算法对上述实验及模拟投影数据进行重建时，=−0.6。为研究不同值对MLEM互补算法重建效果的影响，选取不同值对实验所得厚度=10mm、=20mm、=30mm三种厚度编码板的实测投影值进行互补MLEM算法重建。图5为三种厚度下互补MLEM算法的NRMSE和CC值随修正因子值的变化。

图5 三种厚度铅制编码板NRMSE (a)和CC (b)值随a变化

由图5可知，随着值在(−0.8,0)区间内不断增大，重建图像的质量先随之提高，当达到某一特定值后，重建图像的质量开始下降。对于MLEM算法存在最佳值，使得改进后MLEM算法效果最优。编码板厚度为1 cm、2 cm、3 cm时，最佳值分别为−0.7、−0.61、−0.55。可见随着编码板准直器厚度的增加，MLEM互补算法修正因子的最佳值逐渐减小，最佳值的选取与物质阻止本领有关。由γ射线在物质中的衰减规律，采用上述编码板模型，模拟不同t值下，算法最佳值的变化。图6为最佳值随t的变化。由图6可知，最佳值随t值一同增加，当t值达到某一特定值时，值的变化趋于平缓。计算所得t拟合公式为：

=−0.6+0.12421−0.01252()2(6)

图6 最佳a值随mt值变化拟合曲线

表3为三种厚度下铅制编码板实测投影数据所得最佳值与模拟数据拟合公式计算所得最佳值。

表3 三种厚度铅制编码板实验及模拟数据最佳a值

由表3可知，由于噪声影响，实验数据所得值比模拟数据所得的要大。无法通过模拟数据所得t拟合公式直接给出实测情况下MLEM互补算法最佳值。在使用t拟合公式给出实际成像系统所需最佳值时，需对模拟所得最佳值乘以相应比例系数。

记实验值与模拟值的比值为比例系数，根据t拟合公式变化趋势，对随t值变化进行拟合。图7为比例系数随t值的变化。

图7 比例系数N随最佳mt值变化拟合曲线

计算所得t拟合公式为：

=1.0156+0.4462−0.0651()2(7)

在求MLEM互补算法的最佳值时，可以通过模拟所得t拟合公式计算结果乘以得到。

为验证上述结论，采用1cm钨板进行实验，图8为钨板的NRMSE和CC值随值的变化。由图8可知，此时MLEM互补算法最佳值为0.65。在入射γ射线能量为662keV时，钨的线性衰减系数=1.84cm−1[10]；将=1cm带入拟合公式(6)，得出模拟情况下算法最佳值为0.41，由拟合公式(7)，此时为1.6。最终用于1cm钨制编码板的最佳=1.6×0.41=0.65。与实验测得数据所得最佳值相一致。在实际应用中，可以用上述方法确定MLEM互补算法最佳值。

图8 1cm厚钨制编码板NRMSE (a)和CC (b)值随a变化

3 结语

实验证明，将互补编码板相减去噪思想应用于MLEM算法中是可行的。通过对蒙特卡罗模拟所得投影数据及实测所得投影数据的重建结果可知，相比于MLEM算法，在编码板成像中，MLEM互补算法可以有效抑制噪声的影响，提高图像重建质量；MLEM互补算法中的最佳值的选取与t有关；通过模拟得出t拟合公式所得模拟最佳值与t拟合公式所得的乘积，可以得到实测数据的最佳值，模拟结果与实验结果相一致。采用该方法可以确定MLEM互补算法最佳值；在编码板成像系统中，采用蒙特卡罗模拟方法对其进行分析具有一定的可行性。

1 赵翠兰, 陈立宏, 李勇平. MURA编码辐射成像系统的解码方法[J]. 核技术, 2014, 37(8): 080401. DOI: 10.11889/j.0253-3219.2014.hjs.37.080401.

ZHAO Cuilan, CHEN Lihong, LI Yongping. Decoding process of a radiation imaging system using MURA coded aperture collimator[J]. Nuclear Techniques, 2014, 37(8): 080401. DOI: 10.11889/j.0253-3219.2014.hjs.37.080401.

2 杨娟, 王明泉, 石浪, 等. OSEM重建算法及其改进算法的研究和比较[J]. 计算机工程与设计, 2015, 36(9): 2524‒2526. DOI: 10.16208/j.issn1000-7024.2015.09.040. YANG Juan, WANG Mingquan, SHI Lang,. Research and comparison on OSEM and its improved reconstruction algorithms[J]. Computer Engineering and Design, 2015, 36(9): 2524‒2526. DOI: 10.16208/j. issn1000-7024.2015.09.040.

3 金永杰, 马天予. 核医学仪器与方法[M]. 哈尔滨: 哈尔滨工程大学出版社, 2010.

JIN Yongjie, MA Tianyu. Nuclear medical instrument and method[M]. Harbin: Harbin Engineering University Press, 2010.

4 Barrett H H, Swindell W. 放射成像[M]. 庄天戈, 周颂凯, 译. 北京: 科学出版社, 1988.

Barrett H H, Swindell W. Radiation imaging[M]. ZHUANG Tiange, ZHOU Songkai, trans. Beijing: Science Press, 1988.

5 SheppL A,Vardi Y. Maximum likelihood reconstruction for emission tomography[J]. IEEE Transactions on Medical Imaging, 1982, 1(2):113‒122. DOI: 10.1109/TMI.1982.4307558.

6 张斌, 王英, 艾宪芸, 等. 基于约束的MLEM图像重建算法[J]. 原子能科学技术, 2014, 48(增1): 668‒672. DOI: 10.7538/yzk.2014.48.S0.0668.

ZHANG Bin, WANG Ying, AI Xianyun,. MLEM image reconstruction algorithm based on physical constraints[J]. Atomic Energy Science and Technology, 2014, 48(Suppl 1): 668‒672. DOI: 10.7538/yzk.2014.48. S0.0668.

7 洪俊杰. γ相机中编码孔经准直器设计及数字图像重建[D]. 湖北: 华中科技大学, 2006. HONG Junjie. Design of coded aperture collimator and digital image reconstruction in gamma-ray camera[D]. Hubei: Huazhong University of Science and Technology, 2006.

8 Mu Z P, Liu Y H. Aperture collimation correction and maximum-likelihood image reconstruction for near-field coded aperture imaging of single photon emission computerized tomography[J]. IEEE Transactions on Medical Imaging, 2006, 25(6): 701‒711. DOI: 10.1109/ TMI.2006.873298.

9 何佳伟, 刘东升, 桂志国. 可变有序子集PML算法在PET中的应用[J]. 中北大学学报(自然科学版), 2010, 31(6): 646‒650. DOI: 10.3969/j.issn.1673-3193.2010.06. 021.

HE Jiawei, LIU Dongsheng, GUI Zhiguo. Application of modified subsets PML algorithm to PET image reconstruction[J]. Journal of North University of China (Natural Science Edition), 2010, 31(6): 646‒650. DOI: 10.3969/j.issn.1673-3193.2010.06.021.

10 方杰. 辐射防护导论[M]. 北京: 原子能出版社, 1991.FANG Jie. Introduction to radiation protection[M]. Beijing: Atomic Press, 1991.

Algorithm optimization of MLEM in coded aperture imaging system

LI Hanping WANG Feng AI Xianyun

(State Key Laboratory of NBC Protection for Civilian, Research Institute of Chemical Defense, Beijing 102205, China)

Background: Ingamma-ray imager, reconstruction algorithm directly affects the quality of the reconstructed image. The maximum likelihood expectation maximization (MLEM) iteration algorithm is widely used inmodified uniformly redundant arrays(MURA) coded aperture for the reason of satisfactory performance of suppressing noise, improving signal noise ratio (SNR) and reducing distortion. But MLEM iteration algorithm also has shortcomings, like amplifies the noise with the increase of the iterative times. Purpose: This study aims to improve the ability of suppressing noise for MLEM iteration algorithm, and increase precision and resolution of coded aperture imaging system.Methods:Based on the de-noising method of complementary coded aperture, a correction factorfor MLEM algorithm modification is proposed to control the convergence rate of the complementary MLEM iteration algorithm.Both Monte Carlo simulation and experimental date are used to verify the effectiveness of the complementary MLEM iteration algorithm. Results: The results of Monte Carlo simulation and experiment prove the complementary MLEM iteration algorithm is efficient and practicable in coded aperture image. Its efficiency is relatedto the modifying factor. The fitting formula of-inthe complementary MLEM iteration algorithm is=−0.6+0.12421−0.01252()2. Conclusions: The coded aperture gamma camera is capable of imaging gamma-rays instantly. It can form radioactive 2D-map of different radionuclides. The reconstruction algorithm is important for coded aperture gamma camera. With the development of coded aperture, more and more study of algorithm modification will focus on it.

Coded aperture, Reconstruction of image, Suppress noise, MLEM, Monte Carlo simulation

LI Hanping, male, born in 1992, graduated from Shanghai Jiao Tong University in 2014, master student, major in radiation protection and environmental protection

2016-11-20, accepted date: 2016-12-14

TL812

10.11889/j.0253-3219.2017.hjs.40.020404

李汉平，男，1992年出生，2014年毕业于上海交通大学，现为硕士研究生，辐射防护与环境保护专业

2016-11-20，

2016-12-14