Fast and accurate X-ray fluorescence computed tomography imaging with the ordered-subsets expectation maximization algorithm

Yang, Q.; Deng, B.; Lv, W.; Shen, F.; Chen, R.; Wang, Y.; Du, G.; Yan, F.; Xiao, T.; Xu, H.

doi:10.1107/S0909049511052253

research papers

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 19| Part 2| March 2012| Pages 210-215

doi:10.1107/S0909049511052253

Fast and accurate X-ray fluorescence computed tomography imaging with the ordered-subsets expectation maximization algorithm

Qun Yang,^a Biao Deng,^a Weiwei Lv,^b Fei Shen,^a Rongchang Chen,^c Yudan Wang,^a Guohao Du,^a Fuhua Yan,^b Tiqiao Xiao ^a ^* and Hongjie Xu ^a

^aShanghai Synchrotron Radiation Facility, Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Shanghai 201204, People's Republic of China, ^bDepartment of Radiology, Zhongshan Hospital, Fudan University, Shanghai 200032, People's Republic of China, and ^cINFN, Sezione di Trieste, Trieste 34012, Italy
^*Correspondence e-mail: tqxiao@sinap.ac.cn

(Received 1 August 2011; accepted 3 December 2011; online 5 January 2012)

The ordered-subsets expectation maximization algorithm (OSEM) is introduced to X-ray fluorescence computed tomography (XFCT) and studied; here, simulations and experimental results are presented. The simulation results indicate that OSEM is more accurate than the filtered back-projection algorithm, and it can efficiently suppress the deterioration of image quality within a large range of angular sampling intervals. Experimental results of both an artificial phantom and cirrhotic liver show that with a satisfying image quality the angular sampling interval could be improved to save on the data-acquisition time when OSEM is employed. In addition, with an optimum number of subsets, the image reconstruction time of OSEM could be reduced to about half of the time required for one subset. Accordingly, it can be concluded that OSEM is a potential method for fast and accurate XFCT imaging.

Keywords: X-ray fluorescence; computed tomography; ordered-subsets expectation maximization algorithm; angular sampling interval.

1. Introduction

Synchrotron-radiation-based X-ray fluorescence computed tomography (XFCT), a complement to absorption computed tomography (CT), is a stimulated emission tomography modality. It allows for the reconstruction of element distributions on a virtual section across a sample using various algorithms (Hogan et al., 1991 ; Schroer, 2001 ; La Riviere & Vargas, 2006 ; La Riviere et al., 2010 ; Miqueles & Pierro, 2010 ). Taking advantage of being sensitive to trace-element concentrations and not being prone to introducing contamination during sample preparation, it has been applied in the fields of botany (Kim et al., 2006 ), biology (de Jonge et al., 2010 ), medicine (Takeda et al., 2009 ), materials (Naghedolfeizi et al., 2003 ), mineralogy (Lemelle et al., 2004 ) and so on. Many synchrotron radiation facilities have set up XFCT imaging systems; for example, SPring-8 (Hirai et al., 2007 ), APS (de Jonge et al., 2010), ESRF (Golosio et al., 2004 ), NSLS (Kim et al., 2006) etc. However, XFCT has not been used as a routine tool as yet, especially in biomedical research. The main reason for this is that the long image acquisition time would result in a high radiation dose and may change the element distributions.

To solve this problem a number of detector elements were adopted (McNear Jr et al., 2005 ) to decrease the collection time of each spectrum by expanding the solid angle. A new data acquisition mode (Huo et al., 2009 ) was proposed to accelerate the image acquisition. In the new mode a sheet beam was used, and few translations were needed at each sampling angle. Helical X-ray microtomography has also been combined with XFCT (Golosio et al., 2004). All these methods mentioned above are considered from the aspect of hardware, whereas accelerating the image acquisition from the aspect of algorithms has not been well investigated yet. The aim of this paper is to introduce an algorithm which is able to perform image reconstruction with a smaller number of angular views but keeping the image quality high.

