具有加速因子的OSEM重建算法用于X射線熒光CT研究

孫鵬飛 鄧 彪 楊 群 杜國浩 佟亞軍 肖體喬

具有加速因子的OSEM重建算法用于X射線熒光CT研究

孫鵬飛1,2鄧 彪1楊 群1杜國浩1佟亞軍1肖體喬1,2

1（中國科學院上海應用物理研究所 嘉定園區(qū) 上海 201800）
2（中國科學院大學 北京 100049）

有序子集最大期望值算法(Ordered Subsets Expectation Maximization, OSEM)具有較高的圖像重建質量和較短的計算時間，已經被應用于內源CT（如SPECT、PET、同步輻射X射線熒光CT）的圖像重建中。本文提出了一種具有加速因子的OSEM算法應用于X射線熒光CT的圖像重建，通過引入加速因子h來調制校正因子的步長加快OSEM算法的收斂速度，研究了不同加速因子和不同子集數的AOSEM算法對重建圖像質量的影響。計算機模擬及實驗結果表明，在獲得同等質量重建圖像的同時，具有加速因子的OSEM算法的重建速度是常規(guī)OSEM的兩倍。

同步輻射，X射線熒光CT，圖像重建質量，有序子集最大期望值算法

X射線熒光CT (X-ray fluorescent computed tomography, XFCT)通過探測樣品內某元素受激發(fā)射的特征X射線結合CT重建算法重建元素在樣品內部的空間分布[1]，具有無損、高分辨等優(yōu)點。隨著第三代同步輻射光源和探測技術的發(fā)展，XFCT空間分辨率可達納米量級，在材料、生物醫(yī)學等領域的研究得到了廣泛應用[2?5]。

目前CT重構常用的方法包括濾波反投影(Filter Back Projection, FBP)和迭代算法。迭代重建算法因其在有限的投影數和噪聲影響情況下，能提供比濾波反投影更高質量的重建圖像，被廣泛應用于投影角度較少的內源CT (SPECT、XFCT)的圖像重建中。1982年Shepp等[6]提出了最大似然期望值法(Maximum Likelihood Expectation Maximization, MLEM)，該方法易于實現(xiàn)并可提供很好的重建效果，但收斂速度較慢。為加快收斂速度，Hudson等[7]將有序子集(Ordered Subsets, OS)塊迭代的思想應用到MLEM算法上，由此得到的OSEM算法，其每個子集遍歷投影數據時都要對圖像更新一次，有效地加快了重建速度，被應用于SPECT等內源CT圖像重建中[8?9]。上海光源X射線成像組已將OSEM算法應用于XFCT的圖像重建中，但圖像重建的速度依舊較慢[10?11]，限制了OSEM算法在XFCT中的進一步推廣和應用。

本文將通過引入加速因子h來調制校正因子的步長，以改善OSEM算法的收斂特性，進而提高XFCT重建速度。

1 具有加速因子的OSEM算法原理

X射線熒光CT的MLEM重建算法如下[6]：

式中，Ck+1(j)表示測量得到的投影值經過第k+1次迭代后第i個像素的近似值，系數K(i,j)表示像素i對投影強度yj的貢獻因子。分母為第k次迭代后的投影估計值：

加速因子h作為指數引入到式(1)：

對式(3)進行泰勒一階展開，忽略高次項有：

式(4)給出了具有加速因子MLEM算法(Accelerated Expectation Maximization, AEM)。當h=1時為傳統(tǒng)MLEM算法；當h＞1時可以起到增大校正因子進而增大重建圖像更新步長的作用[12?13]。然而，加速因子h的大小受到圖像更新后像素值正負的影響。根據式(1)，MLEM算法重建后的像素可以避免出現(xiàn)負值，但引入加速因子h以后，當h＞1時不能保證的非負性。因此，需要在每次迭代之后加入約束條件[14]：

對于MLEM式(1)，具有自動保持每次迭代總計數守恒的特性。然而，引入加速因子h以后，AEM重建算法無法保持近似解總計數守恒，會引起圖像的離散化和變形。為保證MLEM算法這一特性，必須在每次迭代之后加入約束條件：

