李漢平 王鋒 艾憲蕓

?

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

李漢平 王鋒 艾憲蕓

（國(guó)民核生化災(zāi)害防護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室防化研究院 北京 102205）

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

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

在γ射線編碼板成像系統(tǒng)中，重建算法的好壞直接影響重建圖像的質(zhì)量。對(duì)于修正均勻冗余陣列(Modified Uniformly Redundant Arrays, MURA)編碼模式下的編碼板成像系統(tǒng)，研究人員對(duì)一些重建算法進(jìn)行過(guò)對(duì)比研究[1]，而最大似然期望最大化(Maximum Likelihood Expectation Maximization, MLEM)算法由于能對(duì)噪聲進(jìn)行很好的抑制，且在投影數(shù)據(jù)不完全時(shí)可以對(duì)圖像進(jìn)行很好的重建而成為主流算法。MLEM算法的缺點(diǎn)是：在編碼板投影圖像所含噪聲較大時(shí)，隨著迭代次數(shù)的增加，噪聲對(duì)重建圖像的影響會(huì)相應(yīng)放大。目前對(duì)于MLEM算法的改進(jìn)主要集中在核醫(yī)學(xué)領(lǐng)域，如有序子集期望最大化(Ordered Subset Expectation Maximization, OSEM)算法、子集序列期望最大化(Subset Sequence Expectation Maximization, SSEM)算法以及計(jì)算調(diào)控有序子集期望最大化(Count Regulated Ordered Subset Expectation Maximization, CROSEM)算法等[2?3]。由于編碼板成像與醫(yī)學(xué)成像在數(shù)據(jù)采集方式上的不同，醫(yī)學(xué)成像中的相關(guān)算法不能直接應(yīng)用于編碼板成像。在編碼板成像中，有關(guān)MLEM算法改進(jìn)的論述較少。

對(duì)于編碼板成像系統(tǒng)，噪聲大小直接影響其圖像質(zhì)量。為了抑制噪聲對(duì)重建圖像的影響，通過(guò)理論證明[4]，可以用兩個(gè)互補(bǔ)的編碼孔經(jīng)重建所得圖像相減來(lái)減小噪聲對(duì)重建圖像的影響。本文基于該思想對(duì)編碼板成像中的MLEM算法進(jìn)行改進(jìn)，提出MLEM互補(bǔ)算法。該算法在傳統(tǒng)MLEM算法公式中加入編碼板修正因子，在迭代算法中進(jìn)行互補(bǔ)編碼板相減計(jì)算。MLEM互補(bǔ)算法中最佳修正因子值的選取與編碼板線性衰減系數(shù)與厚度的乘積有關(guān)。通過(guò)蒙特卡羅模擬給出算法最佳修正因子的t擬合公式，并進(jìn)行實(shí)驗(yàn)驗(yàn)證。實(shí)驗(yàn)結(jié)果證明，與傳統(tǒng)MLEM算法相比，MLEM互補(bǔ)算法在成像時(shí)可以更有效地抑制噪聲。采用模擬所得-擬合公式，可以給出不同情況下MLEM互補(bǔ)算法最佳值。

1 算法簡(jiǎn)介

1.1 MLEM算法介紹

MLEM算法是基于泊松統(tǒng)計(jì)模型的最大貝葉斯后驗(yàn)概率的圖像復(fù)原方法，該算法由Shepp等[5]于1982年提出。MLEM算法主要分為兩步：1) 計(jì)算完全投影數(shù)據(jù)的似然函數(shù)在測(cè)定數(shù)據(jù)及當(dāng)前參數(shù)估計(jì)值下的期望；2) 求出使期望最大化的參數(shù)估計(jì)值。

在γ輻射編碼成像中MLEM算法[6?8]表示為：

1.2 基于互補(bǔ)碼板的MLEM算法

基于互補(bǔ)編碼板相減去除噪聲的思想，在進(jìn)行MLEM迭代時(shí)，對(duì)迭代公式(1)中的編碼函數(shù)(,) 進(jìn)行修正。令(,)=′(,)+h″(,)，其中：′(,)為正常編碼函數(shù)；″(,)為′(,)逆時(shí)針旋轉(zhuǎn)90°所得編碼函數(shù)；為算法修正因子，取值范圍為(?0.8,0)。將(,)帶入迭代公式，可得：

