The reconstruction problem in in-line X-ray Phase-Contrast Tomography is usually approached by solving two independent linearized sub-problems: phase retrieval and tomographic reconstruction. Both problems are often ill-posed and require the use of regularization techniques that lead to artifacts in the reconstructed image. We present a novel reconstruction approach that solves two coupled linear problems algebraically. Our approach is based on the assumption that the frequency space of the tomogram can be divided into bands that are accurately recovered and bands that are undefined by the observations. This results in an underdetermined linear system of equations. We investigate how this system can be solved using three different algebraic reconstruction algorithms based on Total Variation minimization. These algorithms are compared using both simulated and experimental data. Our results demonstrate that in many cases the proposed algebraic algorithms yield a significantly improved accuracy over the conventional L2-regularized closed-form solution. This work demonstrates that algebraic algorithms may become an important tool in applications where the acquisition time and the delivered radiation dose must be minimized.
© 2013 OSA
Quantitative X-ray Phase-Contrast Tomography (X-ray PCT) requires three-dimensional image reconstruction from a series of two-dimensional in-line phase-contrast images acquired under various angles. Within the linear approximation, the problem of image reconstruction in X-ray PCT is typically treated as two separate linear subproblems which are solved sequentially . Namely, the phase retrieval problem, where the projected refraction index (i.e. linear phase and attenuation) of the specimen is retrieved independently for each recorded phase-contrast image and the problem of tomographic reconstruction, where the three-dimensional distribution of the refraction index is computed from the collection of retrieved projections.
Most often linear phase retrieval of an individual tomographic projection is associated with inversion of an ill-posed linear system. In such cases regularization (e.q. L2 or L1 regularization) permits computation of an approximate inverse solution . However, it will often cause artifacts that are subsequently propagated into the tomographic reconstruction of the object.
We believe that a more accurate reconstruction approach can be designed by combining phase retrieval and tomographic reconstruction into a single underdetermined linear problem. This assumption is based on the fact that tomographic projections of the object are in general not independent from each other. According to the so called Helgason-Ludwig consistency conditions , individual projections of the object must be interrelated leading to a certain degree of redundancy within the tomographic data. The theory of Compressed Sensing  suggests an even stronger correlation between the individual tomographic projections for objects that have a sparse representation in a certain domain. For sparse solutions, the projection data shows a low rate of innovation, which implies that a limited number of projections suffices for exact reconstruction . These facts lead us to a conclusion that the accuracy of phase retrieval of a single tomographic projection can benefit from the redundancy contained in the complete tomographic dataset. Regularization techniques based on the assumption of sparsity were shown to improve the quality of image reconstruction in cases when the reconstructed image is not strictly sparse . We will demonstrate how such redundancy can be exploited using an iterative algebraic reconstruction approach for in-line X-ray PCT.
2. Materials and methods
2.1. Single-distance phase retrieval
Let us consider a single-distance phase-retrieval model based on the Contrast Transfer Function (CTF) approach . In this model the projected phase image of the specimen ϕ(s) is assumed to be proportional to the projected attenuation image μ(s). The proportionality ratio σ = ϕ/μ can be calculated as the ratio between the real and imaginary parts of the refractive index of the specimen. This allows us to express the Fourier transform of the phase-contrast image Ĩ(w) as the Fourier transform of the projected attenuation image μ̃(w) multiplied with the CTF that corresponds to the object-to-detector distance D and the X-ray wavelength λ:8, 9]. Obviously, it is impossible to recover the projected attenuation image μ̃(w) from the observed phase-contrast image Ĩ(w) at frequencies w that correspond to the zero-crossings of the CTF (see Fig. 1). Also, in the vicinity of each zero-crossing of the CTF the inverse problem based on Eq. (1) will be ill-posed due to measurement noise in addition to systematic linearization errors. At these frequencies μ̃(w) has to be computed using additional constraints that favor a particular solution based on a-priori knowledge.
Let us assume that the vectors μ and I are defined in the spatial domain on a uniform grid of k values that belong to μ(s) and (I(s) − 1) respectively. A linear operator ℱ will represent the uniform discrete Fourier transform. Element-wise multiplication with the CTF function is denoted by the linear operator 𝒫. Then Eq. (1) can be discretized in the following form:2]. Now the approximate solution to Eq. (1) can be found by solving the corresponding least squares problem: Eq. (2), making sure that terms corresponding to the small P(w) are set to zero: 10] proposed to use a simple constraint (ℱμ) = 0 for |w− w0| < ε in order to solve the phase-retrieval problem. However, we expect that applying sparsity constraints to the joint phase-contrast-tomographic reconstruction will yield a better accuracy. A frequency dependent ε can be defined when it is possible to estimate the power spectrum of the noise and the power spectrum of the reconstructed image modulated with the CTF a-priori. Otherwise a constant factor ε can be chosen, for instance, after evaluating the quality of the direct reconstruction proposed in .
In this subsection we will use coordinate s to describe the transverse coordinate an observed image of the object and w as its counterpart in the Fourier domain. A new coordinate θ will be introduced for the angle at which the image is recorded during the tomographic acquisition. Coordinates (x, y) will be used to describe a 2D tomographic image of the object (see Fig. 1).
Consider a two-dimensional image f(x, y) which is defined on ℝ2 → ℂ as a square integrable function with bounded support. A linear projection of f(x, y) in the direction θ will be defined along the coordinate s as:Eq. (9) that the function f(x, y) can be sampled in 2D Fourier space using 1D Fourier transforms of its own projections p(θ, s). This facilitates direct reconstruction based on the inverse Fourier transform of f̃(u, v) which has to be computed using interpolation. Moreover, it was shown in  that by using appropriate additional constraints, a piece-wise constant object (i.e. an object that has sparse gradient magnitude) can be accurately reconstructed from a severely undersampled Fourier representation f̃(u, v), i.e. using very few projections. The latter has important consequences for regularization of the phase-retrieval problem described in the following subsection.
2.3. Phase-contrast tomographyEq. (10), f̃(u, v) is undetermined for frequencies that correspond to zero-crossings of the CTF. And since the accuracy of the linear approximation and the observations Ĩ(θ, w) is finite, f̃(u, v) can only be computed outside of the circular bands that correspond to |w− w0|< ε (Fig. 1). Thus, the described problem of the tomographic reconstruction based on single-distance phase-contrast data is underdetermined.
Let us introduce a discrete space-domain representation of the tomography problem. Assume that the vector f is composed from m samples of f (x, y) defined on a Cartesian grid. Projections p(θ, s) are sampled on a regular grid of n elements stored in a vector p. Using the Radon transform operator ℛ we can write a discrete representation of Eq. (7):11].Eq. (6). Their matrix representations can be easily constructed by stacking the corresponding “single projection” matrices into a block diagonal matrix of n × n elements.
Before introducing algebraic methods suitable for solving Eq. (12), we would like to introduce two other representations of the PCT problem that, under certain conditions, may be preferable due to faster convergence.Eq. (13) are likely to converge faster than the ones based on Eq. (12), since the inversion of the 𝒫 operator is already solved during the phase retrieval step. However, unlike the approach, where the tomographic reconstruction is computed subsequently after phase retrieval, Eq. (13) allows us to take into account that the information about the attenuation of the object is lost at spatial frequencies |w − w0| < ε. Such an approach can also be used when the number of tomographic projections is limited.
Yet another version of the linear system can be derived from the central slice theorem given by Eq. (10). Now for the case that an accurate inverse Radon transform of the sinogram can be calculated. After changing from polar coordinates to Cartesian coordinates in Eq. (10), we can write down the following discretized system:Equation (14) permits application of algebraic algorithms with additional constraints to the phase retrieval problem while the problem of tomographic reconstruction is solved in a non-iterative manner. This approach can also be faster than calculation based on Eq. (12), since it does not require recalculation of the back and forward Radon transform for each iteration.
2.5. Algebraic methods
An approximate solution of Eq. (12) can be found using least-squares minimization:12]. It is worth noting that the Fourier transform operator ℱ is orthogonal, so ℱT = ℱ−1; operators 𝒵 and 𝒫 are represented by diagonal matrices, so 𝒵T = 𝒵, 𝒫T = 𝒫, and ℛT represents the unfiltered back-projection operator.
We have mentioned previously that in conventional tomography a piecewise constant image f can often be computed accurately from severely underdetermined tomographic system using methods with additional constraints, such as TV minimization. In the TV minimization approach an additional term is added to the objective function:Equation (17) represents a non-smooth convex minimization problem which can be solved using one of the iterative TV minimization methods [13, 14]. In the current investigation we will demonstrate results obtained using so called FISTA-based TV minimization which is described in detail in . The same algorithm was applied by us to several different phase-retrieval problems in . Further on we call the solutions based on minimization of Eq. (12) - full algebraic reconstruction, the ones based on Eq. (13) - algebraic tomographic reconstruction and Eq. (14) - algebraic phase retrieval.
In order to compare the accuracy of tomographic reconstructions based on the proposed reconstruction approaches, we simulated phase-contrast tomographic data using Fresnel propagation. A discretized Shepp-Logan phantom (256 × 256 pixels) was used as a ground truth image to generate the density image f (x, y) (Fig. 2(a)). The wave function of the object T(θ, s) was computed for each tomographic angle θ using the following expression:Fig. 2(c)). Figure 2(e) shows the result of direct application of filtered-back projection to the raw phase-contrast images. The result of the sequential reconstruction is shown on Fig. 2(d) and Fig. 2(f). The retrieved phase images (Fig. 2(d)) suffers from typical artifacts –low frequency noise and fringes around large variations of intensity. These artifacts correspond to spatial frequencies at which the CTF has low amplitude and the reconstruction errors are high. They are propagated into the subsequent tomographic reconstruction (Fig. 2(f)).
All images were generated for monochromatic homogeneous illumination with a wavelength λ = 0.31Å (energy 40 KeV), propagation length D = 1 m and a pixel size of 1 μm. Noise was added to the sinogram after modeling the complete imaging chain. Since the intensity is varying only moderately across the sinogram, an additive Gaussian noise term was assumed to be a good model. Parameters used in various simulations are shown in Table 1.
Simulated data was treated by five different reconstruction techniques: sequential approach (phase retrieval followed by the filtered back-projection), unconstrained algebraic reconstruction based on Eq. (12) and three algorithms based on TV minimization: full algebraic reconstruction - Eq. (12), algebraic tomographic reconstruction - Eq. (13), algebraic phase retrieval - Eq. (14). A known Gaussian OTF was included in all algebraic reconstruction models in order to achieve a sharper reconstruction. Reconstructions were performed on a 512 × 512 pixels grid in order to avoid boundary effects. The same stopping condition was used for all algebraic methods:16]. In this work the regularization parameters were selected empirically so, for a given combination of the linear model and the simulation conditions, a solution with a small RMSE would be obtained. For full algebraic reconstruction λTV varied in the range from 10−7 (for the “blur” simulation) to 10−5 (for the “strong phase” and “realistic” simulations). For algebraic tomographic reconstruction λTV was in the range from 10−10 to 10−8. And for algebraic phase retrieval we used λTV ranging from 10−4 to 10−2. The constant parameter ε was set to 104m−1 for simulations with small errors (“weak phase”, “few projections” and “blur”) and 105m−1 for simulations with large errors (“strong phase”, “noise” and “realistic”). Figure 3 illustrates the error magnitude associated to the resulting reconstructions. The corresponding Root Mean Square Error (RMSE) can be found in Table 2.
In this section we are presenting some preliminary results of the full algebraic reconstruction and the algebraic phase retrieval methods applied to an experimental X-ray PCT dataset. The data was collected at the beamline ID11 of the European Synchrotron Radiation Facility (Grenoble, France). The tomographic scan was acquired in-situ for spherical polycrystalline copper (Makin Metal Powders (UK) Ltd., diameter 50 μm) during sintering at 1050°C. The specimen was placed in a quartz capillary with a 500 μm internal diameter. During the experiment gas shielding (argon: 98% and hydrogen: 2%) was applied. The scan was performed in a continuous 180° mode with 650 projections. Phase-contrast images were acquired using an X-ray beam with a mean energy E = 40 KeV (ΔE/E = 10−3), a source to object distance of 96 m and a propagation distance of 25 cm. The size of each image was 512 × 256 pixels with a pixel size of 1.4 × 1.4μm2.
We compared the standard sequential reconstruction employing L2-regularization with the proposed algebraic methods based on TV-minimization. Figure 4 shows the results obtained using three different reconstruction techniques: sequential approach, full algebraic reconstruction, and algebraic phase retrieval. All reconstruction techniques were applied to the complete dataset of 650 projections and a subset of only 65 projections. The red line on Fig. 4 shows the position of the attenuation profiles that are depicted in Fig. 5. In the sequential approach L2-regularization was used during the phase-retrieval, while the algebraic phase retrieval and the full algebraic reconstruction were performed using a three-dimensional FISTA-based TV minimization with non-negativity constraint. A Gaussian OTF was added to the CTF model in order to account for blurring with FWHM of 4.95 μm. We varied the number of iterations depending on the rough estimate of the convergence speed of a particular algorithm. In full algebraic reconstruction we used 2000 iterations. In algebraic phase retrieval based on 650 tomographic projections 300 iterations were used. Two times more iterations were used in both approaches when applied to 65 tomographic projections.
The results of reconstructions for the simulated data (Fig. 3) demonstrate that a TV minimization approach can yield a nearly flawless tomographic reconstruction based on a single distance X-ray PCT data (full algebraic reconstruction applied to the weak phase case). In most of the demonstrated examples the full algebraic reconstruction approach outperforms the other two approaches: the algebraic tomographic reconstruction and the algebraic phase retrieval.
However, there are certain cases where the technique will not work so well. The algebraic tomographic reconstruction method clearly fails in the case with strong blur applied to the simulated data. That is an expected outcome since the blur is not included in the linear system that is solved algebraically. The algebraic phase retrieval method is failing in the case that the dataset contained only a few projections. It is also expected since the back-projection sub-problem is not solved by the algebraic algorithm but calculated once using filtered back-projection. However, given the advantages of the full algebraic reconstruction, it also suffers from certain disadvantages - it is the most computation-intensive of the methods proposed in this article. It also is likely to converge slower than the other two in terms of the number of iterations.
All reconstructions in this paper were applied to the data (recorded as well as simulated) of objects that comply with the assumption of sparsity, i.e. objects with a piece-wise constant attenuation. Other studies show that TV minimization can improve the accuracy of the reconstruction of images that are not strictly sparse . Further investigation should be carried out in order to test the applicability of the Phase Retrieval Tomography based on TV minimization to non-sparse objects, such as those encountered in soft-tissue imaging.
In the current study we investigated a simple case of tomographic reconstruction for parallel beam geometry combined with a single-distance CTF phase-retrieval model. This combination represents a linear inversion problem, so algebraic reconstruction algorithms can be applied to solve it with minimal adjustments. However, there are a large number of variations to the proposed technique that can be considered during further investigation. A generalized description of the central slice theorem for fan-beam or cone-beam geometries [17, 18] should allow extension of the current approach to these geometries. Some of the phase retrieval models other than the CTF model can be easily incorporated into an algebraic reconstruction similar to the two-dimensional phase retrieval in . Specifically, the problem of tomographic reconstruction based on a multi-distance CTF model  remains linear, whereas incorporation of the so called mixed phase retrieval  will require solving a non-linear problem.
The experimental data was recorded during an in-situ sintering experiment during which acquisition at shorter propagation distance or in attenuation-contrast mode was not possible. Taking into account that the specimen yields strong attenuation and rapid phase variations, this imaging regime leads to significant artifacts when linear phase retrieval or no phase retrieval is used for tomographic reconstruction. However, it seems that the use of TV minimization with a non-negativity constraint leads to a solution with visibly higher contrast and sharper boundaries. A further development of algebraic techniques may facilitate more accurate tomographic reconstructions based on experimental data acquired under similar (suboptimal) conditions.
Another important result is the full algebraic reconstruction based on fewer projections. It can be seen, that this method allows to reduce the artifacts significantly in reconstruction based on a limited number of projections. Both results can justify the high computational cost of the proposed algebraic algorithms in applications where the acquisition time and the radiation dose are highly restricted.
This work was partially supported by the Care4U project with financial support of Point One, an innovation program of the Ministry of Economic Affairs in The Netherlands. K. Joost Batenburg was financially supported by the NWO (the Netherlands Organisation for Scientific Research - The Netherlands, research programme 639.072.005).
References and links
1. R. C. Chen, L. Rigon, and R. Longo, “Quantitative 3D refractive index decrement reconstruction using single-distance phase-contrast tomography data,” J. Phys. D Appl. Phys. 44, 9 (2011) [CrossRef] .
2. A. Kostenko, K.J. Batenburg, H. Suhonen, S.E. Offerman, and L.J. van Vliet, “Phase retrieval in in-line x-ray phase-contrast imaging based on total variation minimization,” Opt. Expr. 21, 710–723 (2013) [CrossRef] .
3. F. Natterer, The Mathematics of Computerized Tomography( New York: Wiley, 1986).
4. E. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory 52, 489–509 (2006) [CrossRef] .
5. M. Vetterli, P. Marziliano, and T. Blu, “Sampling signals with finite rate of innovation,” IEEE Trans. Signal Proc. 50, 1417–1428 (2003) [CrossRef] .
6. X. Bresson and T. F. Chan, “Fast dual minimization of the vectorial total variation norm and applications to color image processing,” Inv. Probl. and Imaging 2, 455–484 (2008) [CrossRef] .
7. L. Turner, B. Dhal, J. Hayes, A. Mancuso, K. Nugent, D. Paterson, R. Scholten, C. Tran, and A. Peele, “X-ray phase imaging: demonstration of extended conditions with homogeneous objects,” Opt. Express 12, 2960–2965 (2004) [CrossRef] [PubMed] .
8. D. Paganin, S. Mayo, T. Gureyev, P. Miller, and S. Wilkins, “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” J. Microsc. 206, 33–40 (2002) [CrossRef] [PubMed] .
10. R. Hofmann, J. Moosmann, and T. Baumbach, “Criticality in single-distance phase retrieval,” Opt. Express 19, 25881–25890 (2011) [CrossRef] .
11. A. C. Kak and M. Slaney, Principles of computerized tomographic imaging (IEEE Press, 1988).
12. L. Armijo, “Minimization of functions having Lipschitz continuous first partial derivatives,” Pacific J. Math. 16, 1–3 (1966) [CrossRef] .
13. A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vision 20, 89–97 (2004) [CrossRef] .
14. J. Dahl, P. C. Hansen, S. H. Jensen, and T. L. Jensen, “Algorithms and software for total variation image reconstruction via first-order methods,” Num. Alg. 53, 67–92 (2010) [CrossRef] .
15. A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Im. Sci. 2, 183–202 (2009) [CrossRef] .
16. G. M. P. van Kemplen and L. J. van Vliet, “The influence of the regularization parameter and the first estimate on the performance of Tikhonov regularized non-linear image restoration algorithms,” J. of Microscopy 198, 63–75 (2000) [CrossRef] .
17. S.-R. Zhao and H. Halling, “A new Fourier method for fan beam reconstruction,” IEEE Nucl. Sci. Symp. Med. and Imaging Conf. 2, 1287–1291 (1995).
18. G.-H. Chen, S. Leng, and C. A. Mistretta, “A novel extension of the parallel-beam projection-slice theorem to divergent fan-beam and cone-beam projections,” Med. Phys. 32, 654–665 (2005) [CrossRef] [PubMed] .
19. P. Cloetens, W. Ludwig, J. Baruchel, D. van Dyck, J. van Landuyt, J. P. Guigay, and M. Schlenker, “Holotomography: Quantitative phase tomography with micrometer resolution using hard synchrotron radiation x-rays,” Appl. Phys. Lett. 75, 2912–2914 (1999) [CrossRef] .
20. J. P. Guigay, M. Langer, R. Boistel, and P. Cloetens, “A mixed contrast transfer and transport of intensity approach for phase retrieval in the Fresnel region,” Opt. Lett. 32, 1617–1619 (2007) [CrossRef] [PubMed] .