加入約束條件的AEM算法，可以通過有效的增大因子h來加快MLEM算法的收斂速度。根據式(6)，下一次迭代重建結合了本次迭代計算所得的投影值。

Hudson等[15]將投影值等分成若干個有序集合，稱為有序子集。一次迭代過程中每個子集都會更新圖像一次，所以一次OSEM迭代相當于MLEM算法對圖像更新了S次，S為有序子集數。XFCT的OSEM重建公式如下：

Sl(l=1, 2, …, L)表示有L個子集數，當L=1時為MLEM算法。

具有加速因子的OSEM算法AOSEM公式可以寫成：

同AEM，AOSEM算法同樣需要在每個子迭代后加入約束條件：

式(8)就是AOSEM算法XFCT迭代重建公式，當加速因子h=1時，為傳統(tǒng)的OSEM算法公式，h>1時可以增大校正因子的步長，加快收斂速度。但根據式(8)，如果h取值太大，每次子迭代更新的像素大小都會成為負值被約束為零，所以h的取值具有一定的上限，否則會導致圖像的離散化。

2 數字模擬

2.1 數字模型

XFCT模擬實驗采用的數字模型由軟組織材料和元素Zn組成[16]。Zn元素亮區(qū)域的質量分數為0.1%，暗區(qū)域為0.05%，其他為軟組織材料。數據規(guī)模為100×100，每個像素大小為5 μm×5 μm。投影數據的獲?。盒D角度采樣范圍180°，采樣間隔為3°；固定投影角下，進行100次平動（步長為5μm），投影數據中加入了泊松噪聲，投影總計數約為6.0×108。分別采用OSEM算法、AOSEM算法進行迭代重建，由于Zn元素熒光能量較低，樣品自身對其吸收效應較強，所以在重建過程中加入了吸收修正，初始化每個圖像的像素大小為0.00001。

2.2 結果與討論

模擬得到的投影數據采用AOSEM算法和傳統(tǒng)的OSEM算法進行不同次數的迭代重建，結果如圖1所示。上排為傳統(tǒng)的OSEM算法（迭代次數分別為2、4、8、16）；下排為加速因子h=2時的AOSEM算法的重建結果（迭代次數分別為1、2、4、8）。

圖1 模擬實驗的重建結果：OSEM算法（上排）2、4、8、16次迭代；AOSEM算法（下排）1、2、4、8次迭代Fig.1 Reconstructed result of the numerical phantom: OSEM (top row) with 2, 4, 8 and 16 iterations and AOSEM (bottom row) using h=2 with 1, 2, 4 and 8 iterations.

為研究具有加速因子的OSEM算法的重建效果，視覺比較和歸一化的均方根誤差方法(Normalized root mean square error, NrMSE)被引入作為評價XFCT圖像重建質量的標準，NrMSE公式定義如下：

通過視覺觀察比較，相同迭代次數下OSEM算法的重建結果要差于本文提出的帶有加速因子AOSEM算法。表1給出了兩種算法不同迭代次數下的NrMSE值，我們發(fā)現(xiàn)OSEM算法迭代2、4、8、16的重建圖像的NrMSE值幾乎與OSEM算法迭代1、2、4、8相對應。而在相同迭代次數下，AOSEM算法的NrMSE值要小于OSEM算法。因此，我們認為AOSEM算法可以提供與傳統(tǒng)OSEM算法相同質量的重建圖像，但是所需的迭代次數僅是OSEM算法的一半。

表1 兩種算法重建圖像的NrMSE值Table 1 NrMSE of the reconstructed images using two algorithms.

在迭代法重建中，介入有序子集是一種加快收斂速度很流行的方法。然而當子集數超過一定值后，重建圖像質量反而會隨著子集數的增加而降低。這是因為子集數目越多，每個子集內包含的投影數越少，導致重建圖像引入過多噪聲[15]。所以選取最佳子集數對重建圖像結果影響甚大。圖2描述了子集數分別為5、10、15、20時，AOSEM算法重建圖像與標準圖像偏差隨迭代次數的變化曲線。相同迭代次數下，隨著子集數的增加重建圖像的偏差減小，圖像收斂速度加快。但當子集數增加到15后，可以看出子集數為15和20時，隨著迭代次數的增加NrMSE值變化相當接近，表明圖像的收斂速度明顯減緩，且子集數繼續(xù)增加反而會降低圖像質量。