在計(jì)算第次放射源強(qiáng)度分布的迭代估計(jì)值和編碼函數(shù)所計(jì)算得出的投影估計(jì)值時(shí)使用互補(bǔ)編碼板相減去噪的方式對(duì)噪聲進(jìn)行抑制。由式(2)可知，由于在計(jì)算時(shí)對(duì)進(jìn)行了噪聲抑制，使得實(shí)際測(cè)定值與修正理論投影值的比值對(duì)重建圖像的修正效果更好。將該值與修正后的(,)做相關(guān)運(yùn)算，使得迭代算法再次對(duì)噪聲進(jìn)行抑制，進(jìn)一步減少噪聲對(duì)重建圖像的影響。

2 模擬及實(shí)驗(yàn)結(jié)果分析

2.1 模擬條件及模擬結(jié)果分析

使用蒙特卡羅模擬軟件MCNP5 (Monte Carlo N Particle Transport Code 5)對(duì)γ輻射編碼板成像過(guò)程進(jìn)行模擬，得到實(shí)驗(yàn)所需編碼板投影數(shù)據(jù)。建立的γ輻射編碼板成像模型示意圖如圖1所示。編碼板以11×11擴(kuò)展模式的MURA編碼板為基礎(chǔ)，幾何尺寸為50mm×50mm×10mm，采用方孔模式，孔徑2mm，材料選用鉛；點(diǎn)源位于視野中心軸線，距編碼板500mm處，能量為662keV；探測(cè)器材料選用鍺酸鉍(Bi2O3-GeO2, BGO)閃爍晶體，幾何尺寸為50mm×50mm×15mm，距編碼板20mm，計(jì)數(shù)時(shí)，將BGO晶體模型按100×100×1進(jìn)行網(wǎng)格劃分，并采用*F8計(jì)數(shù)卡對(duì)每個(gè)網(wǎng)格內(nèi)沉積能量進(jìn)行記錄。通過(guò)蒙特卡羅模擬所得理想點(diǎn)源圖像及其投影圖像如圖2所示。

圖1 編碼板模型示意圖

圖2 理想點(diǎn)源(a)和點(diǎn)源投影(b)圖像

采用蒙特卡羅模擬所得編碼板投影數(shù)據(jù)，分別通過(guò)MLEM算法和MLEM互補(bǔ)算法對(duì)圖像進(jìn)行重建。圖3為MLEM算法與MLEM互補(bǔ)算法(=?0.6)在迭代次數(shù)=5時(shí)的重建圖像。由圖3可知，兩種算法都可對(duì)模擬編碼圖像進(jìn)行有效重建。采用MLEM互補(bǔ)算法對(duì)模擬圖像進(jìn)行重建，其對(duì)噪聲的抑制效果更好，重建圖像更接近模擬的點(diǎn)源圖像。通過(guò)直接觀察法可知，MLEM互補(bǔ)算法的成像效果要優(yōu)于傳統(tǒng)的MLEM算法。

對(duì)重建圖像進(jìn)行評(píng)估時(shí)，主觀評(píng)價(jià)往往很難判斷算法的優(yōu)劣。為進(jìn)一步驗(yàn)證算法的有效性，在對(duì)重建圖像的質(zhì)量進(jìn)行評(píng)價(jià)時(shí)，采用歸一化均方誤差(Normalized Root Mean Square Error, NRMSE)以及相關(guān)系數(shù)(Correlation Coefficient, CC)這兩種定量方法對(duì)重建圖像的質(zhì)量進(jìn)行評(píng)價(jià)。用歸一化均方誤差衡量算法的接近程度，相關(guān)系數(shù)衡量重建圖像與原圖像相似程度[9]。NRMSE與CC評(píng)價(jià)標(biāo)準(zhǔn)的值和定義分別為：

表1 MLEM算法和MLEM互補(bǔ)算法NRMSE與CC值

2.2 實(shí)驗(yàn)條件及實(shí)驗(yàn)結(jié)果分析

為驗(yàn)證模擬結(jié)果的正確性，通過(guò)實(shí)驗(yàn)對(duì)算法的有效性進(jìn)行驗(yàn)證。實(shí)驗(yàn)中編碼板以11×11擴(kuò)展的MURA編碼模式為基礎(chǔ)，編碼板幾何尺寸分別為50mm×50mm×10mm、50mm×50mm×20mm、50mm×50mm×30mm，孔型為方孔，孔徑2mm，材料選用鉛；探測(cè)器由50mm×50mm×10mm的BGO陣列晶體，耦合濱松H8500位置靈敏光電倍增管組成，探測(cè)器距編碼板20mm。點(diǎn)源采用強(qiáng)度3.7×104Bq的137Cs源，位于視野中心軸線，距編碼板500mm處。

