Open Access
Add to BibSonomyAbstract
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 Optical Society of America
1. Introduction
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 [1]. 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 [2]. 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 [3], individual projections of the object must be interrelated leading to a certain degree of redundancy within the tomographic data. The theory of Compressed Sensing [4] 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 [5]. 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 [6]. 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 [7]. 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 λ:
Here w stands for the spatial frequency and δ(w) denotes the delta function at the origin of spatial frequency coordinates. A parallel monochromatic X-ray beam with a uniform illumination and intensity I0 = 1 was assumed for simplicity. This model is valid for chemically homogeneous or quasi-homogeneous objects (i.e. ϕ/μ≈ const) with weak attenuation and slow-varying phase or for objects that are composed from light elements in a limited range of X-ray energies [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.
Fig. 1 Central slice theorem for single-distance X-ray Phase Retrieval Tomography. Spatial domain representation: observations I(θ, s) can be modeled as the projection of the unknown image f(x, y) followed by the convolution with a linear propagator. Fourier domain representation: Fourier transform of the observed image Ĩ(θ, w) can be modeled as the slice of the unknown image f̃(wcosθ, wsinθ)) multiplied with P(w). The Fourier representation of the unknown image f (u, v) is irrecoverable at the frequency bands corresponding to zero-crossings of P(w).
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:
When the matrix notation is considered, 𝒫 can be represented by a square k × k diagonal matrix: where P(w) = 2σsin(πλDw2) − 2cos(πλDw2). A more detailed explanation of the matrix notation is given in our paper on two-dimensional phase retrieval [2]. Now the approximate solution to Eq. (1) can be found by solving the corresponding least squares problem: Here represents the L2 norm of the corresponding functional. To avoid amplification of errors by the terms that correspond to small P(w), we will exclude these terms from the system of equations. To do so we assume that the image μ is irrecoverable within the frequency bands |w − w0| < ε, where w0 is the nearest zero-crossing of the CTF and ε is a small constant that depends on the signal-to-noise ratio. Now an additional linear operator can be introduced: Operator 𝒵 can be represented by a binary matrix and is applied to both sides of Eq. (2), making sure that terms corresponding to the small P(w) are set to zero: The system that we have obtained is underdetermined (it does not allow to determine (ℱμ) for |w − w0| < ε) and does not have a unique solution unless additional constraints are added. Hofmann [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 [10].2.2. Tomography
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:
According to the central slice theorem, if a two-dimensional Fourier transform of the image f(x, y) is computed values of f̃(u, v) that lie on a radial line passing through the center of coordinates under an angle θ (central slice) will correspond to the one-dimensional Fourier transform taken along the s coordinate of the projection p(θ, s): The central slice theorem demonstrates an important relation between the Radon transform and the Fourier transform of the two-dimensional image f(x, y). It follows from 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 [4] 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 tomography
By combining (1) and (9) we can rewrite the central slice theorem for what we call Phase Retrieval Tomography (PRT), i.e. phase retrieval combined with tomographic reconstruction:
Here Ĩ(θ, w) is the one-dimensional Fourier transform of the phase-contrast sinogram. According to Eq. (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):
Here ℛ can be represented by an n × m matrix. This system can be solved using either an approximation of the ℛ−1 (e.g. filtered back-projection) or in least-squares sense using one of the algebraic reconstruction algorithms (e.g. EM, ART, SIRT) [11].For PRT, the tomography model represented by Eq. (11) can be combined with the linear phase-contrast model of Eq. (6) into a single linear inverse problem:
Here, the linear operators 𝒵, 𝒫 and ℱ are applied to vectors of n elements that contain all projections generated by ℛf instead of a single projection of k elements as it is in 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.2.4. Preconditioning
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.
We will assume that for the frequencies |w − w0| > ε, Eq. (1) can be accurately computed using direct inversion. Then Eq. (12) can be rewritten as follows:
where μ is a vector representation of the projected attenuation of the object, which can be calculated from the phase-contrast sinogram I using direct phase retrieval. Algebraic methods based on 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:
where ℛ−1 is the approximated discrete inverse Radon transform (e.g. filtered back-projection), ℱ̂ is the two-dimensional discrete Fourier transform and 𝒫̂ represents an elementwise multiplication with the discrete version of the two-dimensional CTF = 2σsin(πλD(u2 + v2)) −2cos(πλD(u2 + v2)). 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:
where 𝒜 = 𝒵𝒫ℱℛ and Ĩ = 𝒵ℱI. Using gradient descent, f is computed according to the following iterative scheme: where j is the iteration number, 𝒜T is the transpose of 𝒜, so 𝒜T = ℛTℱT𝒫T𝒵T. The constant α represents the step size in the opposite gradient direction. In order to guarantee convergence, the constant α has to be sufficiently small and can either be calculated using an additional line-search step or from the eigenvalues of (𝒜T𝒜) [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:
where , ∇h and ∇v are the horizontal and vertical finite difference operators and λTV denotes the weight of the regularization term. Depending on the magnitude of λTV the solution will be promoted either towards greater conformity with the observed data or towards greater sparsity of the gradient magnitude. 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 [15]. The same algorithm was applied by us to several different phase-retrieval problems in [2]. 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.3. Simulations
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:
were δ denotes the decrement of the complex refractive index, β stands for the attenuation index and f(x, y) denotes a dimensionless normalized attenuation of the digital phantom. Subsequently, the intensity images were generated using the following expression: where the Fresnel propagator is represented by Pλ = eiπλDw2 and the Optical Transfer Function (OTF) of the imaging system is modeled by a Gaussian function with standard deviation a (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)).Fig. 2 Simulated phase-contrast tomography of the Shepp-Logan phantom. (a) Ground truth; (b) projected attenuation; (c) observed phase-contrast sinogram; (d) sinogram after phase-retrieval; (e) ’direct’ filtered-back projection; (f) filtered-back projection after phase retrieval.
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.
Table 1. Simulation parameter sets
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:
The weight of the regularization term λTV had to be adjusted depending on the underlying linear system and the variance of the simulated data. Automated methods for determining the regularization parameters produce a wide spread of results depending on the method [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.Fig. 3 Reconstructions of the simulated data (error magnitude). Columns correspond to the reconstruction methods: sequential, unconstrained algebraic, full algebraic, algebraic tomographic reconstruction and algebraic phase retrieval. Rows correspond to the different simulations: weak phase, few projections, blur, strong phase, noise, realistic.
Table 2. RMSE for six different simulations (rows) and five reconstruction algorithms (columns).
4. Experiments
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.
Fig. 4 X-ray PCT reconstruction of the glass capillary filled with copper spheres (slice from the middle of the volume). Reconstructions based on 650 projections: (a) - sequential approach based on L2-regularized phase retrieval, (c) - full algebraic reconstruction, (e) - algebraic phase retrieval. Reconstructions based on 65 projections: (b) - sequential approach, (d) - full algebraic reconstruction, (f) - algebraic phase retrieval. Red line shows position of the profiles depicted on the next figure.
Fig. 5 Attenuation profiles of the copper sphere and the quartz wall for different reconstruction algorithms. (a) Reconstruction based on 650 projections; (b) reconstruction based on 65 projections.
5. Conclusion
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 [6]. 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 [2]. Specifically, the problem of tomographic reconstruction based on a multi-distance CTF model [19] remains linear, whereas incorporation of the so called mixed phase retrieval [20] 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.
Acknowledgments
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] .
9. X. Wu and A. Yan, “Phase retrieval from one single phase-contrast x-ray image,” Opt. Express 17, 11187 (2009) [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] .