圖2 AOSEM算法子集數分別為5、10、15、20時重建結果與標準圖像的偏差（迭代次數為2，加速因子h=2）Fig.2 NrMSE vs. iterations for AOSEM with various subset numbers of 5, 10, 15 and 20 (Iteration=2, h=2).

模擬XFCT實驗表明，加速因子h具有一定的上限。圖3探究了具有不同加速因子的AOSEM算法在不同迭代次數下重建圖像的收斂速度。AOSEM算法的子集數選擇為10，加速因子h分別取1、1.5、2、2.2、2.5、3，當h=1時為傳統(tǒng)的OSEM算法，重建結果顯示h＞1的AOSEM算法收斂速度要明顯快于OSEM算法，較少迭代次數時加速因子越大重建圖像收斂速度越快，這些曲線在多次迭代后趨于同一值。但當h＞3后重建圖像會出現(xiàn)嚴重的離散化。

圖3 具有不同加速因子h的AOSEM算法重建結果Fig.3 NrMSE vs. iterations for AOSEM with various accelerated factor (h=1, 1.5, 2, 2.2, 2.5 and 3).

圖3中，迭代次數為1時，具有較大加速因子h=2.5、3的AOSEM算法重建結果偏差要大于OSEM算法。這可能是由于在初始迭代，圖像的重建值遠離收斂值，AOSEM算法導致了這一偏差的放大。解決這一問題可以考慮引入一個合適的初始值，比如將重建速度較快的FBP算法作為AOSEM算法的初始迭代值。當迭代次數增加到5次后，收斂速度并不是隨著加速因子而增加，出現(xiàn)了交叉值。取重建圖像上圓形線路徑上的像素與標準模型進行比較。以圓形線右箭頭的起點作為坐標原點，可以得到兩種算法重建圖像的ROI剖線圖，縱坐標為歸一化濃度分布，橫坐標為圓形線與起點之間的路程。圖4(a)為傳統(tǒng)OSEM算法、AOSEM算法分別迭代4次的結果與模擬模型的偏差；圖4(b)給出了傳統(tǒng)OSEM算法迭代8次、AOSEM算法迭代4次與模擬模型的偏差。在同樣4次迭代后，AOSEM算法重建結果更接近于模型真實值，尤其在抑制背景噪聲和提高ROI邊緣像素的更新速度表現(xiàn)更好。而當OSEM算法迭代8次后，與AOSEM算法迭代4次圖像的收斂效果相同。因此，獲得同等質量的重建結果，AOSEM算法可以節(jié)約一半的迭代次數。

圖4 兩種算法不同迭代次數下剖線圖比較(a) OSEM和AOSEM都迭代4次，(b) OSEM迭代8次，AOSEM迭代4次Fig.4 Comparisons of the profiles between the OSEM and AOSEM with various iterations. (a) OSEM and AOSEM both with 4 iterations, (b) OSEM with 8 iterations while AOSEM with 4 iterations

3 實驗研究

模擬實驗結果證明了AOSEM算法具有較好的重建質量和收斂特性。通過測試樣品進一步驗證本文提出的算法。測試樣品為一個直徑10 mm的聚甲基丙烯酸酯圓柱模型。內部含有4個直徑為3 mm的孔洞，其中對稱排列的兩個孔洞被添加了濃度為10?6mol·mL?1的CdCl2溶液[17]，另兩個為空。實驗在上海光源BL13W1束線站進行，入射光能量為32keV的單色光，照射到樣品處的光斑通過狹縫限定為200 μm×200 μm。在數據采集過程中，樣品每次轉動3°，完成180°范圍內的掃描。在固定投影角下，樣品平動步長為200 μm，單幅熒光譜的采集時間為3 s。