采用實(shí)驗(yàn)所得編碼板投影數(shù)據(jù)，分別通過(guò)MLEM算法和MLEM互補(bǔ)算法對(duì)圖像進(jìn)行重建。圖4為MLEM算法和MLEM互補(bǔ)算法在迭代次數(shù)=5時(shí)對(duì)三種厚度編碼板實(shí)測(cè)投影重建圖像。由圖4(a)可知，隨著編碼板厚度的增加，MLEM算法重建圖像的背景噪聲越來(lái)越小，圖像質(zhì)量逐步提高。由圖4(b)可知，采用MLEM互補(bǔ)算法進(jìn)行圖像重建時(shí)，在編碼板厚度相同的情況下，MLEM互補(bǔ)算法的重建結(jié)果比MLEM算法的重建結(jié)果對(duì)噪聲的抑制更明顯。

圖4 MLEM算法(a)與MLEM互補(bǔ)算法(b)實(shí)驗(yàn)數(shù)據(jù)重建圖像

表2為MLEM算法和MLEM互補(bǔ)算法在迭代次數(shù)=5時(shí)對(duì)三種厚度編碼板實(shí)測(cè)投影重建圖像的NRMSE和CC值。由表2可知，在編碼板厚度相同的情況下，MLEM互補(bǔ)算法的NRMSE和CC值皆優(yōu)于傳統(tǒng)的MLEM算法，實(shí)驗(yàn)結(jié)果與模擬結(jié)果相一致。與MLEM算法相比，采用MLEM互補(bǔ)算法對(duì)編碼圖像進(jìn)行重建，可以更有效地抑制噪聲，提高重建圖像的質(zhì)量。

表2 MLEM算法和MLEM互補(bǔ)算法的NRMSE、CC值

2.3值對(duì)MLEM互補(bǔ)算法的影響

通過(guò)MLEM互補(bǔ)算法對(duì)上述實(shí)驗(yàn)及模擬投影數(shù)據(jù)進(jìn)行重建時(shí)，=?0.6。為研究不同值對(duì)MLEM互補(bǔ)算法重建效果的影響，選取不同值對(duì)實(shí)驗(yàn)所得厚度=10mm、=20mm、=30mm三種厚度編碼板的實(shí)測(cè)投影值進(jìn)行互補(bǔ)MLEM算法重建。圖5為三種厚度下互補(bǔ)MLEM算法的NRMSE和CC值隨修正因子值的變化。

圖5 三種厚度鉛制編碼板NRMSE (a)和CC (b)值隨a變化

由圖5可知，隨著值在(?0.8,0)區(qū)間內(nèi)不斷增大，重建圖像的質(zhì)量先隨之提高，當(dāng)達(dá)到某一特定值后，重建圖像的質(zhì)量開(kāi)始下降。對(duì)于MLEM算法存在最佳值，使得改進(jìn)后MLEM算法效果最優(yōu)。編碼板厚度為1 cm、2 cm、3 cm時(shí)，最佳值分別為?0.7、?0.61、?0.55。可見(jiàn)隨著編碼板準(zhǔn)直器厚度的增加，MLEM互補(bǔ)算法修正因子的最佳值逐漸減小，最佳值的選取與物質(zhì)阻止本領(lǐng)有關(guān)。由γ射線在物質(zhì)中的衰減規(guī)律，采用上述編碼板模型，模擬不同t值下，算法最佳值的變化。圖6為最佳值隨t的變化。由圖6可知，最佳值隨t值一同增加，當(dāng)t值達(dá)到某一特定值時(shí)，值的變化趨于平緩。計(jì)算所得t擬合公式為：

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

圖6 最佳a(bǔ)值隨mt值變化擬合曲線

表3為三種厚度下鉛制編碼板實(shí)測(cè)投影數(shù)據(jù)所得最佳值與模擬數(shù)據(jù)擬合公式計(jì)算所得最佳值。

