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

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

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

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

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

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

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

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

本文將通過(guò)引入加速因子h來(lái)調(diào)制校正因子的步長(zhǎng)，以改善OSEM算法的收斂特性，進(jìn)而提高XFCT重建速度。

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

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

式中，Ck+1(j)表示測(cè)量得到的投影值經(jīng)過(guò)第k+1次迭代后第i個(gè)像素的近似值，系數(shù)K(i,j)表示像素i對(duì)投影強(qiáng)度yj的貢獻(xiàn)因子。分母為第k次迭代后的投影估計(jì)值：

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

對(duì)式(3)進(jìn)行泰勒一階展開(kāi)，忽略高次項(xiàng)有：

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

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

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

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

Sl(l=1, 2, …, L)表示有L個(gè)子集數(shù)，當(dāng)L=1時(shí)為MLEM算法。

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

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

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

2 數(shù)字模擬

2.1 數(shù)字模型

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

2.2 結(jié)果與討論

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

圖1 模擬實(shí)驗(yàn)的重建結(jié)果：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算法的重建效果，視覺(jué)比較和歸一化的均方根誤差方法(Normalized root mean square error, NrMSE)被引入作為評(píng)價(jià)XFCT圖像重建質(zhì)量的標(biāo)準(zhǔn)，NrMSE公式定義如下：

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

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

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

圖2 AOSEM算法子集數(shù)分別為5、10、15、20時(shí)重建結(jié)果與標(biāo)準(zhǔn)圖像的偏差（迭代次數(shù)為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實(shí)驗(yàn)表明，加速因子h具有一定的上限。圖3探究了具有不同加速因子的AOSEM算法在不同迭代次數(shù)下重建圖像的收斂速度。AOSEM算法的子集數(shù)選擇為10，加速因子h分別取1、1.5、2、2.2、2.5、3，當(dāng)h=1時(shí)為傳統(tǒng)的OSEM算法，重建結(jié)果顯示h＞1的AOSEM算法收斂速度要明顯快于OSEM算法，較少迭代次數(shù)時(shí)加速因子越大重建圖像收斂速度越快，這些曲線在多次迭代后趨于同一值。但當(dāng)h＞3后重建圖像會(huì)出現(xiàn)嚴(yán)重的離散化。

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

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

圖4 兩種算法不同迭代次數(shù)下剖線圖比較(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 實(shí)驗(yàn)研究

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

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

圖5 測(cè)試樣品驗(yàn)證(a) 測(cè)試樣品實(shí)物圖，兩個(gè)黑色箭頭指向的管中放有CdCl2溶液，(b) Cd元素投影數(shù)據(jù)正弦圖，(c)、(d) OSEM算法迭代4次和16次，(e)、(f) AOSEM算法h=2時(shí)迭代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 結(jié)語(yǔ)

在XFCT重建中，雖然OSEM算法和MLEM算法比較具有較快的重建速度，但收斂速度依然較慢。本文提出了一種快速重建熒光CT的方法——AOSEM，通過(guò)引入加速因子h來(lái)加快OSEM算法的收斂速度，當(dāng)h＞1時(shí)相比同階OSEM算法具有更好的重建圖像質(zhì)量。加速因子h的取值上限也是有所限制的，否則會(huì)因?yàn)榧铀龠^(guò)快導(dǎo)致圖像的離散化。計(jì)算機(jī)模擬和測(cè)試樣品實(shí)驗(yàn)結(jié)果表明，獲得與OSEM算法同等質(zhì)量的XFCT重建圖像，AOSEM算法所需的迭代次數(shù)僅是傳統(tǒng)OSEM算法的一半。通過(guò)比較兩種算法重建圖像的歸一化均方根誤差可以得到相似的結(jié)論。本文的研究結(jié)果表明，具有加速因子的OSEM重構(gòu)算法可有效提高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融合顯像對(duì)肺癌單發(fā)骨轉(zhuǎn)移瘤的診斷價(jià)值[J]. 核技術(shù), 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].中國(guó)圖象圖形學(xué)報(bào): 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ī)學(xué)應(yīng)用研究[D]. 北京: 中國(guó)科學(xué)院大學(xué), 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)、國(guó)家自然科學(xué)基金聯(lián)合基金重點(diǎn)項(xiàng)

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

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

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

2015-01-05，

2015-03-22