Inspired by the image reconstruction in PET (positron emission tomography) and SPECT (single-photon emission computed tomography), the ordered-subsets expectation maximization algorithm (OSEM) is introduced to XFCT here. It has been proven to be accurate in PET reconstruction and to be able to reconstruct with limited angle projection data (Hudson & Larkin, 2002 ). Rust & Weigelt (1998 ) have also concluded that the maximum-likelihood expectation maximization algorithm (ML-EM) is more accurate than the algebraic reconstruction technique. In fact, ML-EM is essentially OSEM with one subset. As an iterative algorithm, OSEM is time-consuming. However, its convergence rate is theoretically proportional to the number of subsets, so the time for image reconstruction could be reduced remarkably with a proper number of subsets. Therefore, OSEM should be a potential candidate for fast and accurate XFCT imaging.

In this paper a numerical simulation has been carried out using a set of angular sampling intervals (ASIs) ranging from 2° to 10°. Comparisons between OSEM and filtered back-projection (FBP) reconstructions were also performed. Then, XFCT experiments, based on OSEM, on an artificial phantom and a cirrhotic liver sample were carried out at the X-ray imaging beamline (BL13W1) of Shanghai Synchrotron Radiation Facility (SSRF) to evaluate the simulation results. Absorption corrections for incident X-ray and X-ray fluorescence were not performed on the reconstructed images.

2. Mathematical model for XFCT

Fig. 1 depicts a schematic diagram of the acquisition geometry for XFCT, which operates on the basis of first-generation CT. The xy-coordinate system is associated with the object and can be rotated anticlockwise by a projection angle θ around the st-coordinate system that is fixed to the laboratory. Taking point P(s, t) as an example, when a pencil beam with an appropriate energy irradiates the object, it stimulates element i to emit X-ray fluorescence isotropically with an intensity that is proportional to the concentration of element i, the intensity of the X-ray that arrives at point P(s, t) and the yield of X-ray fluorescence. The part of the X-ray fluorescence within the range γ_min–γ_max is detected when escaping from the object. Performing a line integral along the penetration path of the incident X-ray, a single projection value I(θ, t) is given by

$[I(\theta,t)=K{I_0}\textstyle\int\limits_{-\infty}^{+\infty}{\Big[{C_i}\Omega(s,t)\,f(\theta,s,t)\,g(\theta,s,t)\Big]}\,{\rm{d}}s,\eqno(1)]$

where

$[\eqalign{f(\theta,s,t) &=\exp\left\{-\textstyle\int\limits_{-\infty}^s\mu^{\rm{I}}\,\left[x(s^\prime,t^\prime),y(s^\prime,t^\prime)\right]\,{\rm{d}}s^\prime\right\}, \cr g(\theta,s,t)&=\textstyle\int\limits_{{\gamma_{\min}}}^{{\gamma_{\max}}}\exp\left(-\textstyle\int\limits_0^L \mu^{\rm{F}}\,{\rm{d}}l\,\right)\,{\rm{d}}\gamma.}\eqno(2)]$

K is the scaling constant and I₀ is the initial intensity of the incident X-rays. C_i represents the concentration of element i. Ω(s,t) denotes the solid angle. μ^I and μ^F are the distributions of the linear attenuation coefficients at the energies of the incident X-rays and X-ray fluorescence, which could be measured by conventional absorption tomography. To obtain the whole projection, the object is scanned through the incident X-ray by translation along the t axis, and then rotated by an angular sampling interval Δθ. The procedure is repeated until a full rotation is completed over 180°.

Figure 1
Schematic diagram of the typical acquisition geometry for XFCT, where (i, j) represents the pixel located at the ith row and jth column of the pixel matrix, (m, n) represents the nth projection line at the mth projection angle, Δd and Δd_q (q = 1, 2, 3, 4) are the vertical distances from the pixel centre and the centre of a quarter section to the line (m, n), respectively. d is the size of the pixel. L denotes the distance from point P(s, t) to the edge of the object. γ_min–γ_max is the range of the solid angle at which point P(s, t) is viewed by the detector.

