Phase retrieval in in-line x-ray phase contrast imaging based on total variation minimization

State-of-the-artechniques for phase retrieval in propagation based X-ray phase-contrast imaging are aiming to solve an underdetermined linear system of equations. They commonly employ Tikhonov regularization − an L2-norm regularized deconvolution scheme − despite some of its limitations. We present a novel approach to phase retrieval based on Total Variation (TV) minimization. We incorporated TV minimization for deconvolution in phase retrieval using a variety of the most common linear phase-contrast models. The results of our TV minimization was compared with Tikhonov regularized deconvolution on simulated as well as experimental data. The presented method was shown to deliver improved accuracy in reconstructions based on a single distance as well as multiple distance phase-contrast images corrupted by noise and hampered by errors due to nonlinear imaging effects. © 2013 Optical Society of America OCIS codes:(110.7440) X-ray imaging; (100.5070) Phase retrieval; (100.3190) Inverse problems. References and links 1. O. Coindreau, C. Mulat, C. Germain, J. Lachaud, and G. L. Vignoles, “Benefits of x-ray CMT for the modeling of C/C composites,” Adv. Eng. Mat. 13,178–185 (2011). 2. F. Cosmi, A. Bernasconi, and N. Sodini, “Phase contrast micro-tomography and morphological analysis of a short carbon fibre reinforced polyamide,” Compos. Sci. Technol. 71,23–30 (2011). 3. M. Herbig, A. King, P. Reischig, H. Proudhon, E. M. Lauridsen, J. Marrow, J. -Y. Buffire, and W. Ludwig, “3-D growth of a short fatigue crack within a polycrystalline microstructure studied using combined diffraction and phase-contrast x-ray tomography,” Acta Mater. 59,590–601 (2011). 4. A. Kostenko, H. Sharma, E. G. Dere, A. King, W. Ludwig, W. van Oel, S. Stallinga, L. J. van Vliet, and S. E. Offerman, “Three-dimensional morphology of cementite in steel studied by x-ray phase-contrast tomography,” Scr. Mater.67, 261-264 (2012). 5. T. Argunova, M. Gutkin, J. H. Je, E. Mokhov, S. Nagalyuk, S. Sergey, and Y. Hwu, “SR phase-contrast imaging to address the evolution of defects during SiC growth,” Phys. Status Solid. A Appl. Mat. 208,819–824 (2011). 6. J. Xu, A. W. Stevenson, D. Gao, M. Tykocinski, D. Lawrence, S. W. Wilkins, G. M. Clark, E. Saunders, and R. S. Cowan, “The role of radiographic phase-contrast imaging in the development of intracochlear electrode arrays,” Otology and Neur. 22,862-868 (2001). #175775 $15.00 USD Received 10 Sep 2012; revised 1 Oct 2012; accepted 4 Oct 2012; published 7 Jan 2013 (C) 2013 OSA 14 January 2013 / Vol. 21, No. 1 / OPTICS EXPRESS 710 7. P. Weiss, L. Obadia, D. Magne, X. Bourges, C. Rau, T. Weitkamp, I. Khairoun, J .M. Bouler, D. Chappard, O. Gauthier , and G. Daculsi, “Synchrotron x-ray microtomography (on a micron scale) provides three dimensional imaging representation of bone ingrowth in calcium phosphate biomaterials,” Biomat. 24, 4591 (2003). 8. P. Coan, F. Gruener, C. Glaser, T. Schneider, A. Bravin, M. Reiser, and D. Habs, “Phase contrast medical imaging with compact x-ray sources at the Munich-centre for advance photonics,” Nucl. Instr. and Meth. Phys. Res. A 608, S44–S46 (2009). 9. E. C. Ismaila, W. Kaabara, D. Garritya, O. Gundogdua, O. Bunkb, F. Pfeifferb, M. J. Farquharsond, and D. A. Bradleya, “X-ray phase contrast imaging of the bone-cartilage interface,” Appl. Rad. and Isot. 68, 767–771 (2010). 10. Q. Tao, D. li, L. Zhang, and S. Luo, “Using x-ray in-line phase-contrast imaging for the investigation of nude mouse hepatic tumors,” PLoS ONE 7, e39936, (2012). 11. A. A. Appel, J. C. Larson, S. Somo, Z. Zhong, P. P. Spicer, F. K. Kasper, A. B. Garson III, A. M. Zysk, A. G. Mikos, M. A. Anastasio, and E. M. Brey, “Imaging of poly( α-hydroxy-ester) scaffolds with x-ray phase-contrast microcomputed tomography,” Tissue Eng. Part C 18, (2012). 12. 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). 13. T.E. Gureyev, T.J. Davis, A. Pogany, S.C. Mayo, and S.W. Wilkins, “Optical phase retrieval by use of first Bornand Rytov-type approximations,” Appl. Opt. 43,2418–2430 (2004). 14. X. Wu, and H. Liu, “A general theoretical formalism for x-ray phase contrast imaging,” J. of X-ray Sci. and Techn.11, 33–42 (2003). 15. 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). 16. M. Langer, F. Peyrin, P. Cloetens, and J. -P. Guigay , “Quantitative comparison of direct phase retrieval algorithms in in-line phase tomography,” Med. Phys. 35, 4556 (2008). 17. 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). 18. 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). 19. X. Wu, and A. Yan, “Phase retrieval from one single phase contrast x-ray image,” Opt. Express 17,11187 (2009). 20. R. Hofmann, J. Moosmann, and T. Baumbach, “Criticality in single-distance phase retrieval,” Opt. Express 19, 25881–25890 (2011). 21. M. Langer, P. Cloetens, and F. Peyrin, “Fourier-wavelet regularization of phase retrieval in x-ray in-line phase tomography,” J. Opt. Soc. Am. A26, 1876-1881 (2009). 22. L.I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D. 60, 259–268 (1992). 23. A. Chambolle, “An algorithm for total variation minimization and applications,” J. Math. Imaging Vision 20, 89-97 (2004). 24. E. Y. Sidky, M. A. Anastasio, and X. Pan, “Image reconstruction exploiting object sparsity in boundary-enhanced x-ray phase-contrast tomography,” Opt. Express 18,10404 (2010). 25. 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). 26. A. W. M. van Eekeren, K. Schutte, and L. J. van Vliet, “Multiframe super-resolution reconstruction of small moving objects,” IEEE Trans. Im. Proc. 19, 2901-2912 (2010). 27. E. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inform. Theory 52,489–509 (2006). 28. X. Bresson, and T. F. Chan, “Fast dual minimization of the vectorial total variation norm and applications to color image processing,” Inverse Probl. and Imag. 2, 455–484 (2008). 29. A. Beck, and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems,” SIAM J. Im. Sci.2, 183-202 (2009). 30. A. Beck, and M. Teboulle, ”Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems,” IEEE Trans. Im. Proc. 18, 2419-2434 (2009). 31. T. Q. Pham, L. J. van Vliet, and K. Schutte, “Robust super-resolution by minimizing a gaussian-weighted l2 error norm,” J. Phys.: Conf. 124, 012037 (2008). 32. R. Barrett, R. Baker, P. Cloetens, Y. Dabin, C. Morawe, H. Suhonen, R. Tucoulou, A. Vivo, and L. Zhang. “Dynamically-figured mirror system for high-energy nanofocusing at the ESRF, Proc. SPIE 12, 813904 (2011). 33. R. Mokso, P. Cloetens, E. Maire,W. Ludwig, and J.-Y. Buffire, “Nanoscale zoom tomography with hard x-rays using Kirkpatrick-Baez optics, Appl. Phys. Lett. 90, 144104 (2007). #175775 $15.00 USD Received 10 Sep 2012; revised 1 Oct 2012; accepted 4 Oct 2012; published 7 Jan 2013 (C) 2013 OSA 14 January 2013 / Vol. 21, No. 1 / OPTICS EXPRESS 711


