編碼板成像系統MLEM算法優(yōu)化

李漢平 王鋒 艾憲蕓

?

編碼板成像系統MLEM算法優(yōu)化

李漢平 王鋒 艾憲蕓

（國民核生化災害防護國家重點實驗室防化研究院 北京 102205）

為增強γ輻射編碼成像系統中最大似然期望最大化(Maximum Likelihood Expectation Maximization, MLEM)算法對噪聲的抑制，提高算法重建質量，基于互補編碼板成像相減去除噪聲的思想，對原有MLEM算法進行改進，提出了MLEM互補算法，并引入修正因子對MLEM互補算法的收斂程度進行控制。通過模擬及實驗方式驗證了MLEM互補算法的有效性，給出了不同情況下MLEM互補算法達到最優(yōu)值時修正因子的擬合曲線。結果表明，MLEM互補算法可有效抑制噪聲，提高重建圖像質量。修正因子的擬合曲線可以有效確定MLEM算法的最佳值。

編碼板，圖像重建，噪聲抑制，最大似然期望最大化法，蒙特卡羅模擬

在γ射線編碼板成像系統中，重建算法的好壞直接影響重建圖像的質量。對于修正均勻冗余陣列(Modified Uniformly Redundant Arrays, MURA)編碼模式下的編碼板成像系統，研究人員對一些重建算法進行過對比研究[1]，而最大似然期望最大化(Maximum Likelihood Expectation Maximization, MLEM)算法由于能對噪聲進行很好的抑制，且在投影數據不完全時可以對圖像進行很好的重建而成為主流算法。MLEM算法的缺點是：在編碼板投影圖像所含噪聲較大時，隨著迭代次數的增加，噪聲對重建圖像的影響會相應放大。目前對于MLEM算法的改進主要集中在核醫(yī)學領域，如有序子集期望最大化(Ordered Subset Expectation Maximization, OSEM)算法、子集序列期望最大化(Subset Sequence Expectation Maximization, SSEM)算法以及計算調控有序子集期望最大化(Count Regulated Ordered Subset Expectation Maximization, CROSEM)算法等[2?3]。由于編碼板成像與醫(yī)學成像在數據采集方式上的不同，醫(yī)學成像中的相關算法不能直接應用于編碼板成像。在編碼板成像中，有關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互補算法的成像效果要優(yōu)于傳統的MLEM算法。

對重建圖像進行評估時，主觀評價往往很難判斷算法的優(yōu)劣。為進一步驗證算法的有效性，在對重建圖像的質量進行評價時，采用歸一化均方誤差(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值皆優(yōu)于傳統的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)區(qū)間內不斷增大，重建圖像的質量先隨之提高，當達到某一特定值后，重建圖像的質量開始下降。對于MLEM算法存在最佳值，使得改進后MLEM算法效果最優(yōu)。編碼板厚度為1 cm、2 cm、3 cm時，最佳值分別為?0.7、?0.61、?0.55?？梢婋S著編碼板準直器厚度的增加，MLEM互補算法修正因子的最佳值逐漸減小，最佳值的選取與物質阻止本領有關。由γ射線在物質中的衰減規(guī)律，采用上述編碼板模型，模擬不同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 金永杰, 馬天予. 核醫(yī)學儀器與方法[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年畢業(yè)于上海交通大學，現為碩士研究生，輻射防護與環(huán)境保護專業(yè)

2016-11-20，

2016-12-14