圖5給出了兩種重建算法不同迭代次數下Cd的分布。沿著圖5(c)中標記的白色線取兩種算法分別迭代4次的剖線圖，結果如圖6所示。由于入射光的能量(32 keV)和Cd元素Kα線能量（約為23 keV）較高，吸收效應影響較弱，從圖6中可以看出重建值比較接近真實值。通過比較圖5(c)和(e)，AOSEM算法明顯具有更高的圖像重建質量。另外，圖6顯示在像素值為零沒有Cd元素分布的地方，AOSEM算法重建的像素值更接近于零，含Cd元素區(qū)域邊界清晰。比較圖5(d)和(f)，常規(guī)OSEM算法16次迭代的重建結果與AOSEM算法迭代8次獲得的圖像幾乎相近。因此，可以得到與模擬實驗相同的結論，具有加速因子的AOSEM算法重建圖像質量更好，收斂速度更快，大約可以節(jié)省一半的迭代次數。

圖5 測試樣品驗證(a) 測試樣品實物圖，兩個黑色箭頭指向的管中放有CdCl2溶液，(b) Cd元素投影數據正弦圖，(c)、(d) OSEM算法迭代4次和16次，(e)、(f) AOSEM算法h=2時迭代4次和8次Fig.5 Test sample for validation. (a) Test sample, (b) The sinogram of Cd projection, (c), (d) Reconstructed image of Cd distribution using OSEM with 4 and 16 iterations respectively, (e), (f) Reconstructed image of Cd distribution using AOSEM with 4 and 8 iterations respectively

圖6 Cd元素重建像素值比較Fig.6 Comparisons of the profiles along white line between the OSEM and AOSEM algorithm.

4 結語

在XFCT重建中，雖然OSEM算法和MLEM算法比較具有較快的重建速度，但收斂速度依然較慢。本文提出了一種快速重建熒光CT的方法——AOSEM，通過引入加速因子h來加快OSEM算法的收斂速度，當h＞1時相比同階OSEM算法具有更好的重建圖像質量。加速因子h的取值上限也是有所限制的，否則會因為加速過快導致圖像的離散化。計算機模擬和測試樣品實驗結果表明，獲得與OSEM算法同等質量的XFCT重建圖像，AOSEM算法所需的迭代次數僅是傳統(tǒng)OSEM算法的一半。通過比較兩種算法重建圖像的歸一化均方根誤差可以得到相似的結論。本文的研究結果表明，具有加速因子的OSEM重構算法可有效提高X射線熒光CT的圖像重建速度。

1 Simionovici A, Chukalina M, Günzler F, et al. X-ray microtome by fluorescence tomography[J]. Nuclear Instruments and Methods in Physics Research Section A, 2001, 467(2): 889?892

2 Boisseau P. Determination of three dimensional trace element distributions by the use of monochromatic X-ray microbeams[D]. Massachusetts Institute of Technology, Department of Physics, 1986

3 DeJonge M D, Holzner C, Baines S B, et al. Quantitative 3D elemental microtomography of Cyclotella meneghiniana at 400-nm resolution[J]. Proceedings of the National Academy of Sciences, 2010, 107(36): 15676?15680

4 Takeda T, Wu J, Huo Q, et al. X-ray fluorescent CT imaging of cerebral uptake of stable-iodine perfusion agent iodoamphetamine analog IMP in mice[J]. Journal of Synchrotron Radiation, 2008, 16(1): 57?62

5 Kim S A, Punshon T, Lanzirotti A, et al. Localization of iron in Arabidopsis seed requires the vacuolar membrane transporter VIT1[J]. Science, 2006, 314(5803): 1295?1298

6 Shepp L A, Vardi Y. Maximum likelihood reconstruction for emission tomography[J]. IEEE Transactions on Medical Imaging, 1982, 1(2): 113?122

7 Hudson H M, Larkin R S. Accelerated image reconstruction using ordered subsets of projection data[J]. IEEE Transactions on Medical Imaging, 1994, 13(4): 601?609

8 Dai T, Ma T, Liu H, et al. A high-resolution small animal SPECT system developed at Tsinghua[J]. Nuclear Science and Techniques, 2011, 22: 344?348

9 常城, 謝文暉, 雷貝, 等. SPECT/CT融合顯像對肺癌單發(fā)骨轉移瘤的診斷價值[J]. 核技術, 2013, 36(9): 090301