Introduction
Recently, the field of X-ray phase-contrast imaging (PCI) has been growing rapidly.X-ray PCI found applications in materials science, ranging from investigating the microstructure of carbon-based materials [1,2] to in-situ measurements of dynamic processes taking place in metal alloys and semiconductors [3][4][5].X-ray PCI is also entering the field of pre-clinical biomedical research, namely, small animal imaging and various ex-vivo/in-vitro studies [6][7][8][9][10][11].The increasing availability of the X-ray PCI techniques over the last years was stimulated by advances in instrumentation and phase retrieval algorithms.
The scope of our current investigation lays within phase retrieval for propagation-based Xray PCI.For more than a decade various algorithms were developed to permit accurate reconstruction of the specimen's phase and attenuation images from phase-contrast data acquired using the propagation-based approach [12].A major effort was aimed at the development of linear approximations to the image formation of PCI that would permit a stable solution of the resulting inverse problem [12][13][14][15][16]. Using these approximations, the phase and attenuation images of the specimen can be calculated from a series of phase-contrast images acquired at different propagation distances.To allow phase retrieval from a single phase-contrast image, methods based on prior information about the specimen's composition were developed [17][18][19].They are referred to as the so called phase-attenuation duality models.
In all linear phase retrieval models the accuracy of the reconstruction is a function of spatial frequency.Depending on the acquisition conditions, the signal-to-noise ratio (SNR), and the fitness of the linear approximation, the phase image of the specimen will be irrecoverable within a particular set of spatial frequencies [20,21].In the case of multi-distance phase retrieval [16] this can lead to large errors at low spatial frequencies while in single-distance approaches [17][18][19] artifacts are produced at middle and high spatial frequencies.In order to avoid large errors in the reconstructed images, most of the phase retrieval approaches rely on a so called L2-norm based regularization also known as Tikhonov regularization.When L2 regularization is used, the solution that has the minimum L2-norm (i.e.Euclidean norm) is promoted.This leads to a suppression of spatial frequencies that are ill-determined by the phase retrieval model or heavily corrupted by noise.Such a solution may not be optimal, especially when it results in a strong suppression of a large band of low frequencies in multi-distance retrieval methods.
Another regularization approach that is currently used in an increasing number of image reconstruction applications is called Total Variation (TV) minimization.It was initially developed for image denoising [22] and recently been introduced in such fields as deblurring, super-resolution, inpainting and tomography [23][24][25][26].The idea underlying TV minimization is to promote a solution that has the sparsest gradient.It was theoretically proven [27] that under certain conditions TV minimization allows exact reconstruction of signals with a sparse gradient from highly incomplete sets of observations.In cases where the gradient of the reconstructed image is not exactly sparse, TV minimization is nevertheless preferred to L2 regularization in many applications [28].
In this paper we introduce a TV minimization approach for solving the inverse problem of phase retrieval in propagation-based X-ray PCI based on various linear models.Any implementation of TV minimization can be chosen from a wide range of algorithms [25].Here we will present only results acquired with an algorithm based on the so called Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) [29].The original implementation of the algorithm was modified to include frequency weighting into the minimization scheme.Frequency weighting permits to account for the frequency-dependent nature of the signal-to-noise ratio and is shown to have a significant influence on the accuracy of the phase retrieval.