For OSEM reconstruction in XFCT, the sample was divided into I × I pixels. M angular views and N translations were employed. The element distribution to be estimated is defined as C(i, j) and the projection values are denoted as I(m, n) (m = 1, 2,…, M; n = 1, 2,…, N). K(i, j, m, n) represents the contribution of pixel (i, j) to I(m, n). S_l (l = 1, 2,…, L) describes the chosen subsets. C^l+1(i, j) and C^l(i, j) are the pixel values corresponding to the subset (l + 1) and subset l. The formula referenced to Hudson & Larkin (2002) was applied,

$[{C^{\,l\,+\,1}}(i,j) = {{{C^{\,l}}(i,j)} \over {\sum\limits_{(m,n)\,\in\,{S_l}} \!\!\!\!\!{K(i,j,m,n)} }}\,\,\sum\limits_{(m,n)\,\in\,{S_l}} \!\!{{{K(i,j,m,n)\,I(m,n)} \over{{p^l}(m,n)}}},\eqno(3)]$

where

$[p^{l}(m,n)=\textstyle\sum\limits_{i,\,j} K(i,j,m,n)\,C^{\,l}(i,j),\qquad(m,n)\in{S_l}.\eqno(4)]$

Here we assumed that the concentration of element i focuses on the centre of each pixel. For the sake of accuracy, each pixel was divided into four parts when calculating the fluorescent intensity,

