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

Journal logoJOURNAL OF
SYNCHROTRON
RADIATION
ISSN: 1600-5775

Fast iterative reconstruction of data in full interior tomography

crossmark logo

aInstitute for Biomedical Engineering, ETH Zurich, 8092 Zurich, Switzerland, and bSwiss Light Source, Paul Scherrer Institute, 5232 Villigen, Switzerland
*Correspondence e-mail: filippo.arcadu@psi.ch

Edited by S. M. Heald, Argonne National Laboratory, USA (Received 18 May 2016; accepted 6 October 2016)

This paper introduces two novel strategies for iterative reconstruction of full interior tomography (FINT) data, i.e. when the field of view is entirely inside the object support and knowledge of the object support itself or the attenuation coefficients inside specific regions of interest are not available. The first approach is based on data edge-padding. The second technique creates an intermediate virtual sinogram, which is, then, reconstructed by a standard iterative algorithm. Both strategies are validated in the framework of the alternate direction method of multipliers plug-and-play with gridding projectors that provide a speed-up of three orders of magnitude with respect to standard operators implemented in real space. The proposed methods are benchmarked on synchrotron-based X-ray tomographic microscopy datasets of mouse lung alveoli. Compared with analytical techniques, the proposed methods substantially improve the reconstruction quality for FINT underconstrained datasets, facilitating subsequent post-processing steps.

1. Introduction

The term interior tomography (INT) refers to the problem of reconstructing an object function, when its support (S) is not a subset of the field of view (FOV) (Natterer & Wübbeling, 2001[Natterer, F. & Wübbeling, F. (2001). Mathematical Methods in Image Reconstruction. Philadelphia: Society for Industrial and Applied Mathematics.]). This imaging modality is broadly applied in medical screening, material science and biology, as it allows high-resolution investigations of small regions of interest (ROIs).

Since INT projections contain information [\notin] FOV, reconstructions with filtered backprojection (FBP) suffer from a DC shift and low-frequency artifacts (Natterer & Wübbeling, 2001[Natterer, F. & Wübbeling, F. (2001). Mathematical Methods in Image Reconstruction. Philadelphia: Society for Industrial and Applied Mathematics.]), that compromise the quantitativeness of the results and make further post-processing, visualization and rendering of the investigated object complicated (Xiao et al., 2007[Xiao, X., De Carlo, F. & Stock, S. (2007). Rev. Sci. Instrum. 78, 063705.]).

Two different analytical reconstruction methods can address this problem: the differentiated backprojection (DBP) and FBP of edge-padded projections (FBP-E). The idea of DBP (Noo et al., 2004[Noo, F., Clackdoyle, R. & Pack, J. D. (2004). Phys. Med. Biol. 49, 3903-3923.]) is to backproject the derivative of the data and, then, recover the original object by Hilbert transform techniques. This method provides quantitative reconstructions, when S is known and only specific geometries are involved: two opposite boundaries of the object are inside the FOV (Noo et al., 2004[Noo, F., Clackdoyle, R. & Pack, J. D. (2004). Phys. Med. Biol. 49, 3903-3923.]); the FOV exceeds the object only from one side (Defrise et al., 2006[Defrise, M., Noo, F., Clackdoyle, R. & Kudo, H. (2006). Inverse Probl. 22, 1037-1053.]); the FOV [\subset S] and the attenuation coefficients are known in a subregion [R\subset] FOV (Ye et al., 2007[Ye, Y., Yu, H., Wei, Y. & Wang, G. (2007). Intl J. Biomed. Imaging, 2007, 63634.]; Kudo et al., 2008[Kudo, H., Courdurier, M., Noo, F. & Defrise, M. (2008). Phys. Med. Biol. 53, 2207-2231.]; Courdurier et al., 2008[Courdurier, M., Noo, F., Defrise, M. & Kudo, H. (2008). Inverse Probl. 24, 65001.]). FBP-E (Seger, 2002[Seger, M. M. (2002). Proceeding of the Symposium on Image Analysis (SSAB02), Lund, Sweden, pp. 33-36.]; Marone et al., 2010[Marone, F., Münch, B. & Stampanoni, M. (2010). Proc. SPIE, 7804, 780410.]; Kyrieleis et al., 2011[Kyrieleis, A. V., Titarenko, V., Ibison, M., Connolley, T. & Withers, P. J. (2011). J. Microsc. 241, 69-82.]) corresponds to FBP of projections, that have been padded with the edge pixels. This strategy prevents the appearance of low-frequency artifacts within the reconstructed FOV, although results are not fully quantitative.