Materials and methods
In the following subsections we will describe how the phase retrieval problem of phase retrieval can be solved using iterative least-squares minimization and TV-minimization.In order to do so, we will introduce a matrix formalism which is uncommon in the field of phase retrieval.To keep our description compact we will refer the reader to the paper of M. Langer [16] for details regarding the theoretical background of X-ray phase retrieval algorithms.

Matrix formalism for phase propagation model
It is essential for the purpose of this paper that the convolution integral can be expressed using matrix formalism.Let us, for instance, define the propagated X-ray wavefront as a matrix product.Let A denote the set of square integrable functions R 2 → C with bounded support.Within the paraxial approximation, the X-ray field H D ∈ A propagated to a distance D from the object can be expressed as the convolution (⋆) of the unpropagated field T ∈ A with the Fresnel propagator P D ∈ A : ( If we denote the Fourier domain representations of T , P D , and H D by T , PD , and HD , respectively, we find that where f denotes the spatial frequency.We now discretize the Fourier domain, representing a spectrum by its values in a discrete set where PD is a diagonal matrix containing a discrete representation of the propagator PD ( f k ) and t is a vector that contains a discrete representation of T ( f k ).Consequently, a discretized representation hD of the propagated field HD can be defined in frequency domain: where the symbol (•) denotes the matrix-vector product.If we introduce a vector t corresponding to the discrete representation of the unpropagated field T , it is possible to construct a Toeplitz matrix P D such that the discretized propagated field h D can be expressed in the spatial domain: As matrix P D is dense, it is hard to compute Eq. ( 5) directly or solve the inverse problem for a large vector t in the spatial domain.However, the fact that matrix P D has a complimentary frequency space representation PD , both play an important role in computing the explicit deconvolution and iterative least-squares inversion.Since PD is a diagonal matrix, its eigenvalues correspond to the discrete representation of the propagator P D ( f k ).That property allows us to construct the well-known explicit deconvolution Eq. ( 8) in Section 2.2 and to analyse the properties of iterative least-squares inversion in Section 2.3.

Linear phase retrieval algorithms: L2-norm regularization
The standard approach to phase retrieval in in-line x-ray PCI relies on the fact that the observed phase-contrast image can be approximated by a linear transformation acting on some unknown image (or a combination of a phase and an attenuation image).This problem can be viewed as a system of linear equations: where A is a matrix that represents one of the linear phase models, b a vector containing the observed phase-contrast images, and x a vector containing the unknown images.Since the observed images are corrupted by noise and inversion of A is usually ill-posed, Eq. ( 6) is often replaced by an L2-norm regularized least-squares problem (a.k.a. a Tikhonov regularization problem): where ... 2 2 denotes the L2-norm and ε the regularization weight.The minimization of the first term guarantees the best fit of the linear phase-contrast model to the observed images, while the second term promotes solutions with the smallest L2-norm thereby suppressing noise and outliers.Given that matrix A has a diagonal representation in frequency domain Ã, Eq. ( 7) has an analytical solution according to which each j th frequency component x j can be calculated using the following expression: where Ã * j, j and | Ã j, j | 2 are respectively the conjugate and the squared magnitude of the diagonal matrix Ã .When in underdetermined or ill-conditioned system | Ã j, j | − → 0, the corresponding frequency component of x will be suppressed by L2 regularization:

Linear phase retrieval algorithms: TV minimization
Similar to the L2-regularized approach Eq. ( 6), Eq. ( 7) can be replaced by a TV-regularized version: where the second term x TV denotes the so called Total Variation norm of x which is defined by the gradient magnitude of the image x where ∇ h and ∇ v are the horizontal and vertical finite difference operators.It has to be noted that the TV norm defined by Eq. ( 11) has the dimension of 1/length (assuming x dimensionless), so ε TV has the dimension of length unlike the dimensionless ε L2 .Therefore the ratio between the optimal ε TV and ε L2 will in general depend on the choice of the pixel size.All values of ε TV given in this article implicitly have the dimension of the pixel size that is used in the calculation.
In contrast with the L2-regularized Eq. ( 7) the second term of Eq. ( 10) promotes solutions with sparse gradient magnitude.Expression (10) represents a non-smooth convex minimization problem and can be calculated numerically using one of the iterative TV minimization algorithms [25].During this study we used an implementation of TV minimization based on FISTA [29] which was introduced by Beck and Teboulle in [30].
The Beck and Teboulle algorithm can be viewed as an extension of dual gradient minimization [28].Minimizing Eq. ( 10) is achieved by splitting each iteration into two sub-problems (or steps), the so called gradient step and denoising step: 1. Gradient step: finding an image x 0 by reducing the unregularized L2 residual term A • x − b 2 2 in the beginning of each iteration, 2. Denoising step: followed by regularizing this intermediate image x 0 by minimizing: Alternating gradient and denoising steps has been shown previously to speed up the convergence without sacrificing accuracy [31].
The gradient step is based on FISTA where, in order to reach a high rate of convergence, the current guess is updated using information from two previous iterations: Here an intermediate vector y is introduced to take into account solution updates from two previous iterations, t n is a scalar that is determined at each iteration as , t 0 = 1, and L denotes a so called Lipschitz constant that can be calculated as the maximum eigenvalue of the product A * • A (see [29] for a detailed description).Operator D L,ε signifies the denoising step which can be implemented using Fast Gradient Projection (FGP) [23].During the denoising step, the TV norm of the current guess is minimized depending on the regularization weight ε TV and the Lipschitz constant L.
If matrix Ã is diagonal, all frequency components of the solution ỹn j can be updated during the gradient step (13) independently from each other: A frequency weighting vector ω j ≤ 1 is introduced to control the convergence of the algorithm.In general, convergence properties might vary among different minimization algorithms.However, the following observation is likely to be correct for algorithms similar to the one described above: frequency components of the solution x j that correspond to small matrix elements | Ã j, j | − → 0 will only be modified in the denoising step.Hence, its final value is determined solely by the TV-term of Eq. (10).The latter ensures a significant difference between how TV minimization and L2 regularization (see Section 2.4) are computing the frequency components of the solution that can not be retrieved from the observations.

Linear phase retrieval algorithms: models
As was explained in the previous section, TV minimization can be used for phase retrieval as long as the phase-contrast model can be expressed as a linear system Eq.( 6).
Let us introduce matrices C D and S D that represent the convolution with the real and imaginary parts of the Fresnel propagator P D respectively.Their frequency domain representations have the following form: where α j = πλ D| f j | 2 , λ stands for the X-ray wavelength and D is the effective propagation distance.Using matrices C D and S D we can construct matrix A and the corresponding update rule for the gradient step for the following linear phase retrieval models: CTF model.The CTF model is widely used for phase retrieval in cases when the specimen yields negligible attenuation and slowly varying phase [12].The Fourier transform of the phasecontrast image ĨD ( f ) is approximated by the Fourier transforms μ( f ) and φ( f ) of the projected attenuation and projected phase images of the specimen: Linear systems that are formed by combining Eqs. ( 16) for a set of m phase-contrast images {I D(1) , ..., I D(m) } can easily be represented in matrix form Eq. ( 6) (in frequency domain for | f j | > 0) by construction of Ã, x and b: Here the unknown vectors μ and φ are concatenated into a single vector x, vector b contains all observed images ĨD(i) and matrix Ã is obtained by concatenating pairs of matrices CD(i) and SD(i) , where each pair corresponds to a particular propagation distance D(i).Using this representation we can find an update rule similar to Eq. ( 14) for each j th frequency component of the unknown vector x.We will separate it into the update rules for μ and φ as follows: Calculation of Eq. ( 18) has to be carried out at each iteration of the gradient and denoising steps.Since the denoising step must be computed in the spatial domain, two inverse Fourier transforms (for μ and φ) must be calculated at each iteration before the denoising step and two Fourier transforms after the denoising step.
Mixed model.The Mixed model [14] is used for phase retrieval in cases with (significant) attenuation.In an approximated version of this model (assuming only the first two terms), the phase-contrast image ĨD ( f ) becomes: Here Ĩ0 ( f ) denotes the Fourier transform of the intensity image at zero distance.Ĩ0 ( f ) is fully determined by the attenuation image of the specimen and can be expressed in the spatial domain as I 0 = e −2µ .The linear system that describes a set of phase-contrast images { ĨD(1) , ..., ĨD(m) } based on Eq. ( 19) can be expressed through the following Ã, x and b: Here we treat the element-wise product of the intensity image and the phase image (I 0 ϕ) as an unknown image independent from the unknown intensity image I 0 .Using this representation we can write down the update rule Eq. ( 14) for the gradient step as: Ĩn−1 0, j + 2 S( j, j) Phase-attenuation duality models.These models can be used when the specimen has a homogeneous composition or, in the limited range of X-ray energies, when the specimen is composed of light elements [17,19].In duality models the phase and attenuation images of a specimen are assumed to be proportional to each other, i.e. σ = ϕ µ , permitting to reduce the number of unknown variables.In this study we consider two duality models: which was derived from the CTF model [18] and : which was derived from the Mixed model [19].Both models can be used to retrieve the projected attenuation image of the specimen based on a single phase contrast image ĨD ( f ) acquired at a suitable distance D > 0 as well as using a set of m phase-contrast images { ĨD(1) , ..., ĨD(m) } recorded at different distances.The update rule for the gradient step based on Eq. ( 22) is then simplified into: where Bi = (σ S( j, j) D(i) − C( j, j) D(i) ).The update rule for the gradient step based on Eq. ( 23) has the following form: where Bi = ( C( j, j) D(i) + (α − σ ) S( j, j) D(i) ).

Simulations
As indicated before, TV minimization permit an accurate solution for a class of inverse problems based on severely incomplete sets of observations.The underlying assumption of all TV minimization approaches is that the unknown signal must have a sparse gradient magnitude.In the field of image reconstruction such assumption is fulfilled when the reconstructed image is piecewise constant, i.e. the intensity only changes in a step-like manner.

Phantom image with sparse gradient magnitude
Our first demonstration of phase retrieval based on TV minimization uses a 'flat' piece-wise constant phantom.Here TV minimization is expected to yield high accuracy.However, in the following subsections we will abandon this restriction and use a 'spheres' phantom to demonstrate the performance of TV minimization in more realistic cases.
Figure 1 shows the comparison between L2-regularized and TV-regularized phase retrieval for the CTF model.The ground truth projected attenuation and phase images (256 × 256 pixels) were computed for a randomly generated composition of overlapping polyethylene disks immersed in a layer of water with a total thickness d = 0.1 mm (Fig. 1(a)).Subsequently phasecontrast images I D 1 and I D 2 are generated using Fresnel propagation for a monochromatic Xray energy of 20KeV and a pixel size of 1 µm (Fig. 1(b, c)).In the current simulation we used propagation distances {0 m, 1 m}.Gaussian noise with standard deviation of 0.02 was added to both images mimicking acquisition with poor SNR.The optimal regularization parameters ε L2 and ε TV were chosen such that the overall Root Mean Square Error (RMSE) is minimal for the listed conditions (see Section 3.3).In practice, they must be derived from an estimate of the SNR of the measured images.No stopping criteria was used in the gradient step of the TV minimization, instead all reconstructions were computed using 1000 iterations in order to guarantee convergence.The resulting solutions are shown in Fig. 1(d, e).The frequency spectrum (we will use this term for the angular average of the magnitude of 2D Fourier transform of the image) of the error magnitude of the solutions can be seen on Fig. 2 (right).The frequency spectrum of the error was normalized against the frequency spectrum of the ground truth image.It is evident from this graph that when the optimal regularization weights are used, TV regularization is significantly more accurate then L2 regularization at all frequencies where the direct phase retrieval problem is undetermined (for frequencies components with | Ã j, j | − → 0).Note that TV minimization can yield higher error within particular bands of frequencies to promote a sparse gradient.We have already pointed out in Section 2.2 and Section 2.3 that there are major differences in how TV minimization treats frequency components of the solution that correspond to | Ã j, j | − → 0. In order to confirm this, we compared the results of the TV method with ε TV = 0.02 and ε TV = 0 (no TV minimization) when applied to noise-free data.Figure 2(left) shows that the frequency spectrum of the normalized difference between the regularized and non-regularized solutions.It is evident that the regularization strength depends heavily on | Ã j, j |.
So far we used scalar regularization weights ε L2 and ε TV in both regularization approaches.One can achieve a significantly improved accuracy by using a frequency dependent regularization weighting instead.In L2 regularization, the scalar ε L2 can be replaced by the frequency dependent factor ε L2, j ∼ SNR −1 ( f j ) (i.e.Wiener deconvolution).In our investigation we as- sumed the SNR of the reconstructed image to be proportional to its frequency spectrum.The latter can be estimated from a preliminary reconstruction based on a scalar regularization weight using either L2 or TV regularization.Such estimation is demonstrated in Fig. 3(left).The estimate of the image SNR was introduced into TV minimization approach using the frequency weighting vector ω j ∼ SNR( f j ) in Eq. ( 14).Resulting L2-and TV-regularized solutions are depicted in Fig. 1(f) and Fig. 1(g) respectively.It can be seen in Fig. 3(right) that the error of TV-regularized solution is significantly lower than the error of L2-regularized solution in a wide frequency range.

Realistic phantom
In order to demonstrate the performance of the TV-minimization in realistic cases we tested L2-regularization and TV-minimization phase retrieval approaches for phantoms with a nonsparse gradient.Figure 4 shows reconstructions of the projected phase image of a composition of randomly sized and positioned polyethylene spheres immersed in water (the rest of simulation parameters match those from Section 3.1).Reconstructed images are obtained using both L2-and TV-regularization based on the CTF model, the Mixed model as well as their phaseattenuation duality modifications.In phase retrieval based on duality models we have used a single simulated phase-contrast image with propagation distance 1 m.The corresponding normalized frequency spectra of the error-magnitude are depicted in Fig. 5.It is evident that the TV-regularized solutions yield lower error in comparison with the L2-regularized ones in a broad range of frequencies.It is also apparent that given the parameters used in the current simulation (low attenuation and moderate phase changes), reconstructions obtained with the CTF model and the Mixed model are virtually indistinguishable, while there is some discrepancy with the models that are based on phase-attenuation duality assumption.

Optimal regularization weights
The accuracy of both regularization methods considered in this paper strongly depend on the choice of the regularization weights ε L2 and ε TV .In practice, these parameters are estimated from the measured data or chosen in some heuristic manner.In order to investigate how the choice of the regularization weight affects the accuracy of phase retrieval, we measured the total RMSE of the reconstructed phase images varying two parameters: standard deviation (STD) of the Gaussian noise and the regularization weight of the phase retrieval algorithm.The phasecontrast images were simulated using the polyethylene spheres phantom using the simulation parameters described in Section 3.1 and 3. L2-regularized solutions with noise STD = 0.01, 0.02, 0.05 and 0.1 respectively.Curves TV-1, TV-2, TV-3 and TV-4 represent the RMSE of the TV-regularized solutions with corresponding noise STD.It can be seen that TV-regularization yields solutions with a 4-7 times lower total RMSE in comparison to those obtained with L2 regularization.It is also evident that estimation of the optimal regularization parameter is more important for TV-regularization then for L2 regularization as it affects the RMSE to a larger extent.
Along with various types of additive noise, non-linear effects that are not taken into account by the phase retrieval models can become an important source of errors.The fraction of signal to error due to non-linearity of the CTF model is, in a certain range of conditions, proportional to | sin(α)| [20].That property allows to treat the non-linearity error as another form of noise.Figure 6(right) demonstrates the values of RMSE of the phase retrieval based on the CTF model applied to noise-free phase-contrast images simulated for phantoms of different thicknesses.Curves L2-1, L2-2, L2-3 and L2-4 show the RMSE of the L2-regularized phase retrieval for phantoms with a total thickness = 0.1 mm, 0.25 mm, 0.5 mm and 1 mm respectively.Curves TV-1, TV-2, TV-3 and TV-4 represent the RMSE of the corresponding TV-regularized solutions.The thickness of the phantom has a dramatic effect on the accuracy of the phase retrieval based on the CTF model since it introduces larger variations in projected attenuation and phase images which leads to greater non-linearity of the observed phase-contrast image.It is evident that TVregularized phase retrieval yields similar accuracy with the L2-regularized solution in the case of thin phantom.However the advantage of TV regularization becomes significant for thicker phantoms.
TV-1 Fig. 6.Influence of the regularization weight on samples with increasing noise levels and on samples of increasing thickness.RMSE of the solution is shown against the magnitude of regularization parameter for L2 regularization (crosses) against TV regularization (plus signs).Left: each curve shows errors for different noise levels, STD = 0.01, 0.02, 0.05, 0.1.Right: each curve shows errors for different specimen thickness, 0.1 mm, 0.2 mm, 0.5 mm, 1 mm.
The focus was used as a point source for projection microscopy [33] giving effective pixel size of 53 nm.Images were acquired at four propagation distances {27.4 mm, 28.3 mm, 31.8 mm, 40.3 mm}.
Figure 7 shows reconstructed phase images of the lines and dots pattern and so called Siemens star pattern.The L2-regularized solutions are shown in the top row, while the bottom row shows the TV-regularized solutions.Both phase-retrieval approaches were based on phase-attenuation duality CTF model.The reconstructions depicted in sub-figures (a, b, e and f) are based on phase-contrast images acquired at four different distances.The reconstructions from sub-figures (c, d, g and h) are based on a single phase-contrast image acquired at a propagation distance 27.4 mm.
It is evident that TV regularization permits a very high quality reconstruction of lines and dots pattern based only on a single phase-contrast image.The high frequency ripple-like artifacts are efficiently suppressed while the pattern is accurately reconstructed.That can be explained by the fact that the reconstructed image has a sparse gradient magnitude and that spatially localized structures such as dots and lines have a broad footprint in frequency space.The problem of phase retrieval based on a single phase-contrast image is ill-posed within a set of spatial frequencies that depends on the acquisition parameters.TV regularization allows to fill-in these gaps in frequency space by applying the constraint of sparse gradient magnitude.Phase retrieval of the Siemens star pattern shows that the structures which are periodic in space are harder to reconstruct.Their footprint in frequency space is localized and might be completely irrecoverable from a single phase-contrast image with the given acquisition parameters.

Conclusion
Phase retrieval in propagation based X-ray PCI can be improved using iterative TV minimization algorithms.Reconstructions based on simulated and experimental data show that phase retrieval based on TV minimization can significantly outperform the current practice, a deconvolution approach with L2 regularization.TV minimization can be used with different linear phase retrieval models including the CTF model, the Mixed model and the phase-attenuation  duality models.Although the method works best for specimen that adhere to the constraints, the method also improves the phase reconstruction for specimen that are not exactly sparse in their gradient magnitude.TV minimization provides an effective regularization instrument for solving an underdetermined linear systems of equations.Analysis of the Fourier spectrum of the error of the reconstructed phase images clearly demonstrates that TV minimization allows partial recovery of the solution within the frequency bands that are undefined by the particular phase-contrast model or corrupted by noise.That feature of TV minimization allows effective suppression of the high frequency artifacts in single-distance phase retrieval based on phaseattenuation duality models and permits more accurate reconstruction of low spatial frequencies in multi-distance approaches.
TV minimization finds the solution in the form of a TV-regularized least-squares fit and it does not require the knowledge of the attenuation part of the specimen.A frequency dependent estimation of the signal-to-noise ratio can be used for each observed image, dramatically improving the accuracy of reconstruction.Simulations show that TV regularization can suppress errors that occur both due to additive noise and due to non-linearity of the phase propagation.Experimental data has shown that TV minimization can also significantly improve the accuracy of phase reconstruction of real specimen that comply with gradient sparsity condition and imaged under realistic circumstances.These results could allow to decrease the number of phase-contrast images that are needed for the accurate image reconstruction in some applications of X-ray PCI.This is particularly important for in situ experiments or to reduce radiation damage to the specimen.

Fig. 1 .
Fig. 1.Phase reconstructions based on the simulated data of the 'flat' phantom for TV and L2 regularization with and without frequency weighting (see Sec. 3.1 for simulation parameters).(a) Ground truth.(b) Intensity image at zero distance with Gaussian noise.(c) Propagated phase-contrast image at 1m with Gaussian noise.(d) L2-regularized solution, no frequency weighting.(e) TV-regularized solution, no frequency weighting.(f) L2regularized solution, with frequency weighting.(g) TV-regularized solution with frequency weighting.

#Fig. 2 .
Fig. 2. The radial frequency spectrum (angular averaged) of the reconstruction.Dotted line shows | Ã j, j | on both graphs.Left: normalized difference between the solutions with and without TV-regularization, i.e. (ε TV = 0.02) vs. (ε TV = 0).Right: normalized errormagnitude of the TV-regularized solution (solid line) and the L2-regularized solution (dashed line).The graphs are normalized against the radial frequency spectrum of the ground truth image.

Fig. 3 .
Fig.3.Radial frequency spectrum (angular averaged) of the retrieval results using the frequency dependent regularization weights (i.e.frequency weighting).Dotted line shows | Ã j, j | on both graphs.Left: estimated frequency spectrum of the specimen (solid line).Right: normalized error-magnitude of the TV-regularized solution (solid line) and L2regularized solution (dashed line) with frequency weighting.

Fig. 4 .
Fig. 4. Reconstructions based on the simulated data of the 'spheres' phantom for different propagation models.All results are obtained using frequency weighting.Top row of images (a, c, e, g) show L2-regularized solutions.The bottom row of images (b, d, f, h) show TVregularized solutions.Solutions (a, b) are based on the CTF model, (c, d) on the Mixed model, (e, f) on the dual-CTF model and (g, h) on the dual-Mixed model.

Fig. 5 .
Fig. 5.The radial frequency spectrum (angular averaged) of the reconstructions for different propagation models.Dotted line shows | Ã j, j | on all graphs.Normalized error magnitude of the TV-regularized solution (solid line) and L2-regularized solution (dashed line).

Fig. 7 .
Fig. 7. Phase retrieval from experimental data using the dual-CTF model.(a-d) Phase retrieval of a 'dots and lines' pattern.(e-h) Phase retrieval of a 'star' pattern.(a, e) L2regularized solution based on 4 images recorded at different propagation distances.(b, f) TV-regularized solution based on 4 images recorded at different propagation distances.(c, g) L2-regularized solution based on a single recorded image.(d, h) TV-regularized solution based on a single recorded image.