$[K(i,j,m,n)= \left\{ \matrix{ (n_0/4)\,\Omega(s,t),\hfill &\quad \Delta{d}\le{d_0},\hfill \cr 0,\hfill & \quad\Delta{d}\,\,\gt\,{d_0},\hfill}\right. \eqno(5)]$

where

$[d_0= \left\{ \matrix{ 2^{1/2}\,d\sin(45^\circ+\theta),\hfill & \quad\theta\le90^\circ,\hfill \cr 2^{1/2}\,d\sin(-45^\circ+\theta), & \quad\theta\,\,\gt\,90^\circ\hfill,}\right.\eqno(6)]$

and n₀ was established by the amount Δd_q that meets the requirement Δd_q $[\le]$ d₀.

In XFCT, all the projection values at the same projection angle form a natural group, so the subsets can be obtained by the division of angular views. Thus S_l was composed of all the projection values at the mth [m = l, l + L, l + 2L,…, l + (M/L − 1)L] projection angle.

3. Simulation

3.1. Numerical phantom

Fig. 2 shows a modified Shepp–Logan phantom, in which the skull is not included, for considering the effect of the skull on the comparison of FBP and OSEM (Shepp & Vardi, 1982 ). Fig. 2(a) shows the distribution of Zn with mass concentrations of 0.01%. Figs. 2(b) and 2(c) show attenuation maps at incident X-ray (12 keV) and Zn Kα line energies. Rotating over 180°, a set of projections was acquired with different ASIs. At each projection angle, 100 translation steps were adopted. The detector was assumed to be 10 mm away from the centre of the phantom.

Figure 2
Numerical phantom with 1024 × 1024 pixels. (a) Zn distribution. (b) Attenuation map at the incident X-ray energy. (c) Attenuation map at the Zn characteristic Kα X-ray energy.

3.2. Reconstruction results

To compare the images reconstructed by OSEM with projections corresponding to different ASIs, a factor F was defined,

$[F=\textstyle\sum\limits_j^J\left(C_j^{\,n+1}-C_j^{\,n}\right)^2/J,\eqno(7)]$

where C_j^{n + 1} and C_jⁿ represent the jth pixel values after n + 1 and n iterations, respectively, and J is the number of pixels. A critical value of F was set in advance to stop the iteration process.

Fig. 3 depicts the Zn distributions reconstructed by FBP and OSEM. It is obvious that all OSEM reconstructions show fewer artifacts than FBP reconstructions. In particular, even when the ASI increases to 6°, the three separate ellipses shown at the bottom of Fig. 2(a) are still discriminable in the OSEM reconstructions. However, they become blurred even when Δθ = 2° in the FBP reconstructions.

Figure 3
Reconstructed Zn distributions by FBP (top row) and OSEM (bottom row) corresponding to different ASIs. All of the images are 100 pixels × 100 pixels, and each pixel is 0.02 mm × 0.02 mm.

OSEM reconstructions corresponding to different ASIs, as shown in Figs. 3(h)–3(l), have almost the same visual features except that the uniformity slightly degrades. It indicates that OSEM could achieve acceptable image quality even when a larger ASI value of 6° is used. This implies that the data-acquisition time can be remarkably reduced by the introduction of OSEM, while maintaining a high image quality.

3.2.1. Quantitative analysis

The root-mean-square error (RMSE) between the reconstructed image and the true one was employed to quantitatively evaluate the reconstruction accuracy. RMSEs for all the images shown in Fig. 3 are given in Fig. 4. It is obvious that in the OSEM reconstruction the RMSE value is much less sensitive to the ASI value than in the case of the FBP reconstruction. This means that a high reconstruction precision could be achieved by OSEM when a large ASI is used. For example, the RMSE value for OSEM reconstruction when Δθ = 10° increases by 9.6% compared with that when Δθ = 2°, while for FBP reconstruction the error increases by 42.0%. In addition, compared with Δθ = 2°, the data-acquisition time could be reduced by a factor of three when Δθ = 6°, while the RMSE value for OSEM reconstruction increases by less than 5% which should be endurable for XFCT.

Figure 4
RMSE of OSEM and FBP reconstructions versus ASI.

3.2.2. Optimum subsets number for OSEM

For OSEM, the number of subsets will affect the reconstruction accuracy and speed simultaneously. In principle, the reconstruction time will decrease when the number of subsets increases, while the accuracy remains. Fig. 5 shows the relationship between the reconstruction time and the number of subsets, for Δθ = 2°. It shows that the optimum number of subsets L_opt is 15, where more than half of the reconstructing time could be saved, while the image accuracy remains almost the same. With a larger subset number the reconstruction time could not be reduced, while the image quality becomes deteriorated. Similar results exist when Δθ = 3°, 4°, 5° and 6°, but L_opt is different. However, Δθ × L_opt was found to be equal or smaller than 30°. This simulation was performed on a computer with IntelCore 2 Duo CPU E7300 at 2.66 Hz and 4 GB Memory. The reconstruction time is about 2 min when Δθ = 2° and 30 iterations with one subset were employed. It may not be necessary to accelerate the reconstruction by an optimum number of subsets. However, if the one-dimensional size of the pixel matrix is 1000 pixels or more, the reconstruction time saved by using an optimized subset number will be considerable.

Figure 5
Relationship between the reconstruction time and the number of subsets when Δθ = 2°. The values in parentheses are the RMSEs. The numbers of iterations are 30 with one subset, 15 with two subsets, ten with three subsets, six with five subsets, five with six subsets, three with nine subsets, three with ten subsets, two with 15 subsets and one with 30 subsets.

4. Experiments

4.1. Experimental set-up

Experiments were carried out at beamline BL13W1 of SSRF. As shown in Fig. 6, a white beam from a 2 T wiggler is monochromized by a double-crystal monochromator cooled by liquid nitrogen. A set of slits is used to collimate and cut the monochromatic beam to a pencil beam. An ionization chamber is placed in front of the sample to monitor the incident X-ray flux. Behind the sample is a CCD camera to help facilitate the beam alignment. An ultra-low-energy Ge detector (LEGe, Canberra industries) placed at 90° to the incident beam is employed to collect the X-ray fluorescence. The sample is supported by a high-precision sample stage with six dimensions.

Figure 6
Schematic diagram for the XFCT imaging system.

4.2. Phantom

An artificial phantom was used to evaluate the OSEM in XFCT. Four holes of diameter 3 mm were distributed symmetrically inside a 10 mm-diameter PMMA pole, of which two non-adjacent holes were filled with Cd (the Kα line energy is approximately 23 keV) solution with a molar concentration of 1 µmol ml⁻¹. The energy of the incident X-ray beam was 32 keV and the incident X-ray beam size was 200 µm (H) × 500 µm (V). Over a range of 180°, the phantom was scanned with 52 translation steps of 200 µm at each projection angle. The data-acquisition time for each translation step was 1 s.

Fig. 7 shows the Cd distributions reconstructed by FBP and OSEM for different ASIs. The profiles along the white line as noted in Fig. 7(a) are shown on the right. In the OSEM reconstructions, one subset and six iterations were employed. As mentioned above, the energies of the incident X-rays and Cd Kα line are high, so the absorption effect is small. Thus, the reconstructed phantoms approximate well the true Cd distribution, especially in Fig. 7(h).

Figure 7
Cd distributions reconstructed by FBP [(a)–(d)] and OSEM [(e)–(h)]. All images are composed of 52 pixels × 52 pixels, and each pixel size is 200 µm × 200 µm.

By comparing Figs. 7(a)–7(c) and Figs. 7(e)–7(g) it is obvious that the OSEM reconstructions have much higher image qualities. In addition, according to Figs. 7(d) and 7(h) the reconstructed pixel values by OSEM are closer to zero where there is actually no Cd.

Figs. 7(d) and 7(h) also reveal that the image quality of the FBP reconstruction deteriorates seriously as Δθ increases from 2° to 6°, while the OSEM reconstruction suppresses this effect efficiently. This result confirms the conclusion achieved by simulation which shows that the RMSE value when Δθ = 6° has less than a 5% increase compared with that when Δθ = 2°. As a result, we can conclude that with OSEM the reconstructed image quality is comparable when Δθ = 2° and Δθ = 6°.

4.3. Pathologic sample

Zinc is an essential mineral to human health and is also a necessary element in the liver. So, a naturally dried cirrhotic liver of a SD rat was chosen as a biomedical sample for reconstructing its Zn distribution to evaluate OSEM in XFCT. The incident X-ray beam size was 100 µm (H) × 500 µm (V) and the photon energy of the monochromatic beam was 12 keV. Experiments with ASIs of 2° and 6° were carried out. Over a range of 180°, at each projection angle, the sample was scanned with 23 translations at a step size of 100 µm. The data-acquisition time for one translation step was 3 s. The total acquisition time corresponding to Δθ = 6° was about one-third of that corresponding to Δθ = 2°.

Fig. 8 shows the experimental results. Figs. 8(b) and 8(c) show the Zn distributions reconstructed by OSEM with one subset and six iterations, and the pixel size is 100 µm. Comparing Fig. 8(b) with Fig. 8(c), there was not a distinct difference in the image quality, which means that the data-acquisition time could be greatly reduced by selecting a larger ASI. In this way, XFCT for tissue samples which are sensitive to the radiation dose may be routinely possible with the introduction of OSEM.

Figure 8
Zn distribution in a naturally dried cirrhotic liver of a SD rat. (a) Slice obtained by X-ray CT. (b), (c) XFCT reconstructions by OSEM with 23 pixels × 23 pixels when Δθ = 2° and Δθ = 6°.

In Figs. 8(b) and 8(c), it may be seen that the Zn distribution is inhomogeneous. At the lower corner the Zn concentration is much lower than that at the upper corner. As mentioned above, these two images have not been corrected for absorption, so we could not determine whether the inhomogeneous distribution is due to biological variations or to the failure to correct for absorption. However, it does not conflict with the conclusion reached above.

5. Conclusions

In this paper, OSEM was introduced to XFCT reconstruction. To evaluate the potential of this method, a digital simulation and experiments were carried out. For comparison, results from FBP-based XFCT are also given. Simulation results show that OSEM-based XFCT could effectively reduce artifacts and achieve good image quality even when using relatively larger ASIs. Accordingly, the data-acquisition time can be effectively reduced, which has been confirmed by the experimental results of an artificial phantom and one cirrhotic liver sample. Therefore, we can conclude that the imaging efficiency of XFCT could be greatly improved by the introduction of OSEM. Combined with the new data-acquisition mode (Huo et al., 2009), OSEM-based XFCT has the potential to become a routine X-ray fluorescence CT imaging method and find broad applications in many fields.

Acknowledgements

This work was supported in part by National Basic Research Program of China grant No. 2010CB834301, External Cooperation Program of Chinese Academy of Sciences grant No. GJHZ09058, National Natural Science Foundation of China grants Nos. 10805071 and 81071154, and the Knowledge Innovation Program of the Chinese Academy of Sciences.

References

Golosio, B., Somogyi, A., Simionovici, A., Bleuet, P., Susini, J. & Lemelle, L. (2004). Appl. Phys. Lett. 84, 2199–2201. Web of Science CrossRef CAS Google Scholar
Hirai, Y., Yoneyama, A., Hisada, A. & Uchida, K. (2007). AIP Conf. Proc. 879, 1345–1348. CrossRef CAS Google Scholar
Hogan, J., Gonsalves, R. & Krieger, A. (1991). IEEE Trans. Nucl. Sci. 38, 1721–1727. CrossRef CAS Google Scholar
Hudson, H. & Larkin, R. (2002). IEEE Trans. Med. Imag. 13, 601–609. CrossRef Web of Science Google Scholar
Huo, Q., Sato, H., Yuasa, T., Akatsuka, T., Wu, J., Lwin, T., Takeda, T. & Hyodo, K. (2009). X-ray Spectrom. 38, 439–445. Web of Science CrossRef CAS Google Scholar
Jonge, M. D. de, Holzner, C., Baines, S. B., Twining, B. S., Ignatyev, K., Diaz, J., Howard, D. L., Legnini, D., Miceli, A., McNulty, I., Jacobsen, C. J. & Vogt, S. (2010). Proc. Natl Acad. Sci. 107, 15676–15680. Web of Science PubMed Google Scholar
Kim, S. A., Punshon, T., Lanzirotti, A., Li, L., Alonso, J. M., Ecker, J. R., Kaplan, J. & Guerinot, M. L. (2006). Science, 314, 1295–1298. Web of Science CrossRef PubMed CAS Google Scholar
La Riviere, P. & Vargas, P. (2006). IEEE Trans. Med. Imag. 25, 1117–1129. Web of Science CrossRef Google Scholar
La Riviere, P., Vargas, P., Xia, D. & Pan, X. (2010). IEEE Trans. Nucl. Sci. 57, 234–241. Web of Science CrossRef PubMed Google Scholar
Lemelle, L., Simionovici, A., Truche, R., Rau, C., Chukalina, M. & Gillet, P. (2004). Am. Mineral. 89, 547–553. CAS Google Scholar
McNear, D. H. Jr, Peltier, E., Everhart, J., Chaney, R. L., Sutton, S., Newville, M., Rivers, M. & Sparks, D. L. (2005). Environ. Sci. Technol. 39, 2210–2218. Web of Science CrossRef PubMed CAS Google Scholar
Miqueles, E. X. & De Pierro, A. R. (2010). Phys. Med. Biol. 55, 1007–1024. Web of Science CrossRef PubMed Google Scholar
Naghedolfeizi, M., Chung, J., Morris, R., Ice, G., Yun, W., Cai, Z. & Lai, B. (2003). J. Nucl. Mater. 312, 146–155. Web of Science CrossRef CAS Google Scholar
Rust, G. & Weigelt, J. (1998). IEEE Trans. Nucl. Sci. 45, 75–88. Web of Science CrossRef CAS Google Scholar
Schroer, C. (2001). Appl. Phys. Lett. 79, 1912–1914. Web of Science CrossRef CAS Google Scholar
Shepp, L. & Vardi, Y. (1982). IEEE Trans. Med. Imag. 1, 113–122. CrossRef CAS Google Scholar
Takeda, T., Wu, J., Thet-Thet-Lwin, Huo, Q., Yuasa, T., Hyodo, K., Dilmanian, F. A. & Akatsuka, T. (2009). J. Synchrotron Rad. 16, 57–62. Web of Science CrossRef CAS IUCr Journals Google Scholar

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 19| Part 2| March 2012| Pages 210-215

doi:10.1107/S0909049511052253

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Search IUCr Journals		doi		Advanced search
Author		volume	page

research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Fast and accurate X-ray fluorescence computed tomography imaging with the ordered-subsets expectation maximization algorithm

1. Introduction

2. Mathematical model for XFCT

3. Simulation

3.1. Numerical phantom

3.2. Reconstruction results

3.2.1. Quantitative analysis

3.2.2. Optimum subsets number for OSEM

4. Experiments

4.1. Experimental set-up

4.2. Phantom

4.3. Pathologic sample

5. Conclusions

Acknowledgements

References

research papers