表3 三種厚度鉛制編碼板實(shí)驗(yàn)及模擬數(shù)據(jù)最佳a(bǔ)值

由表3可知，由于噪聲影響，實(shí)驗(yàn)數(shù)據(jù)所得值比模擬數(shù)據(jù)所得的要大。無(wú)法通過(guò)模擬數(shù)據(jù)所得t擬合公式直接給出實(shí)測(cè)情況下MLEM互補(bǔ)算法最佳值。在使用t擬合公式給出實(shí)際成像系統(tǒng)所需最佳值時(shí)，需對(duì)模擬所得最佳值乘以相應(yīng)比例系數(shù)。

記實(shí)驗(yàn)值與模擬值的比值為比例系數(shù)，根據(jù)t擬合公式變化趨勢(shì)，對(duì)隨t值變化進(jìn)行擬合。圖7為比例系數(shù)隨t值的變化。

圖7 比例系數(shù)N隨最佳mt值變化擬合曲線

計(jì)算所得t擬合公式為：

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

在求MLEM互補(bǔ)算法的最佳值時(shí)，可以通過(guò)模擬所得t擬合公式計(jì)算結(jié)果乘以得到。

為驗(yàn)證上述結(jié)論，采用1cm鎢板進(jìn)行實(shí)驗(yàn)，圖8為鎢板的NRMSE和CC值隨值的變化。由圖8可知，此時(shí)MLEM互補(bǔ)算法最佳值為0.65。在入射γ射線能量為662keV時(shí)，鎢的線性衰減系數(shù)=1.84cm?1[10]；將=1cm帶入擬合公式(6)，得出模擬情況下算法最佳值為0.41，由擬合公式(7)，此時(shí)為1.6。最終用于1cm鎢制編碼板的最佳=1.6×0.41=0.65。與實(shí)驗(yàn)測(cè)得數(shù)據(jù)所得最佳值相一致。在實(shí)際應(yīng)用中，可以用上述方法確定MLEM互補(bǔ)算法最佳值。

圖8 1cm厚鎢制編碼板NRMSE (a)和CC (b)值隨a變化

3 結(jié)語(yǔ)

實(shí)驗(yàn)證明，將互補(bǔ)編碼板相減去噪思想應(yīng)用于MLEM算法中是可行的。通過(guò)對(duì)蒙特卡羅模擬所得投影數(shù)據(jù)及實(shí)測(cè)所得投影數(shù)據(jù)的重建結(jié)果可知，相比于MLEM算法，在編碼板成像中，MLEM互補(bǔ)算法可以有效抑制噪聲的影響，提高圖像重建質(zhì)量；MLEM互補(bǔ)算法中的最佳值的選取與t有關(guān)；通過(guò)模擬得出t擬合公式所得模擬最佳值與t擬合公式所得的乘積，可以得到實(shí)測(cè)數(shù)據(jù)的最佳值，模擬結(jié)果與實(shí)驗(yàn)結(jié)果相一致。采用該方法可以確定MLEM互補(bǔ)算法最佳值；在編碼板成像系統(tǒng)中，采用蒙特卡羅模擬方法對(duì)其進(jìn)行分析具有一定的可行性。

1 趙翠蘭, 陳立宏, 李勇平. MURA編碼輻射成像系統(tǒng)的解碼方法[J]. 核技術(shù), 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ìn)算法的研究和比較[J]. 計(jì)算機(jī)工程與設(shè)計(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ī)學(xué)儀器與方法[M]. 哈爾濱: 哈爾濱工程大學(xué)出版社, 2010.

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

4 Barrett H H, Swindell W. 放射成像[M]. 莊天戈, 周頌凱, 譯. 北京: 科學(xué)出版社, 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]. 原子能科學(xué)技術(shù), 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 洪俊杰. γ相機(jī)中編碼孔經(jīng)準(zhǔn)直器設(shè)計(jì)及數(shù)字圖像重建[D]. 湖北: 華中科技大學(xué), 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 何佳偉, 劉東升, 桂志國(guó). 可變有序子集PML算法在PET中的應(yīng)用[J]. 中北大學(xué)學(xué)報(bào)(自然科學(xué)版), 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 方杰. 輻射防護(hù)導(dǎo)論[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è)于上海交通大學(xué)，現(xiàn)為碩士研究生，輻射防護(hù)與環(huán)境保護(hù)專業(yè)

2016-11-20，

2016-12-14