CHANG Cheng, XIE Wenhui, LEI Bei, et al. SPECT/CT fusion imaging for differential diagnosis of bone solitary metastases in patients with lung cancer[J]. Nuclear Techniques, 2013, 36(9): 090301

10 Yang Q, Deng B, Du G, et al. X-ray fluorescence computed tomography with absorption correction for biomedical samples[J]. X-ray Spectrometry, 2014, 43(5): 278?285

11 Yang Q, Deng B, Lyu W, et al. Fast and accurate X-ray fluorescence computed tomography imaging with the ordered-subsets expectation maximization algorithm[J]. Journal of Synchrotron Radiation, 2012, 19(2): 210?215

12 Zhao X, Hu J, Yang T, et al. An accelerated RAMLA reconstruction algorithm for X-ray cone-beam CT[J]. Insight-Non-Destructive Testing and Condition Monitoring, 2013, 55(5): 237?242

13 Hsiao T, Rangarajan A, Khurd P, et al. An accelerated convergent ordered subsets algorithm for emission tomography[J]. Physics in Medicine and Biology, 2004, 49(11): 2145?2156

14 劉力, 吳朝霞. 具有超松弛因子的OSEM重建算法[J].中國圖象圖形學報: A輯, 2002, 7(8): 814?817

LIU Li, WU Zhaoxia. OSEM reconstruction with overrelaxation[J]. Journal of Image and Graphics: A, 2002, 7(8): 814?817

15 Hudson H, Larkin R. Accelerated image reconstruction using ordered subsets of projection data[J]. IEEE Transactions on Medical Imaging, 2002, 13(4): 601?609

16 楊群. 快速低劑量X射線熒光CT及其生物醫(yī)學應用研究[D]. 北京: 中國科學院大學, 2012

YANG Qun. Fast, low dose X-ray fluorescence CT for the biomedical applications[D]. Beijing: University of Chinese Academy of Sciences, 2012

17 Deng B, Yang Q, Xie H L, et al. First X-ray fluorescence CT experimental results at the SSRF X-ray imaging beamline[J]. Chinese Physics C, 2011, 35(4): 402?404

CLC TL99

An accelerated OSEM reconstruction algorithm using an accelerating factor for X-ray fluorescence tomography

SUN Pengfei1,2DENG Biao1YANG Qun1DU Guohao1TONG Yajun1XIAO Tiqiao1,2
1(Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Jiading Campus, Shanghai 201800, China)
2(University of the Chinese Academy of Sciences, Beijing 100049, China)

Background: The Ordered Subsets Expectation Maximization (OSEM) reconstruction algorithm has been widely applied in emission CT reconstruction, such as SPECT, PET and synchrotron radiation X-ray fluorescence CT (XFCT). The quality of reconstructed image is better than other analytical methods such as filtered-back projection algorithm. However, the convergent rate of the OSEM is slow. Purpose: In order to speed up the convergent rate of OSEM, we want to investigate an improved OSEM algorithm. Methods: In this paper, we present an accelerated OSEM algorithm (AOSEM) by increasing the step size of the correction item and show its convergence characteristics with various subsets and accelerated factors. Results: The AOSEM algorithm is proposed and applied in XFCT image reconstruction. Comparing the two reconstruction algorithms, both the simulations and experimental results showed that AOSEM reached the same image quality as in OSEM but only about half the number of iterations when an accelerated factor was used. Conclusion: AOSEM can further speed up the convergence of OSEM when the accelerated factor h＞1.

Synchrotron radiation, XFCT, Construction image quality, OSEM algorithm

TL99

10.11889/j.0253-3219.2015.hjs.38.060201

No.11275257、No.11105213、No.31100680、No.51274054、No.11375257、No.31300480)、國家自然科學基金聯(lián)合基金重點項

目(No.U1232205/A0802)、CAS-CSIRO合作研究項目(No.GJHZ1303)資助

孫鵬飛，男，1989年出生，2012年畢業(yè)于蘭州大學，粒子物理與原子核物理

鄧彪，E-mail: dengbiao@sinap.ac.cn；肖體喬，E-mail: tqxiao@sinap.ac.cn

2015-01-05，

2015-03-22