Several methods of phase retrieval for in-line phase tomography have already been investigated based on the linearization of the relation between the phase shift induced by the object and the diffracted intensity. In this work, we present a non-linear iterative approach using the Frechet derivative of the intensity recorded at a few number of propagation distances. A Landweber type iterative method with an analytic calculation of the Frechet derivative adjoint is proposed. The inverse problem is regularized with the smoothing L 2 norm of the phase gradient and evaluated for several different implementations. The evaluation of the method was performed using a simple phase map, both with and without noise. Our approach outperforms the linear methods on simulated noisy data up to high noise levels and thanks to the proposed analytical calculation is suited to the processing of large experimental image data sets.
© 2011 OSA
X- ray microtomography is widely used in a large variety of applications such as bone imaging [1–7] and material science [8, 9]. Yet it is difficult to image soft tissue, especially with dense structures. Several techniques to obtain phase contrast have been developed [10–12]. Coupling x-ray microtomography and phase contrast allows to improve the sensitivity of the method [13–15].
If the spatial coherence of the x-rays is sufficient, phase contrast can be achieved by letting the beam propagate in free-space after interaction with the object. In this case, the contrast is related to the 2D Laplacian of the phase shift. The relationship between the phase shift induced by a sample and the intensity recorded at a sample-to-detector distance D relies on Fresnel diffraction theory. Such phase contrast radiographs can be used to define an inverse problem to retrieve the phase shift induced by the object. The phase shift is proportional to a projection through the refractive index distribution in the object, and therefore phase retrieval can be coupled to tomography. The phase shift induced by the object is first retrieved for each projection angle. The reconstruction algorithm then yields a 3D reconstruction of the refractive index.
In the hard x-ray region, phase contrast can give a dramatic increase in the purely refractive part of the sensitivity. In several cases, such as biomaterials, bone and small animal imaging, both dense and soft tissues are present. The x-ray energy has to be chosen so that there will be both strong absorption and phase contrast in the recorded images.
Several phase retrieval methods have already been investigated. Some are based on the Transport of Intensity Equation (TIE) [16, 17] and the use of a series of image measurements obtained at different propagation distances. Theses methods can be refined by other techniques like Gerchberg-Saxton-Fienup algorithms [18, 19]. The others rely on the contrast transfer function (CTF) [10, 20–22] or on a mixed approach between the two former methods . All these approaches are based on a linearized relation between the phase and the intensity, thus involving some approximations to the direct problem.
In this work, a new iterative method is proposed to invert the non-linear function I(φ) relating the measurements to the phase map. Since the inverse operator I −1(φ) is not bounded, the phase retrieval problem considered in this work is ill-posed. Thus Tikhonov regularization is included. The problem is expressed as the minimization of a non-linear regularized functional with the Frechet derivative of I(φ). An analytic expression of the adjoint of the Frechet derivative allows to speed up the algorithm and reduce memory requirements.
Section 2 presents our formal development of the direct problem as well as the linear approaches used in the previous works . In section 3, the problem is expressed as the minimization of a non-linear regularized functional solved by a Landweber-type iterative method. The analytic expression of the adjoint of the Frechet derivative is derived. Section 4 details the simulation procedure used to invert noise-free and noisy data. Results comparing different variants of the algorithm are reported and discussed in section 5.
2. The direct problem and the linearized inverse problem phase contrast imaging
2.1. The direct problem
Considering an object illuminated with partially coherent x-rays of wavelength λ, its interaction with beam can be described by its complex refractive index, usually written as :
The absorption, a(x), and phase shift φ(x) induced by the object can be considered as projections of the absorption and refraction index respectively:
The intensity detected at a distance D is given by the squared modulus of the exit wave:24]. The Fourier transform ℱ of a function g(x) is defined as:
2.2. The linear inverse problem
Most existing phase retrieval approaches are based on a linearization of Eq. (4) valid under some assumptions. In particular, assuming that the absorption is slowly varying it has been shown that the Fourier transforms of the intensity can be written :24]. It will thus be used as a comparison with the non-linear method. To include several distances, a least squares minimization procedure is used.
To yield stable solutions, regularization techniques are required. Several methods of regularization have been tested [24–26] to solve this linear inverse problem, such as classical quadratic Tikhonov regularization and wavelet shrinkage. By introducing the phase-absorption product ψ(x) = I 0(x)φ(x), the phase retrieval problem can be solved iteratively as:
3. The non-linear inverse phase retrieval problem
The problem considered in this work is an ill-posed problem in the sense that the operator is not a bounded linear operator. Prior knowledge about the phase is necessary to retrieve the phase from the intensity measurements. In order to solve the inverse problem, we have searched for a smooth solution because failure of convergence and stagnation of the iterates away from the solution were observed when no prior information was included. Numerous iterative algorithms have been proposed and analysed in the literature, based on regularization functionals defined in a Hilbert space context . In this work, in order to reconstruct the phase in a non-linear framework, we propose an iterative Landweber type method.
3.1. A Landweber type iterative method
In view of Eq. (4), the intensity can be regarded as a function of φ, ID(φ). The operator ID(φ) can be considered as a non-linear operator which is Frechet differentiable in its domain. In the following, we will consider the phase having a Lipschitz bounded support Ω and that the domain 𝒟 [ID(φ)] of the operator ID(φ) belongs to the functional Sobolev space :
The optimality condition is
The step length parameter τk is chosen in order to minimize the Tikhonov’s functional along the descent direction:
This optimality condition defines the descent direction of our steepest descent iterative method, the next phase iterate φk +1 is obtained from the iterate φk with:
Starting from the current phase estimate φk at the iteration k, a linear search procedure is introduced with a variable step τk yielding the following modification of the standard Landweber method
The regularizing parameter is chosen by trial-and-error in order to obtain the best decrease of J(φ). These equations correspond to a Landweber method with a regularizing smoothing term and an adaptative step length. This type of method has been used in several works on non-linear ill-posed problems [28, 29]. In this work, a convergence analysis of the method is not performed. We show that it gives good results for typical phase maps found in in-line phase tomography applications, with various levels of noise.
3.2. Update of iterates
The calculation of the descent direction and of the iterates can be performed with finite differences with an adaptive step length, the implicit filtering method or with an analytic calculation of the adjoint of the diffracted intensity.
3.2.1. Analytical expression of the adjoint of the diffracted intensity derivative
In the case of the phase-intensity relationship, the Frechet derivative of the operator ID(φ) at the point φk is the linear operator Gk defined by the relation:
The linear operator ID′(φk)(ɛ) = Gk(ɛ) can be calculated explicitly as:
Using integrations by parts, it can be shown that the adjoint operator is defined by:
3.2.2. Numerical calculation by implicit filtering
In this work, we retrieve the phase from noise-free and noisy data. It is well known that it is difficult to calculate the gradient of a noisy function. We have thus applied the implicit filtering method described by Kelley et al in [30–32]. In its simplest unconstrained form, implicit filtering is the steepest descent algorithm with finite difference gradients, where the difference increment varies as the iteration progresses. Because the gradient is only an approximation, the computed descent direction may fail to be a descent direction, and the line search may fail. In this case the difference increment is reduced. It has also been shown that the performance of implicit filtering with a central difference gradient is superior to that with forward difference gradient [30–32]. The derivatives in the gradient ID′ (φ) are thus approximated by centered differences formulas. It is noteworthy that this finite difference method requires 2N evaluations of the Tikhonov functional, where N is the phase vector dimension and it leads to a squared matrix of size N 2. In this work, the phase increment is set to 0.05 rad which is the estimated noise level on the phase. This phase increment must be small enough so that the linearization is valid and higher than the noise level.
4.1. Simulation of the diffracted intensity
Following , the imaging system was simulated in a deterministic fashion. Two phantoms were defined, one for the absorption coefficient and one for the refractive index decrement. Theoretical values for the absorption coefficient and the refractive index, δ and β, for different materials at 24 keV were used (λ = 0.5166Å) for the phantom. Propagation in free-space was simulated using Eq. (4). The original phase map to be retrieved is displayed in Fig. 1(a), together with the absorption map in Fig. 1(b) and the corresponding Fresnel diffraction pattern for D=1.4 m without noise in Fig. 1(c). The convolution product was calculated by Fourier transforms and the intensity has been obtained as the squared modulus in the spatial domain of this convolution. Using the free-space propagation equation, images were calculated for the eight propagation distances 0.2 m, 0.4 m, 0.6 m, 0.8 m, 1.0 m, 1.2 m, 1.4 m and 1.6 m. These distances are suitable for testing and comparing the phase retrieval methods.
The intensity and phase values were discretized on a regular grid and they are considered as vectors of a 75 × 75 dimensional real space ℝ75×75. The Frechet derivative Gk calculated with the finite difference method at the point φk is a matrix of ℝ5475×5475. The phase contrast images were all corrupted with additive Gaussian white noise with various peak-to-peak signal to noise ratios (PPSNR), between 24 dB and 0 dB. The peak-to-peak signal to noise ratio is defined by:
4.2. Initialization and stopping conditions
Our phase retrieval algorithms are not globally convergent algorithms. The method will be quantitatively evaluated by measuring the normalized mean square error (NMSE) using the L 2(Ω) norm. If φ * is the phase to be recovered and φk the current estimate, the NMSE is calculated as:
The mixed approximation of the linear problem has been used as the starting point of our simulations. It is displayed in Fig. 1(d). The initial NMSE is 0.147. This a priori guess of the solution ensures the convergence of the algorithms.
To avoid problems at the image boundary, the phase support is assumed to be included in the 70 × 70 inner pixels, and the border pixels have been fixed to zero. Since several intensity maps obtained at different distances are available, the inverse problem may be split into a finite number of sub-problems. In order to take into account more than one intensity map, we propose two variants of a cyclic iteration over the distances. When one steepest descent iteration is performed for each image recording distance, this method will be called a Kaczmarz type method in the following discussion. When a hundred or more iterations are performed for each distance, this methodology will be called a sequential one.
In the following, the iterate φk +1 is accepted if the following two conditions are satisfied:
5. Results and discussion
In the following, we have compared the numerical results for the simple test case introduced above with the following algorithms:
- A 1: Kaczmarz type finite difference Landweber method with a regularization term or without any regularization (α = 0).
- A 2: Sequential type finite difference Landweber method with a regularization term or without regularization term performed with N=100 iterations for each distance.
- A 3: Kaczmarz type analytic Landweber method with the former stopping conditions with a regularization term or without any regularization (α = 0).
These algorithms differ by the number of descent iterations performed for each propagation distances, by the calculation of the adjoint of the Frechet derivative which is based on finite difference or on an analytic expression and by the regularization term. The difference between the algorithms A 1 and A 2 is the way of the cyclic iteration in the intensity maps is performed. The last algorithm A 3 represents the Kaczmarz type analytic Landweber method with the stopping conditions (28), (29) and (30) both without regularization term or with the regularization method.
Figure 2 displays the difference maps between the solution and the true phase to be retrieved. Figure 2(a) is the initial error map obtained from the mixed approach showing that errors are mostly concentrated at the edges. Figures 2(b), 2(c) and 2(d) illustrate the differences maps respectively obtained with the algorithm A 1 (PPSNR=24 dB), A 2 (noisy-free) and A 3 (PPSNR=24 dB). The regularization parameter was set to α = 0.01. In order to have more quantitative information about the convergence rates and to compare the algorithms, we have studied NMSE for the phase shift as a function of the number of iterations. Figure 3 shows the evolution of the NMSE as a function of the number of iterations for the different algorithms on the noise-free and noisy data (PPSNR=24 dB) for α = 0.01.
The errors on the phase have been significantly reduced with all algorithms. The use of several distances improves the reconstruction because it allows a better coverage of the frequency domain and it improves the statistics. Yet, a very slow convergence is obtained with the algorithm A 2. Kaczmarz type methods are thus to be preferred.
It is well known that the Landweber iteration has a regularizing effect, the number of iterations being the regularization parameter. The regularization term is not crucial in that case and improves only slightly the convergence rate. Yet, a divergence of the iterates away from the solution has been observed for other phase maps, if this term is not included in the functional. The regularization parameter has been selected from trial-and-error. As displayed in Figs. 3(a) and 3(b) for noise-free and noisy data, the use of a weak (α = 0.01) smoothing regularization yields good phase retrieval convergence results.
In Fig. 3 it can observed that the algorithm A 3 has good convergence properties. The algorithm A 3 is also much faster since the large scale matrix used in the finite difference methods is replaced by the analytic expression of the adjoint of the Frechet derivative. It should be noted that at the end of the iterations, the condition of Eq. (30) is fulfilled and that we have to stop the iterations considering the noise level.
The phase maps obtained with the algorithm A 3 for the noise-free data and noisy data (PP-SNR=24 dB) with a smooth regularization α = 0.01 are displayed in Figs. 4(a) and 4(b) respectively. The final NMSE are 0.095 for noisy-data (PPSNR=24 dB) and 0.09 for noise-free data.
The phase retrieval error may still be decreased. The drawbacks of the regularization functional ||∇φ||L2 are well-known. An isotropic smoothing effect is obtained and the boundaries are not well preserved. This is obvious in Fig. 2(c). The noise is suppressed but the high values of the gradient are too greatly penalized on the edge. In future works, the gradient ∇φ may be replaced by a non-linear functional of φ as in semi-quadratic regularization or by the bounded variation semi-norm  or by anisotropic terms.
In this work, several non-linear iterative approachs for phase retrieval have been proposed. The methods investigated previously were based on the linearization of the relation between the phase shift induced by the object and the diffracted intensity. They have used the Transport Intensity Equation (TIE), the Contrast Transfer Function (CTF), or mixed approaches. Our nonlinear iterative approaches uses the Frechet derivative of the intensity recorded at small number of propagation distances. We have compared the convergence rates for three algorithms, two based on the finite difference gradient and one, on the analytic expression of the gradient.
The best results are obtained when the inverse problem is regularized with the smoothing L 2 norm of the phase gradient. The best convergence rates are found when the various distances are treated with a Kaczmarz type method where one descent iteration is performed for each distance. The evaluation of the method was performed using a simple phase map, both with and without noise. For the simulated data, the normalized mean square error was measured. Tikhonov regularization based on linear filtering used in this work has some well-known drawbacks since it does not only smooth noise but also blurs important features such as edges. To avoid these shortcomings, non-linear partial differential diffusion equations may be useful. Therefore, other regularization methods will be tested in future work. Our approach outperforms the linear methods on simulated noisy data for PPSNR above 20dB and the nonlinearities of the Fresnel diffraction are well taken into account. The analytic calculation of the adjoint of the Frechet derivative speeds up the calculations and overcomes memory limitations due to the Frechet derivative matrix. Thus this method opens promising perspectives to process experimental data in various applications.
References and links
2. M. Salomé, F. Peyrin, P. Cloetens, C. Odet, A. M. Laval-Jeantet, J. Baruchel, and P. Spanne, “A synchrotron radiation microtomography system for the analysis of trabecular bone samples,” Med. Phys. 26, 2194–2204 (1999). [CrossRef]
3. S. Nuzzo, F. Peyrin, P. Cloetens, J. Baruchel, and G. Boivin, “Quantification of the degree of mineralization of bone in three dimensions using synchrotron radiation microtomography,” Med. Phys. 29, 2672–2681 (2002). [CrossRef]
4. S. Bayat, L. Apostol, E. Boller, T. Brochard, and F. Peyrin, “Quantification of the degree of mineralization of bone in three dimensions using synchrotron radiation microtomography,” Nucl. Instr. Meth. Phys. Res. A 548, 247–252 (2005). [CrossRef]
5. C. Chappard, A. Basillais, L. Benhamou, A. Bonassie, N. Bonnet, B. Brunet-Imbault, and F. Peyrin, “Comparison of synchrotron radiation and conventional X-ray microcomputed tomography for assessing trabecular bone microachitecture of human femoral heads,” Med. Phys. 33, 287–293 (2003).
6. M. Ito, S. Ejiri, H. Jinnai, J. Kono, S. Ikeda, A. Nishida, K. Uesugi, N. Yagi, M. Tanaka, and K. Hayashi, “Bone structure and mineralization demonstrated using synchrotron radiation computed tomography (SR-CT) in animal models: preliminary findings,” J. Bone Miner. Metab. 21, 3568–3577 (2006).
7. G. J. Kazakia, A. J. Burghardt, S. Cheung, and S. Majumdar, “Assessment of bone tissue mineralization by conventional X-ray microcomputed tomography,” Med. Phys. 33, 3170–3179 (2008). [CrossRef]
8. J. Baruchel, E. Marire, P. Merle, and G. Peix, X-ray Tomography in Material Science (Hermes Science Publications, 2000).
9. U. Bonse, “Developments in X-ray tomography II,” Proc. SPIE 3775 (1999).
10. P. Cloetens, R. Barrett, J. Baruchel, J. P. Guigay, and M. Schlenker, “Phase objects in synchrotron radiation hard X-ray imaging,” J. Phys. D.:Appl. Phys. 29, 133–146 (1996). [CrossRef]
11. S. W. Wilkins, T. E. Gureyev, D. Gao, A. Pogany, and A. W. Stevenson, “Phase contrast imaging using polychromatic hard X-rays,” Nature (London) 384, 335–338 (1996). [CrossRef]
12. S. W. Wilkins, T. E. Gureyev, D. Gao, A. Pogany, and A. W. Stevenson, “Quantitative phase imaging using hard X-rays,” Phys. Rev. Lett. 77, 2961–2964 (1996). [CrossRef]
13. P. Cloetens, M. Pateyron-Salome, J. Y. Buffiere, G. Peix, J. Baruchel, F. Peyrin, and M. Schlenker, “Observation of microstructure and damage in materials by phase sensitive radiography and tomography,” J. Appl. Phys. 81(9), 5878–5886 (1997). [CrossRef]
14. A. Momose, T. Takeda, Y. Itai, A. Yoneyama, and K. Hirano, “Phase-contrast tomographic imaging using an X-ray interferometer,” J. Synchrotron. Rad. 5, 309–314 (1998). [CrossRef]
15. T. Weitkamp, C. David, O. Bunk, J. Bruder, P. Cloetens, and F. Pfeiffer, “X-ray phase radiography and tomography of soft tissue using grating interferometry,” Eur. J. Radiol. 68, S13–S17 (2008). [CrossRef]
16. D. M. Paganin, Coherent X-Ray Optics (Oxford University Press, 2006). [CrossRef]
17. T. E. Gureyev and K. A. Nugent, “Phase retrieval with the transport-of-intensity equation. II. orthogonal series solution for nonuniform illumination,” J. Opt. Soc. Am. A 13, 1670–1682 (1996). [CrossRef]
18. T.E. Gureyev, “Composite techniques for phase retrieval in the Fresnel region,” Opt. Commun. 220, 49–58 (2003). [CrossRef]
20. J. P. Guigay, “Fourier transform analyis of Fresnel diffraction patterns in in-line holograms,” Optik 46, 12–125 (1977).
21. S. Zabler, P. Cloetens, J.-P. Guigay, J. Baruchel, and M. Schlenker, “Optimization of phase contrast imaging using hard X-rays,” Rev. Sci. Instrum 76, 1–7 (2005). [CrossRef]
22. 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]
23. M. Born and E. Wolf, Principles of Optics (Cambridge University Press, 1997).
24. M. Langer, P. Cloetens, J. P. Guigay, and F. Peyrin, “Quantitative comparison of direct phase retrieval algorithms in in-line phase tomography,” Med. Phys. 35, 4556–4565 (2008). [CrossRef]
25. M. Langer, P. Cloetens, and F. Peyrin, “Fourier-wavelet regularization of phase retrieval in X-ray in-line phase tomography,” J. Opt. Soc. Am. A 28, 1877–1882 (2009). [CrossRef]
26. M. Langer, P. Cloetens, and F. Peyrin, “Regularization of phase retrieval with phase attenuation duality prior for 3D holotomography,” IEEE Trans. Image Process. 19, 2425–2436 (2010). [CrossRef]
27. O. Scherzer, M. Grasmair, H. Grossauer, M. Haltmeier, and F. Lenzen, Variationnal Methods in Imaging (Springer Verlag, 2008).
28. M. Hanke, A. Neubauer, and O. Scherzer, “A convergence analysis of the landweber iteration for nonlinear ill-posed problems,” Numer. Math. 72, 21–37 (1995). [CrossRef]
29. I. Daubechies, M. Fornasier, and I. Loris, “Accelerated projected gradient method for linear inverse problems with sparsity constraints,” J. Fourier Anal. Appl. 14, 764–792 (2008). [CrossRef]
30. C. T. Kelley and P. Gilmore, “An implicit filtering algorithm for optimization of functions with many local minima,” SIAM J. Optm. 5, 269–285 (1985).
31. D. Stoneking, G. L. Bilbro, R. Trew, P. Gilmore, and C. T. Kelley, “Yield optimization using a gaas process simulator coupled to a physical device model,” IEEE Trans. Microwave Theory Tech. 40, 1353–1363 (1992). [CrossRef]
32. C. T. Kelley, “Iterative methods for optimization,” Frontiers in Applied Mathematics (SIAM, 1999).