Iterative methods have also been proposed for INT datasets, characterized by low photon statistics and/or with a low number of views (e.g. low-dose scan). Usually, these underconstrained datasets cannot be reconstructed with sufficient accuracy with analytical methods. The separable paraboloidal surrogate (SPS) (Xu et al., 2011[Xu, Q., Mou, X., Wang, G., Sieren, J., Hoffman, E. A. & Yu, H. (2011). IEEE Trans. Med. Imaging, 30, 1116-1128.]), the maximum-likelihood expectation-maximization (EM-ML) (Rashed & Kudo, 2013[Rashed, E. A. & Kudo, H. (2013). J. Synchrotron Rad. 20, 116-124.]) and the penalized weighted least square (PWLS) method (Zhang et al., 2014[Zhang, H., Huang, J., Ma, J., Bian, Z., Feng, Q., Lu, H., Liang, Z. & Chen, W. (2014). IEEE Trans. Biomed. Eng. 61, 2367-2378.]) have been modified to deal with truncated projections. In Xu et al. (2011[Xu, Q., Mou, X., Wang, G., Sieren, J., Hoffman, E. A. & Yu, H. (2011). IEEE Trans. Med. Imaging, 30, 1116-1128.]), the regularized SPS is initialized through the projection-onto-convex-sets (POCS) (Defrise et al., 2006[Defrise, M., Noo, F., Clackdoyle, R. & Kudo, H. (2006). Inverse Probl. 22, 1037-1053.]); the knowledge of S, the attenuation coefficients in one or multiple subregions and the position of the background pixels are needed for this reconstruction technique. In Rashed & Kudo (2013[Rashed, E. A. & Kudo, H. (2013). J. Synchrotron Rad. 20, 116-124.]), the standard EM-ML is combined with thresholding, acting on pixels of selected ROIs. In Zhang et al. (2014[Zhang, H., Huang, J., Ma, J., Bian, Z., Feng, Q., Lu, H., Liang, Z. & Chen, W. (2014). IEEE Trans. Biomed. Eng. 61, 2367-2378.]), the PWLS is regularized with a complete high-quality scan of the same object. All these methods work for very specific (mostly medical) applications and require a priori knowledge obtained from previous scans. Iterative reconstruction of INT datasets without knowledge of either the object support or the attenuation coefficients inside specific ROIs (case of `full' interior tomography, abbreviated as FINT) has been initiated by Yu & Wang (2009[Yu, H. & Wang, G. (2009). Phys. Med. Biol. 54, 2791-2805.]). In this work, the authors prove that an iterative scheme working with total variation (TV) regularization yields a unique reconstruction of a FINT dataset, in case the object under study is piecewise constant.

At synchrotron imaging beamlines, time-resolved high-resolution investigations of a large variety of dynamic systems spanning different fields [e.g. biology (Walker et al., 2014[Walker, S. M., Schwyn, D. A., Mokso, R., Wicklein, M., Müller, T., Doube, M. & Stampanoni, M. (2014). PLoS Biol. 12, e1001823.]; Hesse et al., 2015[Hesse, B., Varga, P., Langer, M., Pacureanu, A., Schrof, S., Männicke, N., Suhonen, H., Maurer, P., Cloetens, P., Peyrin, F. & Raum, K. (2015). J. Bone Miner. Res. 30, 346-356.]), material sciences (Aagesen et al., 2010[Aagesen, L. K., Johnson, A. E., Fife, J. L., Voorhees, P. W., Miksis, M. J., Poulsen, S. O., Lauridsen, E. M., Marone, F. & Stampanoni, M. (2010). Nat. Phys. 6, 796-800.]; Campi et al., 2015[Campi, G., Bianconi, A., Poccia, N., Bianconi, G., Barba, L., Arrighetti, G., Innocenti, D., Karpinski, J., Zhigadlo, N. D., Kazakov, S. M., Burghammer, M., Zimmermann, M., v, Sprung, M. & Ricci, A. (2015). Nature (London), 525, 359-362.]), energy research (Ebner et al., 2013[Ebner, M., Marone, F., Stampanoni, M. & Wood, V. (2013). Science, 342, 716-720.]; Lu et al., 2016[Lu, J., Lee, Y. J., Luo, X., Lau, K. C., Asadi, M., Wang, H. H., Brombosz, S., Wen, J., Zhai, D., Chen, Z., Miller, D. J., Jeong, Y. S., Park, J. B., Fang, Z. Z., Kumar, B., Salehi-Khojin, A., Sun, Y. K., Curtiss, L. A. & Amine, K. (2016). Nature (London), 529, 377-382.]), earth and environmental sciences (Berg et al., 2012[Berg, S., Ott, H., Klapp, S. A., Schwing, A., Neiteler, R., Brussee, N., Makurat, A., Leu, L., Enzmann, F., Schwarz, J. O., Ersten, M. K., Irvine, S. & Stampanoni, M. (2012). Proc. Natl Acad. Sci. 110, 3755-3759.]; Ganne et al., 2012[Ganne, J., De Andrade, V., Weinberg, R. F., Vidal, O., Dubacq, B., Kagambega, N., Naba, S., Baratoux, L., Jessell, M. & Allibon, J. (2012). Nat. Geosci. 5, 60-65.])] are becoming routine. Due to the large variability of the examined samples and to the large datasets (several tens of TB) often associated with these experiments, efficient iterative reconstruction algorithms, not bounded to a specific application, are highly demanded.

This paper introduces two fast strategies for iterative reconstruction of FINT data, i.e. the FOV [\subset] S, S is unknown as well as the attenuation coefficients inside specific ROIs of the FOV. In this regard, the presented results can be considered a step further on the line of research launched by Yu & Wang (2009[Yu, H. & Wang, G. (2009). Phys. Med. Biol. 54, 2791-2805.]). Differently from previous studies, this work shows the importance of data pre-processing (either edge-padding or differentiation) for analytical and iterative reconstruction of FINT datasets. The proposed iterative strategies are fast and provide high-quality non-quantitative reconstructions, best suited for subsequent post-processing and analysis steps (e.g. segmentation and morphological studies), independently of the imaged system.

1.1. Contributions

The contributions of this manuscript are summarized as follows:

(i) Analyze the importance of FINT data pre-processing through edge-padding and differentiation for analytical and iterative tomographic reconstruction.

(ii) Introduction of two fast forward gridding projectors for iterative reconstruction of FINT data: one implements the derivative of the Radon transform, the other acts as Radon transform combined with edge-padding (Arcadu et al., 2016a[Arcadu, F., Nilchian, M., Studer, A., Stampanoni, M. & Marone, F. (2016a). IEEE Trans. Image Process. 25, 1207-1218.],b[Arcadu, F., Stampanoni, M. & Marone, F. (2016b). Opt. Express. Submitted.]).

(iii) Introduction of the virtual strategy, as a more efficient alternative to specific forward projectors for INT data. This strategy is independent of the chosen iterative scheme and regularization.

(iv) Validation of the proposed methods within the framework of the alternate direction method of multipliers plug-and-play (ADMP) (Venkatakrishnan et al., 2013[Venkatakrishnan, S. V., Bouman, C. A. & Wohlberg, B. (2013). Proceeding of the 2013 IEEE Global Conference on Signal and Information Processing (Globalsip 2013), Austin, TX, USA, 3-5 December 2013, pp. 945-948.]), on simulated and real datasets of full interior X-ray tomographic microscopy.

The edge-padding and virtual strategies studied in this work are neither bounded to the ADMP nor applicable only to the case of piecewise constant objects, but they can be utilized within any kind of iterative reconstruction algorithm and regularization scheme.

2. Reconstruction artifacts in interior tomography

2.1. Preliminaries

This work focuses on the reconstruction of a two-dimensional slice from line projections acquired in parallel beam geometry. However, results can be generalized to more complex tomographic configurations (e.g. fan- and cone-beam).

The collection of line projections, measured at different angles [\theta\in[0,\pi)], is called a sinogram. The object to be reconstructed is a finite integrable real function, whose support, S, is assumed to be compact.

In interior tomography, the object support lies outside the FOV, i.e. [S\not\subseteq] FOV and [S\,\,\cap] FOV [\neq\emptyset]. Different INT configurations, listed by Defrise et al. (2006[Defrise, M., Noo, F., Clackdoyle, R. & Kudo, H. (2006). Inverse Probl. 22, 1037-1053.]) and displayed in Fig. 1[link], have been studied in the literature. The case of full interior tomography, often characterizing microtomographic scans (Fig. 1d[link], FOV [\subset S]), is considered in this work.

[Figure 1]
Figure 1
Different interior tomography configurations. The black line represents the boundary of the object support S; the red line delimits the FOV. The striped area identifies [S\,\cap] FOV. (a) Standard tomography. (b) INT with two opposite sides of [\partial S\subset] FOV (Noo et al., 2004[Noo, F., Clackdoyle, R. & Pack, J. D. (2004). Phys. Med. Biol. 49, 3903-3923.]; Zhang et al., 2014[Zhang, H., Huang, J., Ma, J., Bian, Z., Feng, Q., Lu, H., Liang, Z. & Chen, W. (2014). IEEE Trans. Biomed. Eng. 61, 2367-2378.]). (c) INT with the FOV exceeding S only from one side (Defrise et al., 2006[Defrise, M., Noo, F., Clackdoyle, R. & Kudo, H. (2006). Inverse Probl. 22, 1037-1053.]). (d) Full INT (FINT), where FOV [\subset S] (Courdurier et al., 2008[Courdurier, M., Noo, F., Defrise, M. & Kudo, H. (2008). Inverse Probl. 24, 65001.]; Xu et al., 2011[Xu, Q., Mou, X., Wang, G., Sieren, J., Hoffman, E. A. & Yu, H. (2011). IEEE Trans. Med. Imaging, 30, 1116-1128.]; Rashed & Kudo, 2013[Rashed, E. A. & Kudo, H. (2013). J. Synchrotron Rad. 20, 116-124.]; Yu & Wang, 2009[Yu, H. & Wang, G. (2009). Phys. Med. Biol. 54, 2791-2805.]). This work focuses on the latter configuration.

2.2. Case of analytical reconstruction methods

FBP reconstruction of a FINT sinogram is affected by a DC shift and low-frequency artifacts (Seger, 2002[Seger, M. M. (2002). Proceeding of the Symposium on Image Analysis (SSAB02), Lund, Sweden, pp. 33-36.]; Marone et al., 2010[Marone, F., Münch, B. & Stampanoni, M. (2010). Proc. SPIE, 7804, 780410.]; Kyrieleis et al., 2011[Kyrieleis, A. V., Titarenko, V., Ibison, M., Connolley, T. & Withers, P. J. (2011). J. Microsc. 241, 69-82.]), due to the effect of the ramp filter on truncated projections (Seger, 2002[Seger, M. M. (2002). Proceeding of the Symposium on Image Analysis (SSAB02), Lund, Sweden, pp. 33-36.]). To clarify this point, the example discussed by Seger (2002[Seger, M. M. (2002). Proceeding of the Symposium on Image Analysis (SSAB02), Lund, Sweden, pp. 33-36.]) is reproposed here.

Fig. 2[link] shows a complete and truncated projection of a homogeneous circle before and after ramp filtering. In standard tomography (Fig. 2[link], left), the filtered projection has a constant, positive profile in the middle with symmetric negative tails. Once all projections in [[0,\pi)] are smeared back to the image grid, the circle is filled with constant pixel values and the negative tails of each projection are zeroed-out by the positive contributions from all other projections (Fig. 3b[link]). In interior tomography, the filtered projection (Fig. 2[link], right) is instead everywhere positive, characterized by a bowl-shaped profile and tails with very high values. After backprojection, these positive tails are still present, because no negative compensation in the reconstruction process is possible, yielding the bowl artifact displayed in Fig. 3(c)[link]. The DC shift and bowl-shaped profile, shown in Fig. 3(d)[link], can severely compromise the interpretation/morphological analysis of reconstructed FINT datasets.

[Figure 2]
Figure 2
Illustrative example of the effect of the ramp filter on a single projection in the case of standard and interior tomography. Notation: Ft ([{\cal F}^{\,{-1}}_{{\omega}}]) is the (inverse) Fourier transform along the variable t(ω); [|\omega|] is the ramp filter; [P_{{\theta}}(t)] is the projection in parallel beam geometry acquired at angle θ.
[Figure 3]
Figure 3
(a) Simulated homogeneous circle: the background pixels are set to 0.0, the pixels inside the circle are set to 1.0. The image size is 512 × 512 pixels. The dashed red circle identifies the FOV for the interior tomography case. (b) FBP reconstruction of a sinogram with 800 views × 512 pixels in standard tomography. (c) FBP reconstruction of a sinogram with 800 views × 200 pixels in interior tomography with FOV depicted in (a). The image is zoomed to have the same dimensions as (a) and (b). (d) Comparison of the line profiles along the segment D, indicated in (a), for the reconstructions shown in (b) (blue line) and (c) (green line).

2.3. Case of iterative reconstruction methods

Although iterative algorithms do not involve explicit ramp filtering, reconstructions of FINT data are affected by the same artifacts as for analytical methods. Fig. 4[link] shows an iterative reconstruction with SPS of the same sinogram used for Fig. 3(c)[link] and the corresponding line profile along the central row. The image is once again affected by a DC shift and a bowl artifact (Fig. 4b[link]). The reconstruction with any other iterative algorithm would present the same problems.

[Figure 4]
Figure 4
(a) Iterative reconstruction with SPS of the FINT sinogram used also for Fig. 3(c)[link]. (b) Profile of the central line in the iterative reconstruction (a).

This observation can be explained by the fact that iterative methods `mimic' the effect of the ramp filter: the reconstruction after the first iteration corresponds to a blurred version of the object, similar to the result of backprojection when omitting the filtering step; after several iterations, the object becomes sharper, until a similar spatial resolution as for the FBP reconstruction is reached. This `mimicking' behavior leads to the same problems characterizing FBP, when dealing with FINT data.

3. Proposed approach

The proposed approach builds upon existing methods, extended and combined to address the FINT problem. Although here for illustration purposes, specific projectors, iterative scheme and regularization have been chosen, the proposed strategy is more general and not bounded to these choices.

3.1. Flexible iterative reconstruction scheme

The alternate direction method of multipliers plug-and-play (ADMP) (Venkatakrishnan et al., 2013[Venkatakrishnan, S. V., Bouman, C. A. & Wohlberg, B. (2013). Proceeding of the 2013 IEEE Global Conference on Signal and Information Processing (Globalsip 2013), Austin, TX, USA, 3-5 December 2013, pp. 945-948.]) offers an optimization framework, where the reconstructive and the regularization tasks are neatly separated in two subproblems. This structure allows the direct use of any forward projector for the reconstructive subproblem and any denoising method for the regularization subproblem. The ADMP can be viewed as a particular case of the classical formulation of the alternate direction method of multipliers (Boyd et al., 2011[Boyd, S., Parikh, N., Chu, E., Peleato, B. & Eckstein, J. (2010). Mach. Learn. 3, 1-122.]).

In parallel beam geometry, the tomographic problem for a piecewise constant object has the following form,

[\eqalign{ & \tilde{{x}} \, = \, \mathop {\rm{argmin}}_{\bf{x}}\, R({{\bf{x}}}) \,\big|\, {\bf{A}}{{\bf{x}}}_{i} = {\bf{b}}_{i} \quad {\rm{for}} \,\, i = 0,\ldots,n_{z}-1, \cr& {{\bf{x}}} = \left({{\bf{x}}}_{0}^{T},\ldots,\,{{\bf{x}}}_{n_{z}-1}^{T}\right)^{T} \in {{\bb{R}}}^{n_{z}\times n\times n}, \quad {{\bf{x}}}_{i} \in {\bb{R}}^{n \times n}, \cr& \quad {\bf{A}}\in{\bb{R}}^{m \times n,n \times n}, \quad {\bf{b}}_{i}\in{\bb{R}}^{m \times n},} \eqno(1)]

where m is the number of views, nz is the number of slices, n is the number of pixels along a row or a column (assuming the image grid to be square), [{\bf{x}}] is the unknown three-dimensional object, [{\bf{x}}_{i}] is the unknown ith slice, [{\bf{A}}] is the forward projection operator ([{\bf{A}}^{\dagger}] is the adjoint operator or backprojector), [{\bf{b}}] is the sinogram and R is the functional enforcing the object to be piecewise constant. The dual variable [{\bf{u}}] is introduced to transform the constrained optimization problem (1)[link] into

[\eqalignno{ & \mathop {\rm{argmin}}_{{\bf{x}},{\bf{u}}} \,F({\bf{x}},{\bf{u}}) = \mathop {\rm{argmin}}_{{\bf{x}},{\bf{u}}} \left\{ {{1}\over{2}}\,\sum\limits_{i\,=\,0}^{n_{z}-1}\left\| {\bf{A}}{\bf{x}}_{i} - {\bf{b}}_{i} \right\|_{2}^{2} \,\,+\,\, \lambda R({\bf{u}}) \right\} \cr& \quad\quad{\rm{subject\,\,to}}\,\,{\bf{u}} = {\bf{x}}. &(2)}]

The constraint [{\bf u}] = [{\bf x}] is incorporated in the functional by augmenting F with an array of multipliers [\boldgamma]. In this way, the following Lagrangian is obtained,

[\eqalignno{ & \mathop {\rm{argmin}}_{{\bf{x}},{\bf{u}},\boldgamma} \,{\cal L}({\bf{x}},{\bf{u}},{\boldgamma}) = \cr& \mathop {\rm{argmin}}_{{\bf{x}},{\bf{u}},{\boldgamma}} \left\{ {{1}\over{2}} \sum\limits_{i\,=\,0}^{n_{z}-1}\left\| {\bf{A}}{\bf{x}}_{i} - {\bf{b}}_{i} \right\|_{2}^{2} \,\,+\,\,\lambda R({\bf{u}}) + {{\mu}\over{2}}\| {\bf{x}} - {\bf{u}} + {\boldgamma} \|_{2}^{2} \right\}. &(3)}]

The original problem (1)[link] is, therefore, mapped into the minimization of the Lagrangian (3)[link]. The ADMP solves (3)[link] by cyclically minimizing two subproblems with respect to the variable [{\bf{x}}] and [{\bf{u}}] and by updating the multipliers [{\boldgamma}] until convergence,

[\left\{ \eqalign{ & {\bf{x}}^{(k+1)} \quad\longleftarrow\quad \mathop {\rm{argmin}}_{{\bf{x}}}\,{\cal L} \left({\bf{x}},{\bf{u}}^{(k)},{\boldgamma}^{(k)}\right)_{\vphantom{\Big|}} \cr & {\bf{u}}^{(k+1)} \quad\longleftarrow\quad \mathop {\rm{argmin}}_{{\bf{u}}}\,{\cal L} \left({\bf{x}}^{(k+1)},{\bf{u}},{\boldgamma}^{(k)}\right)_{\vphantom{\Big|}} \cr& {\boldgamma}^{(k+1)} \quad\longleftarrow\quad {\boldgamma} ^{(k)} + \left({\bf{x}}^{(k+1)} - {\bf{u}}^{(k+1)} \right) } \right. \eqno(4)]

where the superscript (k) refers to the kth iteration of a selected variable. The first and third terms of (3)[link] will contribute to the [{\bf x}]-subproblem, whereas the second and the third terms will contribute to the [{\bf u}]-subproblem. To explicitly derive the form of the [{\bf x}]-subproblem, the gradient of L with respect to [{\bf x}] is calculated,

[\eqalignno{ {\bb{R}}^{n \times n} \ni 0 & = \nabla_{{\bf{x}}_{i}} \, {\cal L} \left({\bf{x}},{\bf{u}}^{(k)},{\boldgamma}^{(k)}\right) \cr & = \nabla_{{\bf{x}}_{i}} \left\{ {{1}\over{2}} \sum\limits_{i\,=\,0}^{n_{z}-1} \left(\left\| {\bf{A}}{\bf{x}}_{i} - {\bf{b}}_{i} \right\|_{2}^{2}\,\, +\,\, \mu\big\| {\bf{x}}_{i} - {\bf{u}}_{i}^{(k)} + {\boldgamma}_{i}^{(k)} \big\|_{2}^{2} \right)\right\} \cr & = {\bf{A}}^{\dagger}\left({\bf{A}}{\bf{x}}_{i}-{\bf{b}}_{i}\right) + \mu \left({\bf{x}}_{i} - {\bf{u}}_{i}^{(k)} + {\boldgamma}_{i}^{(k)} \right) \cr & = \left({\bf{A}}^{\dagger}{\bf{A}} + \mu{\bf{I}} \right) {\bf{x}} - \left[{\bf{A}}^{\dagger}{\bf{b}} + \mu \left({\bf{u}}_{i}^{(k)} - {\boldgamma}_{i}^{(k)} \right) \right] \cr & = \tilde{{\bf{A}}}{\bf{x}} - \tilde{{\bf{b}}} .&(5)}]

Since [\tilde{{\bf{A}}}] is symmetric and positive-definite, the conjugate gradient (CG) technique (Hestenes & Stiefel, 1952[Hestenes, M. R. & Stiefel, E. (1952). J. Res. Natl. Bur. Stan. 49, 409-436.]) is adopted to iteratively find the solution of the system [\tilde{{\bf{A}}}{\bf{x}}] = [\tilde{{\bf{b}}}]. For the [{\bf{u}}]-subproblem, it results that

[\eqalignno{ \mathop{\rm{argmin}}_{{\bf{u}}}\, {\cal L} \left({\bf{x}}^{(k+1)},{\bf{u}},{\boldgamma}^{(k)}\right) & = \mathop {\rm{argmin}}_{{\bf{u}}} \, \Big\{ \lambda R({\bf{u}}) + {{\mu}\over{2}}\| {\bf{x}}^{(k+1)} \cr& \quad\,\,- {\bf{u}} + {\boldgamma}^{(k)} \|_{2}^{2} \Big\}&(6) \cr & = \mathop{\rm{argmin}}_{{\bf{u}}} \, \left\{ {{1}\over{2}}\| {\bf{u}} - \tilde{{\bf{u}}} \|_{2}^{2} + {{\lambda}\over{\mu}} R({\bf{u}}) \right\}. }]

The last term in (6)[link] corresponds exactly to a denoising problem, where [\tilde{{\bf u}}] represents the input noisy image and [\tau] = [\lambda/\mu] is the regularization strength. The form of the [{\bf u}]-subproblem gives the `plug-and-play' qualification to the ADMP, as the type of denoising can be changed without altering the structure of the entire algorithm.

The iterative procedure is stopped when the L2-norm of the relative difference between reconstructions of subsequent iterations, [{\bf x}^{{(k)}}] and [{\bf x}^{{(k+1)}}], is smaller than a certain threshold, i.e. [\|{\bf x}^{{(k+1)}}-{\bf x}^{{(k)}}\|_{{2}}^{{2}}/\|{\bf x}^{{(k)}}\|_{{2}}^{{2}}\,\,\lt\,\,\varepsilon] = 0.01. One of the main advantages of ADMM-type methods is that satisfactory results in terms of image quality can be achieved after very few iterations.

For the experiments presented in §4[link], ADMP with split Bregman total variation (SBTV) (Goldstein & Osher, 2009[Goldstein, T. & Osher, S. (2009). SIAM J. Imaging Sci. 2, 323-343.]) as plug-and-play regularization has been used. τ strongly determines the quality of the iterative reconstruction since it controls the trade-off between spatial resolution and noise removal.

3.2. Fast forward projectors for interior tomography

The DC shift compromises the quantitativity of FINT reconstructions, but, except for an offset constant throughout the tomographic slice, the information is preserved. The bowl artifact leads, instead, to non-quantitative reconstructions characterized by varying grey level values within homogeneous regions and prevents, therefore, any reliable morphological analysis, unless `decupping' algorithms or sophisticated segmentation approaches are utilized.

Two analytical methods can be adopted to avoid this bowl artifact: DBP and FBP-E. Reconstructions with these techniques are non-quantitative, i.e. exact up to an unknown constant, but can be safely used for image analysis if the knowledge of the attenuation coefficients in an absolute sense is not needed, as is often the case. Fig. 5[link] shows the reconstruction of a FINT sinogram created from a Shepp–Logan (SL) phantom using FBP, DBP and FBP-E. The original SL with 2048 × 2048 pixels (Fig. 5a[link]) is forward projected over 800 angles in [[0,\pi)] and only the central 512 pixels of the sinogram are retained. The bowl artifact is visible at the corners of the FBP reconstruction (Fig. 5c[link]) and in a line profile (Fig. 5f[link], red). This artifact is instead not present in the DBP and FBP-E reconstructions (Figs. 5d, 5e; Fig. 5f[link], green, black). For the iterative reconstruction of FINT data, the proposed idea is to use tomographic forward operators based on the principles of DBP and FBP-E, i.e. a projector that implements the derivative of the Radon transform and a Radon transform acting on edge-padded projections, respectively.

[Figure 5]
Figure 5
(a) Complete SL phantom. (b) SQRES of the SL phantom. (c) SQRES of the FBP reconstruction. (d) SQRES of the DBP reconstruction. (e) SQRES of the FBP-E reconstruction. (f) Line profiles along segment D. The line profiles in (f) show that DBP and FBP-E successfully remove the bowl artifact. The line profiles are manually shifted along the vertical axis to better distinguish one from another.

To guarantee fast reconstructions, the forward gridding projector (FGP) (Arcadu et al., 2016a[Arcadu, F., Nilchian, M., Studer, A., Stampanoni, M. & Marone, F. (2016a). IEEE Trans. Image Process. 25, 1207-1218.]), that works in the Fourier space and has a complexity of O(N2log2N), has been selected for this work. Studies conducted by Arcadu et al. (2016a[Arcadu, F., Nilchian, M., Studer, A., Stampanoni, M. & Marone, F. (2016a). IEEE Trans. Image Process. 25, 1207-1218.]) showed that the FGP is significantly faster than standard space-based projectors [complexity of O(N3)] (Toft, 1996[Toft, P. (1996). PhD thesis, Denmark Technical University, Denmark.]) and its accuracy results comparable with that of very sophisticated operators, when used in iterative schemes.

In standard tomography, the FGP works on an oversampled grid created with zero-padding. For FINT data, the following modifications of the original forward gridding projector are proposed: (i) the differentiated FGP (DFGP) (Arcadu et al., 2016b[Arcadu, F., Stampanoni, M. & Marone, F. (2016b). Opt. Express. Submitted.]), that implements the derivative of the Radon transform (still works on an oversampled grid created with zero-padding) and computes a differential sinogram; (ii) the FGP combined with edge-padding (FGP-E), i.e. the oversampled grid is created by edge-padding of the object. If these two operators are used within the ADMP framework, two iterative schemes can be defined: ADMP-D implementing the DFGP and ADMP-E implementing FGP-E.

As a proof of principle, the experiment in Fig. 5[link] is repeated with the standard ADMP, ADMP-D and ADMP-E. The bowl artifact, appearing in the reconstruction with the standard ADMP (Fig. 6c[link]), is not present in the images computed with ADMP-D (Fig. 6d[link]) and ADMP-E (Fig. 6e[link]).

[Figure 6]
Figure 6
(a) Complete SL phantom. (b) SQRES of the SL phantom. (c) SQRES of the ADMP reconstruction. (d) SQRES of the ADMP-D reconstruction. (e) SQRES of the ADMP-E reconstruction. (f) Line profiles along segment D. The line profiles in (f) show that ADMP-D and ADMP-E successfully remove the bowl artifact. The line profiles are manually shifted along the vertical axis to better distinguish one from another.

Edge-padding has two main advantages with respect to the differentiated operator. First, the computation of the differential sinogram by finite differentiation is very likely to enhance the noise affecting the original sinogram; if a sophisticate technique is used instead [like the Savitzky–Golay filter (Savitzky & Golay, 1964[Savitzky, A. & Golay, M. J. E. (1964). Anal. Chem. 36, 1627-1639.])], additional parameters, ruling the trade-off between noise and spatial resolution, are added to the reconstruction problem. This argument is valid for both analytical and iterative reconstructions. Second, the number of sub-iterations required by the CG step in the ADMP is related to the conditioning number (CN) of the operator [{\bf A}]: for ADMP-D, [{\bf A}] corresponds to the derivative of the Radon transform; for ADMP-E, [{\bf A}] is the Radon transform itself, having a lower CN than its derivative. It has been experimentally observed that to reach convergence with the same number of iterations the CG loop of ADMP-D requires at least 15 sub-iterations, whereas four are enough for the CG loop of ADMP-E. For these reasons, only ADMP-E is utilized for the experiments discussed in §4[link].

3.3. Iterative virtual method

For the iterative reconstruction of FINT data without a priori knowledge of the object support or of the attenuation coefficients inside specific ROIs, an alternative four-step strategy is also proposed: the iterative virtual method (abbreviated as ADMP-V). The steps of ADMP-V are (Fig. 7[link]): (i) data reconstruction by an analytical method, like DBP or FBP-E; (ii) zeroing of all pixels outside the reconstruction circle; (iii) forward projection of this newly computed image to obtain a `virtual' sinogram, simulating a non-FINT dataset; (iv) reconstruction with ADMP using projectors for standard tomography (FGP, in this case) and using physical constraints (e.g. zeroing all pixels outside the reconstruction circle at each iteration). After steps (i), (ii) and (iii), the initial FINT dataset is transformed into a standard tomographic dataset with known support, that can be reconstructed with any analytical or iterative algorithm.

[Figure 7]
Figure 7
Diagram showing the steps of the virtual strategy (ADMP-V). (i) Reconstruction with FBP-E. (ii) Zero-out all pixels [\notin] FOV. (iii) Forward projection to obtain the virtual sinogram ([m_{\rm{v}}] = [n\pi/2]). (iv) The virtual sinogram is used as input for ADMP. The [\alpha^{{\rm{(an)}}}]-padding used for FBP-E is coloured in light blue, whereas the [\alpha^{{\rm{(int)}}}]-padding inside ADMP is coloured in light red.

4. Experiments

4.1. Analysis framework

Metrics to assess the reconstruction accuracy are computed over the square inscribed inside the reconstruction circle (SQRES). The contrast-to-noise ratio (CNR) is defined as

[{\rm{CNR}} = {{ |S_{{r_{{i1}}}}-S_{{r_{{i2}}}}| }\over{ \sigma_{{r_{{i1}}}}+\sigma_{{r_{{i2}}}} }}, \eqno(7)]

where Sri and [\sigma_{{r_{{i}}}}] are the mean and standard deviation of the ROI ri, assumed to be homogeneous. For the computation of the CNR, ri1 and ri2 must be neighbouring. The values of CNR, reported in the tables of §4[link], represent averages over multiple ROIs at different distances from the image centre.

In the case of simulated data, the peak signal-to-noise ratio (PSNR) (Huyn-Thu & Ghanbaru, 2008[Huyn-Thu, Q. & Ghanbaru, M. (2008). Electron. Lett. 44, 800-801.]) and the mean structural similarity index (MSSIM) (Wang et al., 2004[Wang, Z., Bovik, A. C., Sheikh, H. R. & Simoncelli, E. P. (2004). IEEE Trans. Image Process. 13, 600-612.]) are also used as figures of merit. Since results are not quantitative due to the DC shift, a linearly regressed version, [I_{{\rm{regr}}}], of the reconstruction, I, is used as input for the analysis. Calling O the phantom, [I_{{\rm{regr}}}] is computed as

[I_{{\rm{regr}}} = a\cdot I+b \quad {\rm{such\,\,that}} \, \left\Vert I_{{\rm{regr}}}-O\right\Vert^{{2}}_{{2}} \,\, {\rm{is\,\,minimized}}, \eqno(8)]

where [a,b\,\in{\bb{R}}]. In this way, PSNR and MSSIM scores are not biased by the DC shift. Edge profiles and histograms are also displayed to provide additional information about the spatial resolution and the segmentability of the images. Reconstructions are mapped to the interval [0, 1] to facilitate the display and direct comparison of profiles and histograms.

4.2. Data and settings

The phantoms used to generate the simulated data are displayed in Fig. 8[link]. The first dataset is the SL phantom (Shepp & Logan, 1974[Shepp, L. & Logan, B. F. (1974). IEEE Trans. Nucl. Sci. 21, 2143.]), which is often utilized to validate tomographic reconstruction algorithms. The second dataset is a segmented reconstruction of mouse lung tissue and it is labelled MLT. This phantom is evidently characterized by a high structural complexity and is related to the real data used in §4.7[link]. The procedure to create FINT sinograms from a reference image has been described in §3.2[link].

[Figure 8]
Figure 8
Phantoms used to generate the simulated FINT sinograms. (a) SL phantom. (b) MLT phantom.

The real datasets have been acquired at the TOMCAT beamline of the Swiss Light Source in the framework of the SNF project `in vivo study of lung physiology with fast X-ray tomographic microscopy' (Lovric et al., 2013[Lovric, G., Barré, S. F., Schittny, J. C., Roth-Kleiner, M., Stampanoni, M. & Mokso, R. (2013). J. Appl. Cryst. 46, 856-860.]). The three sinograms correspond to scans of mouse lung tissue at the micrometre scale with an effective detector pixel size of 2.9 µm (MOUSE-1, MOUSE-2) and 1.1 µm (MOUSE-3). The first two datasets consist of 901 Paganin phase-retrieved (Paganin et al., 2002[Paganin, D., Mayo, S., Gureyev, T., Miller, P. & Wilkins, S. W. (2002). J. Microsc. 206, 33-40.]) projections and 2016 pixels, whereas the third one has 501 projections × 2016 pixels.

The regularization strength τ (the only free parameter in ADMP) is tuned so that reconstructions with ADMP-E and ADMP-V look visually as similar as possible.

4.3. Optimal edge-padding length for analytical and iterative reconstructions

The edge-padding length rules the trade-off between the reconstruction accuracy and computational efficiency for both analytical and iterative reconstruction methods. An insufficient amount of padding fails at removing the bowl artifact, whereas an excessive padding substantially increases the computational cost with insignificant gain in the reconstruction quality.

In this work, analytical reconstructions of FINT data are performed by means of the ramp-filtered adjoint FGP-E. The method, labelled GRID-E, is equivalent to FBP-E. To improve the signal-to-noise ratio (SNR) of the retrieved image, projections are additionally windowed with a Hamming function superimposed on the ramp filter. GRID-E is also used in the first step of ADMP-V.

When performing analytical reconstructions with GRID-E, the edge-padding required for the INT dataset and the oversampling needed by the gridding backprojector coincide. The edge-padding or oversampling factor, i.e. the ratio between the number of pixels of edge-padded and original projections, used for GRID-E is indicated with [\alpha^{{\rm{(an)}}}] (the superscript stands for `analytical'). Fig. 9[link] shows the PSNR and MSSIM scores as a function of [\alpha^{{\rm{(an)}}}] for the reconstructions of a SL and MLT sinograms with 1500 views × 512 pixels created from initial phantoms with 4096 × 4096 pixels. The smallest [\alpha^{{\rm{(an)}}}] corresponding to the highest values of both PSNR and MSSIM in Fig. 9[link] is 2.32 (marked with a dashed red line).

[Figure 9]
Figure 9
PSNR and MSSIM of the reconstructions of a FINT sinogram of SL and MLT phantoms computed by GRID-E as a function of [\alpha^{{(\rm{an})}}], the edge-padding factor. The dashed red line locates the optimal oversampling factor, [\alpha^{{(\rm{an})}}] = 2.32.

An iterative reconstruction method utilizing the FGP-E as forward projector (like ADMP-E) depends on two distinct edge-padding factors: [\alpha^{{\rm{(int)}}}] (the superscript stands for `internal'), corresponding to the oversampling factor required by each call of the gridding implementations of [{\bf A}] and [{\bf A}^{{\dagger}}] for both standard and interior tomography; [\alpha^{{\rm{(ext)}}}] (the superscript stands for `external'), which is the edge-padding factor required to address the reconstruction of FINT datasets, can be considered a simple data extrapolation approach and is used for the entire duration of the iterative procedure. A sinogram [{\bf b}], reconstructed by ADMP-E, is first padded by a factor [\alpha^{{\rm{(ext)}}}], becoming [{\bf b}^{{\prime}}] [[\beta] = (number of pixels b′)/(number of pixels b)]. Then, every time [{\bf A}^{{\dagger}}] (analogously it works for [{\bf A}]) is invoked inside the CG, [{\bf b}^{{\prime}}] is padded again by a factor [\alpha^{{\rm{(int)}}}], becoming [{\bf b}^{{\prime\prime}}]. The second padding is only temporary (since it is a requirement of the gridding operators) and, once [{\bf A}^{{\dagger}}] ([{\bf A}]) has completed its own calculation, the resulting image (sinogram) is cropped to remove the additional [\alpha^{{\rm{(int)}}}]-padding. The [\alpha^{{\rm{(ext)}}}]-padding, instead, remains for the entire run of ADMP-E. The double-padding mechanism required by ADMP-E is shown in Fig. 10[link].

[Figure 10]
Figure 10
Difference between the two edge-paddings required by ADMP-E. The edge-padding by a factor [\alpha^{{\rm{(ext)}}}] [[\alpha^{{\rm{(ext)}}}]-padding, indicated in the scheme with a light green colour] is carried out before the start of ADMP-E and remains for the entire duration of the iterative reconstruction. The edge-padding by a factor [\alpha^{{\rm{(int)}}}] [[\alpha^{{\rm{(int)}}}]-padding, indicated in the scheme with a light red colour] is required every time the forward gridding projector, [{\bf A}], and its adjoint operator, [{\bf A}^{{\dagger}}], are called inside ADMP-E. The [\alpha^{{\rm{(int)}}}]-padding is summed up to the [\alpha^{{\rm{(ext)}}}]-padding when making use of [{\bf A}] and [{\bf A}^{{\dagger}}].

Fig. 11[link] shows the PSNR and MSSIM maps as a function of [\alpha^{{\rm{(int)}}}] and [\alpha^{{\rm{(ext)}}}] for the reconstruction of the same SL sinogram as in the experiment of Fig. 9[link]. The maps clearly point out that the reconstruction accuracy strongly varies with [\alpha^{{\rm{(ext)}}}] and only in a weaker way with [\alpha^{{\rm{(int)}}}]. We selected [\alpha^{{\rm{(int)}}}] = 1.7 and [\alpha^{{\rm {(ext)}}}] = 1.87 as optimal edge-padding parameters, because they are the smallest values (highest computational efficiency) guaranteeing simultaneously maximum accuracy in terms of both PSNR and MSSIM.

[Figure 11]
Figure 11
PSNR and MSSIM maps of the reconstructions of a FINT sinogram of a SL phantom computed by ADMP-E as a function of [\alpha^{{(\rm{int})}}] and [\alpha^{{(\rm{ext})}}].

The optimal edge-padding factors determined here are used for the following reconstructions with ADMP-V and ADMP-E. For ADMP-V, [\alpha^{{\rm{(an)}}}] = 2.32 is utilized inside GRID-E and [\alpha^{{\rm{(int)}}}] = 1.7 inside ADMP. ADMP-E utilizes [\alpha^{{\rm{(ext)}}}] = 1.87 as external and [\alpha^{{\rm {(int)}}}] = 1.7 as internal edge-padding factors.

4.4. Validation for different zoom-in factors

The reconstruction accuracy of ADMP-E and ADMP-V has been tested for different zoom-in factors (ZIFs), defined as the ratio between the number of pixels of the sinogram in standard tomography and the corresponding FINT one. For example, given a sinogram in standard tomography with 4096 pixels, ZIF = 4 describes a FINT sinogram consisting of the central 1024 pixels of the original dataset.

Fig. 12[link] shows reconstructions computed by the two iterative methods for increasing ZIFs. The original MLT sinogram has 1500 views × 4096 pixels. At visual inspection the reconstructions are not affected by FINT artefacts and look almost identical, showing that both methods can satisfactorily tackle the reconstruction of FINT data for different ZIFs.

[Figure 12]
Figure 12
Reconstructions computed by ADMP-E and ADMP-V of MLT sinograms with increasing ZIF: 2, 4 and 8. (a) ADMP-E reconstruction: ZIF = 2. (b) ADMP-E reconstruction: ZIF = 4. (c) ADMP-E reconstruction: ZIF = 8. (d) ADMP-V reconstruction: ZIF = 2. (e) ADMP-V reconstruction: ZIF = 4. (f) ADMP-V reconstruction: ZIF = 8.

The reconstructions with ZIF = 8 [Figs. 12(c), 12(f)[link]] are overall characterized by a sort of background pattern. For the simulated data used here, this problem, starting to be barely visible with ZIF = 8, can be neglected up to ZIF ≃ 32. This does not represent a problem for the real datasets used in §4.7[link], as the highest ZIF is roughly 9.

4.5. Validation for different conditions of asymmetry

The accuracy of ADMP-E and ADMP-V has been tested also for the reconstruction of FOVs placed at various distances from the phantom centre. In this way, the iterative methods can be validated for different conditions of asymmetry, i.e. not symmetrical distribution of the attenuation coefficients around the selected FOV.

An example of such a test, performed with the SL phantom, is shown in Fig. 13[link]. The sinograms for FOV-1 and FOV-2 both have 1500 views × 512 pixels. Once again, the reconstructions are not affected by FINT artefacts and look practically identical, proving that both methods can handle FINT cases with pronounced feature asymmetry around the FOV.

[Figure 13]
Figure 13
Reconstructions of two different FOVs computed by ADMP-E and ADMP-V. (a) SL phantom with FOV-1 and FOV-2. (b) ADMP-E reconstruction: FOV-1. (c) ADMP-V reconstruction: FOV-1. (d) ADMP-E reconstruction: FOV-2. (e) ADMP-V reconstruction: FOV-2.

4.6. Reconstruction of underconstrained sinograms

Iterative tomographic algorithms are mostly designed to address underconstrained datasets, providing insufficient direct information for accurate reconstruction with analytical methods. An underconstrained sinogram is either undersampled, noisy or a combination of both factors. A sinogram in standard tomography is considered undersampled if m < [n\,\pi/2], where m is the number of views and n is the number of detector pixels (Kak & Slaney, 2001[Kak, A. C. & Slaney, M. (2001). Principles of Computerized Tomographic Imaging. Philadelphia: Society for Industrial and Applied Mathematics.]). An FINT sinogram is undersampled when m < [n_{\rm{s}}\,\pi/2] with [n_{\rm{s}}] = [n\,{\rm{ZIF}}]: [n_{\rm{s}}] represents the number of detector cells required to `cover' the entire object in standard tomography and n the number of available detector cells from the FINT scan (Xiao et al., 2007[Xiao, X., De Carlo, F. & Stock, S. (2007). Rev. Sci. Instrum. 78, 063705.]). Since low-dose fast scans of FINT tomographic microscopy are usually characterized by [m\ll n\ll n_{\rm{s}}], undersampling combined with local tomography artifacts result in a sort of `background texture', that can severely affect FINT reconstructions. Iterative reconstruction in standard tomography can greatly reduce the amount of this background texture through the usage of constraints and regularization. Since no constraints are available for the FINT datasets, the removal of the background texture relies entirely on the regularization: for the same amount of noise in the data, the regularization strength τ should be bigger, the higher the undersampling factor, defined as UF = [[1-m/(n\,\pi/2\,{\rm{ZIF}})]\times100\%]. If τ is too low, the jump associated with the background texture can be regarded by the TV as a collection of faint edges to be preserved and the regularization will only remove the noise component.

ADMP-E and ADMP-V are here tested for the reconstruction of underconstrained datasets against GRID-E. The regularization strength, τ, is separately chosen for ADMP-E and ADMP-V to guarantee a higher quality reconstruction in terms of CNR (while maintaining a similar spatial resolution) compared with GRID-E.

In the first experiment, a FINT sinogram of the SL phantom with ZIF = 4, 200 views × 512 pixels, UF = 94% and corrupted by Poisson noise with [\sigma] = 2.5% of the sinogram mean is considered. Fig. 14[link] shows that the reconstructions computed by ADMP-E and ADMP-V look more similar to the phantom compared with the GRID-E reconstruction. They have higher MSSIM, PSNR and CNR (increased by a factor of 4.4) scores, as reported in Table 1[link]. The edge profiles (Fig. 15a[link]) indicate that the proposed iterative approaches lead to improved CNR without deterioration of the spatial resolution. Moreover, the peaks corresponding to the different phases of the SL phantom can be clearly identified in the histograms for the iterative reconstructions in Figs. 15(c) and 15(d)[link], whereas only two peaks appear in Fig. 15(b)[link] for the analytical reconstruction.

Table 1
PSNR, MSSIM and CNR computed for the reconstructions of the SL (Fig. 14[link]) phantom

  GRID-E ADMP-E ADMP-V
MSSIM 0.026 0.047 0.045
PSNR 14.74 24.69 24.43
CNR 0.66 2.92 2.85
[Figure 14]
Figure 14
Reconstructions computed by GRID-E (a), ADMP-E (b) and ADMP-V (c) of the SL sinogram with 200 views × 512 pixels, additional Poisson noise with [\sigma] = 2.5% of the sinogram mean and ZIF = 4. The green segment in (a) shows the position of the line profiles displayed in Fig. 15(a)[link]. (a) SL phantom. (b) GRID-E reconstruction. (c) ADMP-E reconstruction. (d) ADMP-V reconstruction.
[Figure 15]
Figure 15
Edge profiles (a) and histograms (b, c, d) for the reconstructions with GRID-E, ADMP-E and ADMP-V shown in Fig. 14[link].

For the second experiment, a FINT sinogram of the MLT phantom with the same condition of undersampling as considered in the previous case (200 views × 512 pixels, ZIF = 4) and a larger amount of Poisson noise (3.5% of the sinogram mean) is used. ADMP-E and ADMP-V substantially reduce the noise in the reconstructions displayed in Figs. 16(c) and 16(d)[link] compared with GRID-E (Fig. 16b[link]) making the light structures clearer against the dark ones. The metric scores in Table 2[link] are largely higher for the iterative reconstructions, with an improvement of the CNR by a factor of nearly 4.4. The edge profiles (Fig. 17a[link]) for the three different reconstructions coincide almost perfectly at the edge position, demonstrating the efficacy of the split Bregman total variation regularization in removing noise while preserving edges. The higher quality of the iterative reconstructions is also highlighted by the histograms in Figs. 17(c) and 17(d)[link] showing clear peaks corresponding to the two phases of the MLT phantom.

Table 2
PSNR, MSSIM and CNR computed for the reconstructions of the MLT (Fig. 16[link]) phantom

  GRID-E ADMP-E ADMP-V
MSSIM 0.101 0.135 0.135
PSNR 8.39 12.68 12.63
CNR 0.44 1.98 1.87
[Figure 16]
Figure 16
Reconstructions computed by GRID-E (a), ADMP-E (b) and ADMP-V (c) of the MLT sinogram with 200 views × 512 pixels, additional Poisson noise with [\sigma] = 3.5% of the sinogram mean and ZIF = 4. The green segment in (a) shows the position of the line profiles displayed in Fig. 17(a)[link]. (a) MLT phantom. (b) GRID-E reconstruction. (c) ADMP-E reconstruction. (d) ADMP-V reconstruction.
[Figure 17]
Figure 17
Edge profiles (a) and histograms (b, c, d) for the reconstructions with GRID-E, ADMP-E and ADMP-V shown in Fig. 16[link].

For the third experiment, the dataset of MOUSE-1, characterized by ZIF ≃ 3.2, 901 views × 2016 pixels, UF = 91%, a pronounced asymmetry and highly absorbing structures (e.g. ribs) outside the FOV, is considered. Due to the feature size and the pattern complexity, the different quality of the reconstructions displayed in Fig. 18[link] can be best appreciated in the blow-ups below the images. Features in the iterative reconstructions can be more easily identified thanks to the reduced noise and to the CNR increased by a factor of 5.3 (Table 3[link]). The line profiles at the edge position (Fig. 19a[link]) for the ADMP-E and ADMP-V results match that for the GRID-E reconstruction. Moreover, the double peak in Figs. 19(c) and 19(d)[link] shows that the two main phases of the lung tissue are better separated in terms of grey level in the iterative reconstructions, whereas a single peak characterizes the histogram for the analytical reconstruction in Fig. 19(b)[link].

Table 3
CNR scores computed for the reconstructions of the sinograms of MOUSE-1 (Fig. 18[link])

  GRID-E ADMP-E ADMP-V
CNR 0.16 0.83 0.85
[Figure 18]
Figure 18
Reconstructions computed by GRID-E (a), ADMP-E (b) and ADMP-V (c) of the MOUSE-1 sinogram with 900 views × 2016 pixels. The green segment in (a) shows the position of the line profiles displayed in Fig. 19(a)[link]. The blow-ups below each reconstruction are zoom-ins for two different ROIs.
[Figure 19]
Figure 19
Edge profiles (a) and histograms (b, c, d) for the reconstructions with GRID-E, ADMP-E and ADMP-V shown in Fig. 18[link].

In the fourth experiment, the reconstructed dataset of MOUSE-2 (901 views × 2016 pixels, ZIF ≃ 3.2) has a very similar pattern complexity as MOUSE-1, whereas the morphology, e.g. the shape of the small structures, is quite different. Reconstructions with the ADMP methods in Fig. 20[link] clearly show higher quality compared with the analytical one, thanks to the decreased noise level and the higher CNR (increased by a factor of 4.6, Table 4[link]). Small dark features are well identifiable in Figs. 20(b) and 20(c)[link], whereas they are mainly covered by noise and undersampling/local tomography artefacts in Fig. 20(a)[link]. The edge profiles in Fig. 21[link] for the iterative reconstructions overlap almost exactly at the edge position with that for the analytical result, indicating also in this case that the total variation regularization operates with negligible loss in terms of spatial resolution. The histograms corresponding to the reconstructions with ADMP-E and ADMP-V in Figs. 21(c) and 21(d)[link] show the presence of two phases. This is not the case in the histogram for GRID-E in Fig. 21(b)[link].

Table 4
CNR scores computed for the reconstructions of the sinograms of MOUSE-2 (Fig. 20[link])

  GRID-E ADMP-E ADMP-V
CNR 0.26 1.21 1.19
[Figure 20]
Figure 20
Reconstructions computed by GRID-E (a), ADMP-E (b) and ADMP-V (c) of the MOUSE-2 sinogram with 900 views × 2016 pixels. The green segment in (a) shows the position of the line profiles displayed in Fig. 21(a)[link]. The blow-ups below each reconstruction are zoom-ins for two different ROIs.
[Figure 21]
Figure 21
Edge profiles (a) and histograms (b, c, d) for the reconstructions with GRID-E, ADMP-E and ADMP-V shown in Fig. 20[link].

For the fifth experiment, the dataset of MOUSE-3, characterized by 500 views × 2016 pixels and ZIF ≃ 9.0, is used. In this case, UF = 98%, but features are on average much larger compared with MOUSE-1 and MOUSE-2, since projections were acquired with a higher magnification. The iterative reconstructions in Fig. 22[link] offer a clearer vision of the object; nevertheless, all structures recognizable in Figs. 22(b) and 22(c)[link] can be visually identified in Fig. 22(a)[link] as well. The CNR is improved in this example by a factor of 3.7 (Table 5[link]). In Fig. 23(a)[link] an edge can hardly be recognized for the reconstruction with GRID-E, whereas those corresponding to ADMP-E and ADMP-V, practically identical, do show a flank. Histograms in Figs. 23(b) and 23(c)[link] show a distinction between the two main phases of the lung tissue, whereas a single peak dominates in the histogram for the analytical reconstruction in Fig. 23(a)[link].

Table 5
CNR scores computed for the reconstructions of the sinograms of MOUSE-3 (Fig. 22[link])

  GRID-E ADMP-E ADMP-V
CNR 0.24 0.88 0.89
[Figure 22]
Figure 22
Reconstructions computed by GRID-E (a), ADMP-E (b) and ADMP-V (c) of the MOUSE-3 sinogram with 500 views × 2016 pixels. The green segment in (a) shows the position of the line profiles displayed in Fig. 23(a)[link]. The blow-ups below each reconstruction are zoom-ins for a ROI.
[Figure 23]
Figure 23
Edge profiles (a) and histograms (b, c, d) for the reconstructions with GRID-E, ADMP-E and ADMP-V shown in Fig. 22[link].

These five experiments with simulated and real sinograms show that: (i) ADMP-E manages to substantially improve the image quality compared with an analytical method like GRID-E when tackling the reconstruction of underconstrained FINT datasets with different ZIF, noise level and feature complexity; (ii) the virtual strategy incorporated in ADMP-V can provide comparable results with those achievable with the ADMP-E.

4.7. Computational cost

To give an idea of the superior computational performance of ADMP-V compared with ADMP-E, the time required for a single iteration, [\Delta t], has been measured for two different datasets on an Intel(R) Core(TM) i7-3520M CPU 2.90 GHz. For a sinogram with 800 views × 504 pixels (convergence reached after eight iterations), [\Delta t](ADMP-E) = 46.7 s and [\Delta t](ADMP-V) = 1.45 s; for a sinogram with 1584 views × 1008 pixels (convergence reached after nine iterations), [\Delta t](ADMP-E) = 173.2 s and [\Delta t](ADMP-V) = 5.7 s. Thus, ADMP-V for small and medium datasets is approximately 30 times faster with respect to ADMP-E.

5. Conclusion

This work introduces two novel strategies for iterative reconstructions of datasets in full interior tomography (FINT), when the FOV is completely inside the object support and no a priori knowledge regarding the support itself and the distribution of the attenuation coefficients in certain ROIs is available. FINT represents a very general case, frequently encountered in many tomographic applications like synchrotron-based X-ray tomographic microscopy.

One strategy works with an edge-padding forward projector. The second is a four-step procedure, requiring the creation of an intermediate virtual sinogram, simulating a full tomography dataset; this sinogram is, then, reconstructed by a standard iterative algorithm, while enforcing a tight constraint outside the reconstruction circle. Both strategies, although not bounded to a specific iterative scheme, have been implemented in this work inside the alternate direction method of multipliers plug-and-play (ADMP), that offers a versatile optimization framework, where different forward projectors and regularizers can be easily plugged in, without altering the structure of the iterative solver. The two resulting iterative methods have been labelled ADMP-E (edge-padding strategy) and ADMP-V (virtual strategy).

The forward gridding projector with minimal oversampling (FGP) is used as standard and edge-padding forward operator for the ADMP. The FGP guarantees fast iterative reconstructions, while keeping the same accuracy of the results achieved with more sophisticated, but much slower, implementations of the Radon transform.

ADMP-E and ADMP-V have been, first, validated for the reconstruction of FINT datasets with different zoom-in factors and asymmetry conditions around the FOV. The methods have, then, been tested on underconstrained simulated and real FINT datasets. Results show that both iterative techniques yield reconstructions of higher quality compared with a standard analytical method: the CNR is greatly improved (on average four times higher), while preserving the spatial resolution, and small features can be more easily identified. The reconstruction quality achieved with the two proposed iterative strategies is comparable. ADMP-V provides, though, superior computational efficiency (about 30 times faster), since it requires a much smaller grid for the computations inside the iterative procedure.

References

First citationAagesen, L. K., Johnson, A. E., Fife, J. L., Voorhees, P. W., Miksis, M. J., Poulsen, S. O., Lauridsen, E. M., Marone, F. & Stampanoni, M. (2010). Nat. Phys. 6, 796–800.  Web of Science CrossRef CAS Google Scholar
First citationArcadu, F., Nilchian, M., Studer, A., Stampanoni, M. & Marone, F. (2016a). IEEE Trans. Image Process. 25, 1207–1218.  Web of Science CrossRef Google Scholar
First citationArcadu, F., Stampanoni, M. & Marone, F. (2016b). Opt. Express. Submitted.  Google Scholar
First citationBerg, S., Ott, H., Klapp, S. A., Schwing, A., Neiteler, R., Brussee, N., Makurat, A., Leu, L., Enzmann, F., Schwarz, J. O., Ersten, M. K., Irvine, S. & Stampanoni, M. (2012). Proc. Natl Acad. Sci. 110, 3755–3759.  Web of Science CrossRef Google Scholar
First citationBoyd, S., Parikh, N., Chu, E., Peleato, B. & Eckstein, J. (2010). Mach. Learn. 3, 1–122.  CrossRef Google Scholar
First citationCampi, G., Bianconi, A., Poccia, N., Bianconi, G., Barba, L., Arrighetti, G., Innocenti, D., Karpinski, J., Zhigadlo, N. D., Kazakov, S. M., Burghammer, M., Zimmermann, M., v, Sprung, M. & Ricci, A. (2015). Nature (London), 525, 359–362.  Web of Science CrossRef CAS Google Scholar
First citationCourdurier, M., Noo, F., Defrise, M. & Kudo, H. (2008). Inverse Probl. 24, 65001.  Web of Science CrossRef Google Scholar
First citationDefrise, M., Noo, F., Clackdoyle, R. & Kudo, H. (2006). Inverse Probl. 22, 1037–1053.  Web of Science CrossRef Google Scholar
First citationEbner, M., Marone, F., Stampanoni, M. & Wood, V. (2013). Science, 342, 716–720.  Web of Science CrossRef CAS PubMed Google Scholar
First citationGanne, J., De Andrade, V., Weinberg, R. F., Vidal, O., Dubacq, B., Kagambega, N., Naba, S., Baratoux, L., Jessell, M. & Allibon, J. (2012). Nat. Geosci. 5, 60–65.  Web of Science CrossRef CAS Google Scholar
First citationGoldstein, T. & Osher, S. (2009). SIAM J. Imaging Sci. 2, 323–343.  Web of Science CrossRef Google Scholar
First citationHesse, B., Varga, P., Langer, M., Pacureanu, A., Schrof, S., Männicke, N., Suhonen, H., Maurer, P., Cloetens, P., Peyrin, F. & Raum, K. (2015). J. Bone Miner. Res. 30, 346–356.  Web of Science CrossRef CAS PubMed Google Scholar
First citationHestenes, M. R. & Stiefel, E. (1952). J. Res. Natl. Bur. Stan. 49, 409–436.  CrossRef Web of Science Google Scholar
First citationHuyn-Thu, Q. & Ghanbaru, M. (2008). Electron. Lett. 44, 800–801.  Google Scholar
First citationKak, A. C. & Slaney, M. (2001). Principles of Computerized Tomographic Imaging. Philadelphia: Society for Industrial and Applied Mathematics.  Google Scholar
First citationKudo, H., Courdurier, M., Noo, F. & Defrise, M. (2008). Phys. Med. Biol. 53, 2207–2231.  Web of Science CrossRef PubMed Google Scholar
First citationKyrieleis, A. V., Titarenko, V., Ibison, M., Connolley, T. & Withers, P. J. (2011). J. Microsc. 241, 69–82.  Web of Science CrossRef CAS Google Scholar
First citationLovric, G., Barré, S. F., Schittny, J. C., Roth-Kleiner, M., Stampanoni, M. & Mokso, R. (2013). J. Appl. Cryst. 46, 856–860.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationLu, J., Lee, Y. J., Luo, X., Lau, K. C., Asadi, M., Wang, H. H., Brombosz, S., Wen, J., Zhai, D., Chen, Z., Miller, D. J., Jeong, Y. S., Park, J. B., Fang, Z. Z., Kumar, B., Salehi-Khojin, A., Sun, Y. K., Curtiss, L. A. & Amine, K. (2016). Nature (London), 529, 377–382.  Web of Science CrossRef CAS Google Scholar
First citationMarone, F., Münch, B. & Stampanoni, M. (2010). Proc. SPIE, 7804, 780410.  CrossRef Google Scholar
First citationNatterer, F. & Wübbeling, F. (2001). Mathematical Methods in Image Reconstruction. Philadelphia: Society for Industrial and Applied Mathematics.  Google Scholar
First citationNoo, F., Clackdoyle, R. & Pack, J. D. (2004). Phys. Med. Biol. 49, 3903–3923.  Web of Science CrossRef PubMed Google Scholar
First citationPaganin, D., Mayo, S., Gureyev, T., Miller, P. & Wilkins, S. W. (2002). J. Microsc. 206, 33–40.  Web of Science CrossRef PubMed CAS Google Scholar
First citationRashed, E. A. & Kudo, H. (2013). J. Synchrotron Rad. 20, 116–124.  Web of Science CrossRef CAS IUCr Journals Google Scholar
First citationSavitzky, A. & Golay, M. J. E. (1964). Anal. Chem. 36, 1627–1639.  CrossRef CAS Web of Science Google Scholar
First citationSeger, M. M. (2002). Proceeding of the Symposium on Image Analysis (SSAB02), Lund, Sweden, pp. 33–36.  Google Scholar
First citationShepp, L. & Logan, B. F. (1974). IEEE Trans. Nucl. Sci. 21, 2143.  Google Scholar
First citationToft, P. (1996). PhD thesis, Denmark Technical University, Denmark.  Google Scholar
First citationVenkatakrishnan, S. V., Bouman, C. A. & Wohlberg, B. (2013). Proceeding of the 2013 IEEE Global Conference on Signal and Information Processing (Globalsip 2013), Austin, TX, USA, 3–5 December 2013, pp. 945–948.  Google Scholar
First citationWalker, S. M., Schwyn, D. A., Mokso, R., Wicklein, M., Müller, T., Doube, M. & Stampanoni, M. (2014). PLoS Biol. 12, e1001823.  Web of Science CrossRef Google Scholar
First citationWang, Z., Bovik, A. C., Sheikh, H. R. & Simoncelli, E. P. (2004). IEEE Trans. Image Process. 13, 600–612.  Web of Science CrossRef PubMed Google Scholar
First citationXiao, X., De Carlo, F. & Stock, S. (2007). Rev. Sci. Instrum. 78, 063705.  Web of Science CrossRef Google Scholar
First citationXu, Q., Mou, X., Wang, G., Sieren, J., Hoffman, E. A. & Yu, H. (2011). IEEE Trans. Med. Imaging, 30, 1116–1128.  Web of Science CrossRef Google Scholar
First citationYe, Y., Yu, H., Wei, Y. & Wang, G. (2007). Intl J. Biomed. Imaging, 2007, 63634.  Google Scholar
First citationYu, H. & Wang, G. (2009). Phys. Med. Biol. 54, 2791–2805.  Web of Science CrossRef PubMed Google Scholar
First citationZhang, H., Huang, J., Ma, J., Bian, Z., Feng, Q., Lu, H., Liang, Z. & Chen, W. (2014). IEEE Trans. Biomed. Eng. 61, 2367–2378.  Web of Science CrossRef Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.

Journal logoJOURNAL OF
SYNCHROTRON
RADIATION
ISSN: 1600-5775
Follow J. Synchrotron Rad.
Sign up for e-alerts
Follow J. Synchrotron Rad. on Twitter
Follow us on facebook
Sign up for RSS feeds