Sep 15

3D Model-Based Iterative Reconstruction for Synchrotron Tomography

Introduction and Theory

Synchrotron based X-ray tomography is widely used for three dimensional imaging of material samples. Typically, reconstructions are done using analytical reconstruction techniques such as filtered back projection (FBP) or Fourier reconstruction methods like gridrec. However, direct reconstruction of synchrotron data typically result in strong ring and streak artifacts in the reconstruction. The ring artifacts are caused by non-idealities in the measurement system such as dusty scintillator screens, faulty detector elements, etc. The streak artifacts are caused by sudden blinding of detector pixels by high energy photons (called "zingers"). Typically, several pre/post-processing steps are to used to correct for these artifacts.

Our approach to solving the reconstruction problem is to use a model-based iterative reconstruction algorithm (MBIR) [1] to reconstruct the synchrotron data without any pre/post processing steps. Our MBIR algorithm automatically produces artifact free reconstructions by modeling the non-idealities in the measurement system. MBIR uses a model for the measurement system (i.e., a forward model), a model for the noise statistics and a model for the object (i.e., a prior model) to achieve high quality reconstructions. In MBIR, the reconstruction is typically formulated as the following minimization problem
(\hat{x},\hat{\phi}) = argmin_{x,\phi} \{ -log p(y|x,\phi) - log p(x)\}
where y is the projection data, x is the vector of all voxel values, \phi represents the unknown calibration parameters, p(y|x,\phi) is the likelihood function (i.e., forward model), p(x) is the prior model, \hat{x} is the reconstruction and \hat{\phi} is the estimated value of the unknown calibration parameters.

Let y_{n,i} be the projection value at the n^{th} view and i^{th} detector element. Then, it has been shown that the non-idealities which cause ring artifacts are modeled by a detector dependent offset d_i at every view [1, 2]. To model the zinger measurements, we differentially weight the zinger measurements from the normal measurements using a generalized Huber function \beta_{T,\delta}(x). Zinger measurements result in a much larger data mismatch error than the normal measurements. So, to account for zingers, we use the Huber function to penalize the zinger measurements with a large data mismatch error. The generalized Huber function is given by
\beta_{T,\delta}(z) = \begin{cases}
z^{2} & |z| < T \\
2\delta T|z| + T^{2}(1-2\delta) & |z| \geq T
where T,\delta are parameters of the Huber function. We can then show that the log-likelihood term which models the measurement system is given by,
-\log p(y|x,d,\sigma) =
\frac{1}{2}\displaystyle \sum_{n=1}^{N} \sum_{i=1}^{M}
+ MN\log(\sigma) + \tilde{f}(y)
where d=[d_1 \cdots d_M], A_{n,i,*} is the i^{th} row of the forward projection matrix A_n, \Lambda_{n} is the noise matrix, \sigma is a proportionality constant for the noise and \tilde{f}(y) is a constant which is ignored in the subsequent optimization. Here e_{n,i}=y_{n,i}-A_{n,i,*}x_j-d_i is the data mismatch error corresponding to the measurement y_{n,i}. We use a special case of the q-generalized Gaussian Markov Random Field (qGGMRF) [3] as a prior model p(x) for the voxels. The reconstruction \hat{x} is then given by solving the following optimization problem,
(\hat{x},\hat{d},\hat{\sigma}) = argmin_{x,d,\sigma} \{ -\log p(y|x,d,\sigma) - log p(x)\}

Figure 1. Comparison of reconstructions using different methods.

(a) is reconstruction using the gridrec algorithm without any pre/post processing. (b) is reconstruction using conventional MBIR algorithm without any model for the measurement non-idealities. (c) is gridrec with pre/post-processing steps to get rid of ring and streak artifacts. (d) is using the proposed MBIR reconstruction method which models the measurement non-idealities and zingers.


In Figure 1, we compare reconstructions of a ceramic composite material imaged using synchrotron radiation [4]. The object is reconstructed from a dataset consisting of 128 view projections. The gridrec [5] reconstruction without any pre/post-processing steps is shown in Fig. 1 (a).  In Fig. 1(b), we can see the reconstruction using the conventional MBIR algorithm without a model correcting for the ring and streak artifacts. We can see that the MBIR reconstruction has much higher quality than gridrec reconstructions. However, it also has strong ring artifacts and streaks (caused by zingers). In Fig. 1(c), we can see the gridrec reconstructions with pre/post processing steps to correct for the ring and streak artifacts. In Fig. (d), we can see the reconstructions using the proposed MBIR algorithm which models the measurement non-idealities and zingers. We can see that the proposed MBIR reconstruction method produces the best reconstructions without any artifacts. 


[1] K. A. Mohan, S. V. Venkatakrishnan, L. F. Drummy, J. Simmons, D. Y. Parkinson, and C. A. Bouman, “Model-based iterative reconstruction for synchrotron X-ray tomography,” in Acoustics, Speech and Signal Processing (ICASSP), 2014 IEEE International Conference on, May 2014. 

[2] J. Hsieh, Computed Tomography: Principles, Design, Artifacts, and Recent Advances,. SPIE Press, Nov. 2009. 

[3] J. B. Thibault, K. D. Sauer, C. A. Bouman, and J. Hsieh, “A three- dimensional statistical approach to improved image quality for multislice helical CT,” Medical Physics, vol. 34, no. 11, pp. 4526–4544, 2007.

[4] H. A. Bale, A. Haboub, A. A. MacDowell, J. R. Nasiatka, D. Y. Parkinson, B. N. Cox, D. B. Marshall, and R. O. Ritchie, “Real- time quantitative imaging of failure events in materials under load at temperatures above 1,600◦ C,” Nature Materials, vol. 12, no. 1, 2013. 

[5] F. Marone and M. Stampanoni, “Regridding reconstruction algorithm for real-time tomographic imaging,” Journal of Synchrotron Radiation, vol. 19, no. 6, pp. 1029–1037, Nov 2012.