基于OSEM算法的X射線熒光CT重建算法研究

劉亞楠

（重慶工商職業(yè)學(xué)院，重慶 400052）

X射線熒光CT（X-ray Fluorescence Computed Tomography，XFCT）[1-2]，是將X射線熒光分析（X-ray Fluorescence analysis，XRF）[3]和X射線CT（X-ray Computed Tomography，X-CT）[4]技術(shù)相結(jié)合的產(chǎn)物，它是一種新近發(fā)展起來的無損檢測與分析手段，能以非侵入、無損的方式測量樣品中元素的分布和含量[5]。鑒于X射線熒光CT可對樣品中元素的分布和含量進(jìn)行無損分析，能夠彌補(bǔ)單一X射線熒光分析和X射線CT之不足，XFCT在生物醫(yī)學(xué)、植物學(xué)、藥學(xué)、地球科學(xué)等領(lǐng)域的應(yīng)用逐漸增多，成為該領(lǐng)域研究的熱點(diǎn)。本文以XFCT圖像重建為出發(fā)點(diǎn)，基于有序子集—期望最大化（Ordered-Subsets Expectation Maximization，OSEM）方法[6]，研究了基于OSEM的XFCT重建算法，并進(jìn)行相關(guān)仿真和驗證，獲得良好的效果。

1 X射線熒光CT原理

通常，X射線熒光CT是以平移—旋轉(zhuǎn)方式掃描獲得投影數(shù)據(jù)。建立如圖1所示的兩個直角坐標(biāo)系，xy坐標(biāo)系為固定于樣品上的旋轉(zhuǎn)坐標(biāo)系，uv坐標(biāo)系為實驗室坐標(biāo)系。在掃描成像中，樣品繞坐標(biāo)原點(diǎn)作逆時針旋轉(zhuǎn)。因此，兩個坐標(biāo)系的關(guān)系可表示為：

考慮其中一束X射線在樣品逆時針旋轉(zhuǎn)θ后的投影。

（1）假定入射X射線束強(qiáng)度為I0, 其經(jīng)過樣品()時被吸收，在到達(dá)Q點(diǎn)之前X射線束的強(qiáng)度可表示為：

其中，μI(x,y)為在入射X射線光子能量下的吸收系數(shù)分布，Q點(diǎn)的坐標(biāo)為(u,v)。

（2）若Q處某元素此時能夠被激發(fā)產(chǎn)生X射線熒光，則產(chǎn)生的熒光強(qiáng)度與f(θ,u,v)，熒光產(chǎn)額ω，光電吸收系數(shù)μph以及此時元素的濃度ρ有關(guān)，微元Δu受激發(fā)產(chǎn)生的熒光并被探測器所探測到的強(qiáng)度為：

圖1 筆束X射線熒光CT結(jié)構(gòu)

強(qiáng)度可表示為：

其中：

可見，Ii是與元素濃度ρ(u,v)有關(guān)的熒光CT一個投影，X射線熒光CT就是根據(jù)獲得的所有投影數(shù)據(jù)重建出元素分布ρ(u,v)的圖像。

2 OSEM重建算法

OSEM算法中，將熒光CT投影數(shù)據(jù)分為T個經(jīng)過排序的子集{S1,S2,ST}，即有序子集，對每個投影數(shù)據(jù)依次使用標(biāo)準(zhǔn)的EM算法來最大化似然函數(shù)，重建的結(jié)果作為下一個子集的初值。OSEM算法可簡寫為：

每次重建時使用一個子集內(nèi)的投影數(shù)據(jù)同時對各像素進(jìn)行校正，重建圖像更新一次，完成一次迭代（所有子集都對像素校正一次）重建圖像已經(jīng)更新了T次，從而大大降低重建時間。

3 實驗與討論

為驗證OSEM算法用于熒光CT重建的效果，本文采用如圖2(a)所示的數(shù)值模體。其中，子集個數(shù)為5，角度采樣間隔Δθ分別為2°，4°，6°，8°以及10°，重建結(jié)果如圖2所示。不難看出，隨著采樣間隔的增加，重建圖像質(zhì)量有所降低，但是OSEM算法在角度抽樣間隔Δθ=8°圖像質(zhì)量，沒有明顯下降。這說明，該算法在采用大角度采樣間隔降低投影數(shù)據(jù)獲取時間的同時，圖像質(zhì)量也可以得到保持。

為進(jìn)一步衡量圖像的重建精度均方根誤差（Root Mean Square Error，RMSE）、重建時間與子集個數(shù)之間的關(guān)系，我們比較了當(dāng)角度抽樣間隔為1°時，圖像重建時間隨不同子集的變化。本文采用Intel(R) Core(TM)i3-2120 CPU@3.30 GHz 型CPU、內(nèi)存4G，計算平臺是Matlab2015a。計算可得當(dāng)子集數(shù)為3時，圖像重建時間已經(jīng)下降至40%，而重建圖像的RMSE值，卻無明顯變化。當(dāng)子集數(shù)大于18時，圖像重建時間隨子集數(shù)增加已無明顯變化，重建質(zhì)量精度開始逐漸降低，由此認(rèn)為，此時最佳子集數(shù)為18。當(dāng)前情況下，重建圖像大小為128×128，共計使用約7 min。顯然，當(dāng)重建圖像較大，像素矩陣接近或超過5個數(shù)量級時，選擇合適的子集數(shù)可以大大地減少圖像重建時間。

4 結(jié)語

本文在詳細(xì)推導(dǎo)XFCT成像原理的基礎(chǔ)上，針對XFCT圖像重建，基于OSEM方法，研究了基于OSEM的XFCT重建算法，進(jìn)行相關(guān)仿真和驗證，討論了重建質(zhì)量、重建時間與子集劃分之間的關(guān)系，優(yōu)選了最優(yōu)子集劃分個數(shù)，這對于指導(dǎo)XFCT重建具有重要的意義。

圖2 不同采樣角度下OSEM-TV算法重建圖像

[參考文獻(xiàn)]

[1]YUASA T，AKIBA M，TAKEDA T，et al.Reconstruction method for fluorescent X-ray computed tomography by least-squares method using singular value decomposition[J].IEEE Transaction Nuclear Science，1997（1）：54-62.

[2]JONES B L，CHO S H.The feasibility of polychromatic cone-beam X-ray fluorescence computed tomography（XFCT）imaging of gold nanoparticle-loaded objects: a Monte Carlo study[J].Physics in Medicine Biology，2011（12）：3719-3730.

[3]RICKETTS K，GUAZZONI C，CASTOLDI A，et al.A bench-top K X-ray fluorescence system for quantitative measurement of gold nanoparticles for biological sample diagnostics[J].Nuclear Instruments and Methods in Physics Research Section A: Accelerators，Spectrometers，Detectors and Associated Equipment，2016（816）：5-32.

[4]莊天戈.CT原理與算法[M].上海：上海交通大學(xué)出版社，1992.

[5]YANG Q，DENG B，LV 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（2）：210-215.

[6]JIANG S H，HE P，DENG L Z，et al.Monte Carlo simulation for polychromatic x-ray fluorescence computed tomography with sheet beam geometry[J].International Journal of J Biomedical Imaging，2017（8）：79.