Phase retrieval in X-ray phase-contrast imaging suitable for tomography

doi:10.1364/OE.19.010359

« Show journal navigation

Phase retrieval in X-ray phase-contrast imaging suitable for tomography

Anna Burvall, Ulf Lundström, Per A. C. Takman, Daniel H. Larsson, and Hans M. Hertz »View Author Affiliations

Optics Express, Vol. 19, Issue 11, pp. 10359-10376 (2011)
http://dx.doi.org/10.1364/OE.19.010359

View Full Text Article

Acrobat PDF (3492 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Article Navigation

▸ Article Outline
– Abstract
1. Introduction
2. Phase-retrieval methods
3. Choosing your method
4. Numerical results
5. Discussion and conclusions
– References and links

Article Tools
Full Text
Article Info
References
Cited By
Figures
Metrics
Download PDF

Viewing Equations in the
HTML Article

Optimal Viewing Recommendations

Full Text
Article Info
References (39)
Cited By
Figures (9)
Metrics

Abstract

In-line phase-contrast X-ray imaging provides images where both absorption and refraction contribute. For quantitative analysis of these images, the phase needs to be retrieved numerically. There are many phase-retrieval methods available. Those suitable for phase-contrast tomography, i.e., non-iterative phase-retrieval methods that use only one image at each projection angle, all follow the same pattern though derived in different ways. We outline this pattern and use it to compare the methods to each other, considering only phase-retrieval performance and not the additional effects of tomographic reconstruction. We also outline derivations, approximations and assumptions, and show which methods are similar or identical and how they relate to each other. A simple scheme for choosing reconstruction method is presented, and numerical phase-retrieval performed for all methods.

1. Introduction

Classical medical three-dimensional X-ray imaging systems (computed tomography, CT) rely on absorption for contrast [1

See, e.g., P. Suetens Fundamentals of Medical Imaging (Cambridge Univ Press, 2009). [CrossRef]

]. Phase-sensitive detection shows promise for improved contrast where absorption is insufficient [2

R. Fitzgerald, “Phase-sensitive X-ray imaging,” Phys. Today 53(7), 23–26 (2000). [CrossRef]

, 3

R. A. Lewis, “Medical phase contrast X-ray imaging: current status and future prospects,” Phys. Med. Biol. 49, 3573–3583 (2004). [CrossRef] [PubMed]

] and, thus, for improved biomedical imaging. Unfortunately the simplest phase-contrast imaging arrangements do not automatically provide quantitative phase data suitable for tomographic reconstruction, making the application of phase-retrieval algorithms necessary. In the present paper we evaluate and compare the phase-retrieval methods applicable to in-line phase-contrast tomography and prove that they all, despite different origins, follow the same mathematical pattern. Thus, the applicability of the algorithms for different imaging situations can be compared in an unbiased way.

If the refractive index of a material is n = 1 − δ + iβ, the imaginary part β describes the absorption while the real part δ describes the phase shift introduced by the material. The phase is not directly observable and must be turned into phase contrast, i.e., intensity variations. Since the early publications on X-ray phase-contrast imaging [4

U. Bonse and M. Hart, “An X-ray interferometer,” Appl. Phys. Lett. 6, 155–156 (1965). [CrossRef]

–11

D. Chapman, W. Thomlinson, R. E. Johnston, D. Washburn, E. Pisano, N. Gmür, Z. Zhong, R. Menk, F. Arfelli, and D. Sayers, “Diffraction enhanced X-ray imaging,” Phys. Med. Biol. 42, 2015–2025 (1997). [CrossRef] [PubMed]

], different methods for doing this have been developed, as thoroughly reviewed by Nugent [12

K. A. Nugent, “Coherent methods in the X-ray sciences,” Adv. Phys. 59, 1–99 (2010). [CrossRef]

]. Currently four methods are in use: interferometry, diffraction-enhanced imaging, in-line phase contrast, and grating interferometry (Talbot imaging) [2

R. Fitzgerald, “Phase-sensitive X-ray imaging,” Phys. Today 53(7), 23–26 (2000). [CrossRef]

, 13

F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nat. Phys. 2, 258–261 (2006). [CrossRef]

, 14

A. Momose, “Recent advances in X-ray phase imaging,” Jpn. J. Appl. Phys. 44, 6355–6367 (2005). [CrossRef]

]. Out of those methods, the last two show distinct advantages. They can both be used with broad-band sources. A version of the grating interferometry (Talbot-Lau interferometry) can also be performed with extended sources as the demands on spatial coherence are lowered, at the cost of losing a portion of the flux. In-line phase contrast (in-line holography, propagation-based phase contrast) is the simplest arrangement and therefore the least sensitive to, e.g., misalignments of the optical components. It requires a small source of high spatial coherence and has therefore been limited to synchrotron or low-power microfocus sources, but has recently been demonstrated with table-top liquid-jet-anode sources of higher power [15

T. Tuohimaa, M. Otendal, and H. M. Hertz, “Phase-contrast X-ray imaging with a liquid-metal-jet-anode micro-focus source,” App. Phys. Lett. 91, 074104 (2007). [CrossRef]

]. None of the X-ray phase-contrast methods yield direct measurements of the phase, so data processing is required. The simplest is for interferometry, which only requires unwrapping of the 2π phase shifts. The most demanding is for in-line phase contrast which requires phase retrieval; the experimental simplicity is paid for by more complex data processing.

In this paper, we consider the method of in-line phase contrast, illustrated in Fig. 1. Because of its simplicity and relatively low sensitivity to misalignments it is suitable for tomography. The arrangement is similar to absorption X-ray imaging. The main difference is the propagation distance: in absorption, a projection through the object is captured at a plane close to the sample. In in-line phase contrast, the X-rays propagate a distance through free space after leaving the sample, before being registered by the detector. This distance gives time for the phase introduced by the sample to develop into phase contrast. If the propagation distance is chosen properly, the phase-contrast will consist of edge enhancement, as shown in the example in Fig. 2 where (a) is the object and (b) the simulated phase-contrast image. While the edges in (b) are clear and small structures can be detected, the quantitative relation between phase-contrast image and object is not intuitive. Fig. 2(c) shows instead the phase of the object in (a), a physical quantity directly related to object properties. Figure 2(c) has been obtained from Fig. 2(b) using phase retrieval, a method for reconstructing phase from measured intensity. Phase retrieval was first used for optical wavelengths [16

R. A. Gonsalves, “Phase retrieval by differential intensity measurements,” J. Opt. Soc. Am. A 4, 166–170 (1987). [CrossRef]

] and then brought into the X-ray regime when the need arose. The paper by Nugent [17

K. A. Nugent, “X-ray noninterferometric phase imaging: a unified picture,” J. Opt. Soc. Am. A 24, 536–546 (2007). [CrossRef]

] is an excellent summary of phase-retrieval techniques in X-ray imaging, while Gureyev et al. [18

T. E. Gureyev, S. C. Mayo, D. E. Myers, Ya. Nesterets, D. M. Paganin, A. Pogany, A. W. Stevenson, and S. W. Wilkins, “Refracting Röntgen’s rays: propagation-based X-ray phase contrast for biomedical imaging,” J. Appl. Phys. 105, 102005 (2009). [CrossRef]

] give a readable overview of in-line phase-contrast imaging, including practical considerations for designing a system. For in-line phase tomography, the phase is retrieved for all images, and three-dimensional tomographic reconstruction is performed on the retrieved phase. In this paper, we consider the phase retrieval but not the final step, i.e., that of three-dimensional reconstruction.

Fig. 1 Illustration of in-line phase contrast. Extending the propagation distance after passage through the object lets intensity differences due to refraction develop.

Fig. 2 (a) Phase contribution of an object consisting of polystyrene cylinders and spheres in air. (b) Simulated phase-contrast image of the object. (c) Phase retrieved from (b), using the single-material phase-retrieval method. Color scale is linear from minimum (black) to maximum (white).

In general, phase retrieval requires at least two measurements of the intensity, taken at two different distances from the source [17

K. A. Nugent, “X-ray noninterferometric phase imaging: a unified picture,” J. Opt. Soc. Am. A 24, 536–546 (2007). [CrossRef]

]. (The same effect can also be achieved by using different amounts of defocus in, e.g., phase-contrast microscopy [19

A. Barty, K. A. Nugent, D. Paganin, and A. Roberts, “Quantitative optical phase microscopy,” Opt. Lett. 23, 817–819 (1998). [CrossRef]

], or by taking images at different wavelengths.) This relates to the uniqueness of the retrieved phase. If only one image is captured, at a certain distance from the source, we could find a phase object with constant absorption that would give rise to the registered intensity. But we could also find a pure absorption object that would produce the exact same intensity. Considering that most objects have both varying absorption and phase, there are an infinite number of solutions to the problem. Consequently, the phase retrieved from only one image will not be unique. But the intensity generated by absorption and the intensity generated by phase will propagate differently, so if two images are taken at different distances, the phase and absorption properties can be untangled. As an example, a contact image taken just after the sample contains only the contributions from absorption, and once the absorption is known the phase information can be retrieved from an image taken further away from the object. So with two intensity measurements at different distances from the source, the retrieved phase will be unique. (Unique in a practical sense, as more exotic objects like phase vortices are not unambiguously recovered by phase retrieval [17

K. A. Nugent, “X-ray noninterferometric phase imaging: a unified picture,” J. Opt. Soc. Am. A 24, 536–546 (2007). [CrossRef]

].)

In tomography, hundreds of images will be taken in fast succession while either the object or the system is rotated. Taking two images, at two different distances and at each rotation angle, causes problems: it is difficult to arrange experimentally, it increases the dose delivered to the sample, and it increases the time spent making the measurements. The dose increment could be avoided by using only half the dose in each image, but this would instead increase the noise levels. Thus, taking only one image at each angle is preferable. In most cases, the investigated samples are not completely unknown, and prior knowledge can be used to reduce the number of required images. For example, in a homogeneous object (contains only one material and air) β and δ are known constants and the projected thickness can be found from only one image. The two most common assumptions on material properties are either that the absorption is constant and thus can be neglected, or that absorption and phase coefficients β and δ are proportional to each other. The second case will, for the remainder of this paper, be referred to in short as absorption proportional to phase.

In this paper, we review in-line phase retrieval methods applicable to tomography. First, that means we consider only methods that require one image at each projection angle. Second, we consider only analytical methods as we have hundreds of images to process, and iterative methods take longer. In the literature, we have identified seven methods that fulfil these criteria. (Sometimes identical or nearly identical methods have been derived in different ways by various authors, and are then considered as one method.) Surprisingly, despite being derived in different ways with different approximations and assumptions along the way, all these methods follow the same pattern and are numerically implemented in the same manner. This pattern is outlined in Sec. 2.1. This general pattern makes it very easy to compare them, and also to implement them numerically. As the methods differ only at specific points, the same numerical code can be used for all methods with an option for the user to choose the relevant method at each run of the program. In Sec. 2.2 we list and compare the assumptions and approximations made in the derivations, and thus identify the different situations where some methods should work better than others. Section 3 contains a scheme for choosing the most suitable phase-retrieval method in a particular situation, based on the information in Sec. 2.2. Finally, in Sec. 4 we test the methods on simulated and experimental phase-contrast images, and then in Sec. 5 discuss the results and conclude which methods are preferable in different situations.

The considered methods are (1) the Bronnikov [20

A. V. Bronnikov, “Reconstruction formulas for phase-contrast imaging,” Opt. Commun. 171, 239–244 (1999). [CrossRef]

] method that assumes no absorption, (2) the modified Bronnikov algorithm by Groso et al. [21

A. Groso, R. Abela, and M. Stampanoni, “Implementation of a fast method for high resolution phase contrast tomography,” Opt. Express 14, 8103–8110 (2006). [CrossRef] [PubMed]

] that allows for small absorption, (3) the phase-attenuation duality algorithm by Wu et al. [22

X. Wu and H. Liu, “X-Ray cone-beam phase tomography formulas based on phase-attenuation duality,” Opt. Express 13, 6000–6014 (2005). [CrossRef] [PubMed]

] for absorption proportional to phase, (4) the method for homogeneous object or single material by Paganin et al. [23

D. Paganin, S. C. Mayo, T. E. Gureyev, P. R. Wilkins, and S. W. Wilkins, “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” J. Microsc. 206, 33–40 (2002). [CrossRef] [PubMed]

], (5) the two-material method by Beltran et al. [24

M. A. Beltran, D. M. Paganin, K. Uesugi, and M. J. Kitchen, “2D and 3D X-ray phase retrieval of multi-material objects using a single defocus distance,” Opt. Express 18, 6423–6436 (2010). [CrossRef] [PubMed]

], and (6) the Fourier method with the Born approximation or (7) the Rytov approximation derived by Gureyev et al. [25

T. E. Gureyev, T. J. Davis, A. Pogany, S. C. Mayo, and S. W. Wilkins, “Optical phase retrieval by use of first Born- and Rytov-type approximations,” Appl. Opt. 43, 2418–2430 (2004). [CrossRef] [PubMed]

]. Number 6, the Fourier method with Born approximation, has also been derived by Zabler [26

S. Zabler, P. Cloetens, J.-P. Guigay, and J. Baruchel, “Optimization of phase contrast imaging using hard x-rays,” Rev. Sci. Intstrum. 76, 073705 (2005). [CrossRef]

], Guigay [27

J. P. Guigay, “Fourier transform analysis of Fresnel diffraction patterns and in-line holograms,” Optik 49, 121–125 (1977).

], and Turner [28

L. D. Turner, B. B. Dhal, J. P. Hayes, A. P. Mancuso, K. A. Nugent, D. Paterson, R. E. Scholten, C. Q. Tran, and A. G. Peele, “X-ray phase imaging: Demonstration of extended conditions with homogeneous objects,” Opt. Express 12, 2960–2965 (2004). [CrossRef] [PubMed]

] under somewhat different approximations.

2. Phase-retrieval methods

2.1. General pattern

All seven methods listed in the introduction follow the same pattern. The steps of this pattern are outlined in Fig. 3, assuming the object is illuminated by a plane wave. (At the end of Sec. 2.2 it is shown how to extend the theory to the more practical cone-beam case, where the source is placed at a finite distance from the object and the illumination is a spherical wave.) As input, the algorithms take the intensity I(r _⊥) in the image plane registered as a function of transverse position r _⊥ = (x, y) on the detector. First, a function g(I(r _⊥)) is calculated. This function, which varies between the different methods, is rather simple and sometimes just a normalization. Then the quantity g(r _⊥) is filtered in the frequency domain. Intuitively, this filtering can be perceived as a deconvolution of a diffraction integral (like the Guigay equation [27

J. P. Guigay, “Fourier transform analysis of Fresnel diffraction patterns and in-line holograms,” Optik 49, 121–125 (1977).

]) or as a Fourier transform solution of a wave equation (like the transport of intensity equation [9

K. A. Nugent, T. E. Gureyev, D. F. Cookson, D. Paganin, and Z. Barnea, “Quantitative phase imaging using hard X-rays,” Phys. Rev. Lett. 77, 2961–2964 (1996). [CrossRef] [PubMed]

, 29

M. R. Teague, “Deterministic phase retrieval: a Green’s function,” J. Opt. Soc. Am. 73, 1434–1441 (1983). [CrossRef]

]). The filtering is performed by first taking the Fourier transform of g(I), multiplying by a filter function H _p(w) where w = (u, v) is the spatial frequency, and then taking the inverse Fourier transform to get the filtered quantity g_F (r _⊥). Again, the filter function H _p(w) depends on the chosen method. Finally, a function f (g_F ) is taken to yield the 2D phase distribution φ(r _⊥) at a plane just after the object, i.e., at the contact plane. The function f (g_F ) is rather simple and often an identity, f (g_F ) = g_F . The procedure can be written as one equation,

φ (r_{⊥}) = f (ℱ^{- 1} {H_{p} \cdot ℱ [g (I)]})

(1)

where ℱ denotes the 2D Fourier transform with respect to r _⊥.

Fig. 3 The process of phase retrieval followed by tomographic reconstruction.

As shown in Fig. 3, the process can be continued to retrieve not just the 2D projections of the phase, but the 3D phase distribution of the object. Then for each tomographic angle, the 2D phase is filtered again in the Fourier plane, this time by a function H _fbp(w). For the standard method of 3D reconstruction, namely filtered backprojection, H _fbp(w) = |u|. Other methods are in use, in particular some methods where the filter has been tailored specifically to suit phase-contrast images [20

A. V. Bronnikov, “Reconstruction formulas for phase-contrast imaging,” Opt. Commun. 171, 239–244 (1999). [CrossRef]

, 22

X. Wu and H. Liu, “X-Ray cone-beam phase tomography formulas based on phase-attenuation duality,” Opt. Express 13, 6000–6014 (2005). [CrossRef] [PubMed]

]. After the filtering, the image is back-projected onto the 3D volume of interest [30

A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (Siam, 2001). [CrossRef]

The process is identical for all seven phase-retrieval methods, except that different functions g(I), H _p(w), f (g_F ), and sometimes H _fbp(w) are used. Each method is completely described by those four functions, and by comparing the functions we can compare the methods. This is done in Sec. 2.2.

In Fig. 3 we see that before and after taking the function f (g_F ), there is a Fourier transform and an inverse Fourier transform. If f (g_F ) = g_F or some other very basic function, the two transforms need not be performed. Then the two filters H _p(w) and H _fbp(w) can be combined into a single filter H _tot(w), and the procedure in Fig. 3 can be simplified into that of Fig. 4. The new process starts by taking the function g(I(r)) for each tomographic angle, followed by a Fourier transform and multiplication in frequency space by the combined filter H _tot(w). After an inverse Fourier transform, each image is back-projected to the 3D volume. After processing images at all angles, a 3D phase-retrieved volume has been reconstructed. This reconstruction is hardly more time-consuming than filtered back-projection alone. The first combined method was derived by Bronnikov [20

A. V. Bronnikov, “Reconstruction formulas for phase-contrast imaging,” Opt. Commun. 171, 239–244 (1999). [CrossRef]

Fig. 4 The process of phase retrieval and tomographic reconstruction performed together.

While the combined approach in Fig. 4 efficiently saves computational time, the process in Fig. 3 has the advantage of flexibility. For example, several different 3D reconstruction techniques can be tested on the same set of phase-retrieved images, without having to re-calculate the phase retrieval in between. It is possible to choose a custom-made algorithm for the phase retrieval, and combine it with commercial, fast software for the 3D reconstruction. We also note that some phase-retrieval methods are suitable for the combined approach in Fig. 4, while others are not. A slightly more complicated function f (g_F ), containing, e.g., a logarithm, makes the combined approach impossible to use.

As the intention of this paper is to compare phase retrieval methods, we will concentrate on the first three functions g(I), H _p(w), and f (g_F ). While comparing different H _fbp would be interesting, it lies outside the scope of the present survey. For methods that combine phase retrieval and 3D reconstruction, we have done the mathematics necessary to extract only the phase-retrieval part as described by the three functions.

2.2. Comparison of different methods

In the derivations of the seven phase-retrieval methods, different assumptions have been made. Some concern the material, while others deal with the wave propagation. One common assumption is to start from the Fresnel diffraction integral, another to use the paraxial (or parabolical) wave equation. As the Fresnel diffraction integral is actually the solution of the paraxial wave equation [31

R. Grella, “Fresnel propagation and diffraction and paraxial wave equation,” J. Opt. (Paris) 13, 367–364 (1982).

], these two assumptions are the same and will be referred to as the Fresnel approximation. Sometimes the Transport of Intensity Equation (TIE) [29

M. R. Teague, “Deterministic phase retrieval: a Green’s function,” J. Opt. Soc. Am. 73, 1434–1441 (1983). [CrossRef]

] is used. It is derived from the paraxial wave equation, with an additional assumption identical to assuming a large Fresnel number, N_F = a ² /λd ≫ 1 [32

T. E. Gureyev, Y. I. Nesterets, A. W. Stevenson, P. R. Miller, A. Pogany, and S. W. Stevenson, “Some simple rules for contrast, signal-to-noise and resolution in in-line phase-contrast imaging,” Opt. Express 16, 3223–3241 (2008). [CrossRef] [PubMed]

]. This expression relates the smallest feature size a of the object to the wavelength λ and the propagation distance d, and holds for relatively small propagation distances. Another approximation used in several derivations is the slowly varying phase (SVP) approximation, also referred to as moderate phase variation, |φ (r _⊥ + λd w _m ) – φ(r _⊥)| ≪ 1 which implies that the changes in phase in the object are not too rapid. The frequency |w _m | is the highest spatial frequency of the object, or at least the highest frequency to be included in the reconstruction.

For all methods, the object is characterized by its three-dimensional refractive index n(r) = 1 – δ(r) + iβ (r), where δ represents the phase and β the absorption, and r is position in three dimensions. The quantity β is related to the linear absorption coefficient μ by μ = 4π/λ · β [33

D. M. Paganin, Coherent X-Ray Optics (Oxford Science Publications, 2006). [CrossRef]

]. Under the projection approximation, the intensity I(r _⊥) in the contact plane is given by [33

D. M. Paganin, Coherent X-Ray Optics (Oxford Science Publications, 2006). [CrossRef]

]

I (r_{⊥}) = I_{in} exp [- \int d z μ (r_{⊥}, z)]

(2)

where r _⊥ is the coordinates in the projection plane, perpendicular to the projection direction, z the coordinate along the projection direction, and the incident intensity I _in is assumed constant. If the material is homogeneous, this turns into Beer’s law I(r _⊥) = I _in exp[–μT(r _⊥)] where T(r _⊥) is the projected thickness of the material. The phase is similarly given by

φ (r_{⊥}) = - \frac{2 π}{λ} \int d z δ (r_{⊥}, z)

(3)

which for homogeneous objects becomes φ(r _⊥) = −δT(r _⊥) · 2π/λ. In all derivations except the one for two materials, one out of two approximations is made: either the material has no absorption, μ(r _⊥) = 0, or the absorption is proportional to the phase, μ ∝ δ.

Below follows a short description of each method and its derivation. The functions g(I), H _p(w), and f (g_F ) of all methods are given in Table 1, while the different approximations are listed in Table 2. For the methods of single material and two materials, the final result will be given as projected thickness of the sample, while for the others it is given as phase. For the single-material method, the result can be turned into phase using the simplified version of Eq. (3). For the same two methods [23

, 24

], our expressions differ slightly from the ones given in the paper, as we use a different definition of the Fourier transform. We define the Fourier transform f̂ (w) of a function f (r _⊥) as

\hat{f} (w) = \int \int_{- \infty}^{\infty} d^{2} r_{⊥} f (r_{⊥}) exp (- 2 π r_{⊥} \cdot w)

(4)

and its inverse as

f (r_{⊥}) = \int \int_{- \infty}^{\infty} d^{2} u \hat{f} (w) exp (2 π w \cdot r_{⊥}) .

(5)

This definition matches the common definition of the Discrete Fourier Transform (DFT) and thus of the Fast Fourier Transform (FFT) which is generally used for numerical evaluation of the transforms.

Table 1 Methods of Phase Retrieval Suitable for In-Line Phase-Contrast Tomography, and Their Properties ¹

Method		g(I)	H _p(u, v)	f (g_F )	Derivation
Bronnikov		$\frac{I}{I_{in}} - 1$	[2πλd\|w\|²]⁻¹	g_F	Assume μ = 0. Use weak focusing condition.
Modified Bronnikov		$\frac{I}{I_{in}} - 1$	[2πλd\|w\|² + α]⁻¹	g_F	Assume μ ≈ 0. Solve TIE.
Phase-atten. duality		$\frac{I}{I_{in}}$	${[2 π \frac{r_{e} λ^{2} d}{σ_{KN}} \| w \|^{2} + 1]}^{- 1}$	$\frac{λ r_{e}}{σ_{KN}} ln g_{F}$	Assume μ ∝ δ. Propagate Wigner distribution function for small Fresnel numbers.
Single material		$\frac{I}{I_{in}}$	${[4 π^{2} d \frac{δ}{μ} \| w \|^{2} + 1]}^{- 1}$	$- \frac{1}{μ} ln g_{F}$	Assume known δ and μ (implies μ ∝ δ). Solve TIE.
Two materials		$\frac{I exp [μ_{2} A (r_{⊥})]}{I_{in}}$	${[4 π^{2} d \frac{δ_{1} - δ_{2}}{μ_{1} - μ_{2}} \| w \|^{2} + 1]}^{- 1}$	$- \frac{1}{μ_{1} - μ_{2}} ln g_{F}$	Assume two materials and SVP. Solve TIE.
Fourier method, Born type	a	$\frac{1}{2} (\frac{I}{I_{in}} - 1)$	[sin (πλd\|w\|²)]⁻¹	g_F	Use Born approximation and Fresnel propagation. Finally, assume (a) μ = 0 or (b) μ ∝ δ.
Fourier method, Born type	b	$\frac{1}{2} (\frac{I}{I_{in}} - 1)$	[γ cos (πλd\|w\|²) + sin (πλd\|w\|²)]⁻¹	g_F
Fourier method, Rytov type	a	$\frac{1}{2} ln \frac{I}{I_{in}}$	[sin (πλd\|w\|²)]⁻¹	g_F	Use Rytov approximation and Fresnel propagation. Finally, assume μ = 0 (a) or μ ∝ δ (b).
Fourier method, Rytov type	b	$\frac{1}{2} ln \frac{I}{I_{in}}$	[γ cos (πλd\|w\|²) + sin (πλd\|w\|²)]⁻¹	g_F

¹ The functions g(I), H _p(u, v) and f (g_F ) are introduced in Sec. 2.1, and give a full description of the methods.

Table 2 Approximations Made in the Derivation of the Methods of Table 1 ¹

Method		μ ≈ 0	μ ∝ δ	2 mat.	Fresnel	$\frac{a^{2}}{d λ} ≫ 1$	SVP	Born	Rytov
Bronnikov		✓			✓	✓
Modified Bronnikov		✓			✓	✓
Phase-att. duality			✓		✓	✓	✓
Single material			✓		✓	✓
Two materials				✓	✓	✓	✓
Fourier (Born)	a	✓			✓			✓
Fourier (Born)	b		✓		✓			✓
Fourier (Rytov)	a	✓			✓				✓
Fourier (Rytov)	b		✓		✓				✓

¹ The first three approximations apply to the material, which is assumed to have no absorption (μ ≈ 0), to have absorption proportional to phase (μ ∝ δ), or to consist of two known materials and air. The following apply to the method of propagation and are, in order: the Fresnel approximation, the assumption that the Fresnel number is large, the slowly varying phase approximation (SVP), the first Born approximation, and the first Rytov approximation.

Bronnikov [20

A. V. Bronnikov, “Reconstruction formulas for phase-contrast imaging,” Opt. Commun. 171, 239–244 (1999). [CrossRef]

]. This was the first method derived. The absorption is assumed to be zero, and weak focusing conditions [34

J. M. Cowley, Diffraction Physics (Elsevier, 1995).

] applied. The weak focusing condition formula is similar to the TIE, and is derived from the Fresnel diffraction integral for large Fresnel numbers. The derivation goes all the way to the 3D reconstruction, so to get the intermediate step of the phase we have applied Fourier methods to Eq. (5) of Ref. [20

A. V. Bronnikov, “Reconstruction formulas for phase-contrast imaging,” Opt. Commun. 171, 239–244 (1999). [CrossRef]

]. The method is limited to very thin samples, because of the assumption of no absorption.

Modified Bronnikov [21

A. Groso, R. Abela, and M. Stampanoni, “Implementation of a fast method for high resolution phase contrast tomography,” Opt. Express 14, 8103–8110 (2006). [CrossRef] [PubMed]

]. The object is assumed to be of weak and almost homogeneous absorption, and the TIE is applied and solved. A result identical to Bronnikov’s is derived, and afterwards a small term α is added to the filter function to account for what little absorption there might be. We note that although the method of derivation differs, the initial assumptions are the same as for the Bronnikov method: no absorption, Fresnel approximation, and large Fresnel number. The absorption requirement is relaxed when α is introduced.

Phase-attenuation duality [22

X. Wu and H. Liu, “X-Ray cone-beam phase tomography formulas based on phase-attenuation duality,” Opt. Express 13, 6000–6014 (2005). [CrossRef] [PubMed]

]. It is assumed that the phase-attenuation duality applies, where phase and absorption are both proportional to the electron density and thus proportional to each other, μ ∝ δ. This holds where Compton scattering is the main contributor to μ, such as for light materials and photon energies of 60 – 500 keV. Additionally, the SVP approximation is made. Then the Wigner distribution function is found and propagated to yield a solution for large fresnel numbers. The propagation of the Wigner function is valid in the Fresnel approximation [35

M. J. Bastiaans, “Application of the Wigner distribution function to partially coherent light,” J. Opt. Soc. Am. A 3, 1227–1237 (1986). [CrossRef]

], so the applied approximations are absorption proportional to phase, Fresnel approximation, large Fresnel number, and SVP approximation. The proportionality constant between phase and absorption is δ/μ = λ ² r_e /2πσ_KN where λ is the wavelength, r_e the classical electron radius, and σ_KN the total cross section for X-ray photon Compton scattering from a single free electron [22

X. Wu and H. Liu, “X-Ray cone-beam phase tomography formulas based on phase-attenuation duality,” Opt. Express 13, 6000–6014 (2005). [CrossRef] [PubMed]

Single material [23

]. It is assumed that the object is homogeneous, and that δ and μ are known. This assumption implies that μ ∝ δ. Under this condition the TIE is solved to yield the thickness of the material. Thus, the assumptions are absorption proportional to phase, Fresnel approximation, and large Fresnel number.

Two materials [24

]. The object is assumed to consist of air and two other materials, one embedded in the other, of known β ₁ and δ ₁ (embedded material) and β ₂ and δ ₂ (encasing material). It is also assumed that the total projected thickness of the object A(r _⊥) (i.e., the sum of the projected thickness of the two materials) is known. The TIE is solved assuming that the thickness of the encasing material varies slowly, an assumption very similar to that of slowly varying phase. So the assumptions are two known materials of known total thickness, the Fresnel approximation, large Fresnel number, and slowly varying phase.

Fourier method with Born approximation [25

]. Starting form the Fresnel diffraction integral, it is assumed that the Born approximation applies. In this context it means the scattered field ψ(r _⊥) is small both in amplitude and phase, |ψ (r _⊥)| ≪ 1. The Fresnel integral is propagated to yield the results. Finally, either μ = 0 or μ = γδ· 4π/λ is assumed, where γ is a proportionality constant. The assumptions are thus the Fresnel approximation, the Born approximation, and either no absorption or absorption proportional to phase. In the case of no absorption, the Born approximation is not needed and the same expression can be derived assuming only the Fresnel approximation and slowly varying phase (SVP) [26

S. Zabler, P. Cloetens, J.-P. Guigay, and J. Baruchel, “Optimization of phase contrast imaging using hard x-rays,” Rev. Sci. Intstrum. 76, 073705 (2005). [CrossRef]

]. Similarly if the absorption is proportional to the phase, the same expression can be derived under less stringent assumptions [27

J. P. Guigay, “Fourier transform analysis of Fresnel diffraction patterns and in-line holograms,” Optik 49, 121–125 (1977).

, 28

], namely Fresnel approximation, homogeneous object, weak absorption, and slowly varying phase (SVP).

Fourier method with Rytov approximation [25

]. Starting from the Fresnel diffraction integral, it is assumed that the Rytov approximation applies,

| \nabla_{⊥} ψ (r_{⊥}) |^{2} ≪ | \nabla_{⊥}^{2} ψ (r_{⊥}) |

. This approximation is often equivalent to the Born approximation [25

]. The Fresnel diffraction integral is propagated to yield the results. Finally, either μ = 0 or μ = γδ · 4π/λ is assumed. The assumptions are thus the Fresnel approximation, the Rytov approximation, and either no absorption or absorption proportional to phase.

Comparing the results of the derivations, as shown in Table 1, shows for example that if the somewhat arbitrary constant α of the modified Bronnikov method is set to zero, it is identical to the Bronnikov method. Whereas if the constant is set to λ/2π · μ/δ, it turns instead into something very close to the single-material method. For a homogeneous object at high energies, δ/μ = λ ² r_e /2πσ_KN [22

X. Wu and H. Liu, “X-Ray cone-beam phase tomography formulas based on phase-attenuation duality,” Opt. Express 13, 6000–6014 (2005). [CrossRef] [PubMed]

] and the phase-attenuation duality method is identical to the single-material method. If the approximation of a large Fresnel number is used on the two Fourier methods, they turn into something very similar to the single-material method, and for the Born approximation with no absorption this yields the Bronnikov method. The Fourier methods in the Born and Rytov approximations become identical if the contrast is low enough that ln(I/I ₀) ≈ I/I ₀ − 1 [25

]. In short, the seven methods are similar and can sometimes be turned into each other by choices of parameters, approximations, or just a change of notation.

In Table 2, the approximations of the different methods are listed. First, we note that all methods apply the Fresnel approximation. This is reasonable, since development of phase contrast requires a certain propagation distance, so it is highly unlikely that the Fresnel approximation would not be valid. The possible exception would be for large cone-beam angles. All methods also assume the projection approximation of Eqs. (2) and (3), although this is not noted in the table. We can also see that the first five methods are restricted to shorter propagation distances as the Fresnel number is high, whereas the last two can be used for longer propagation distances. The first assumes no absorption at all, while the others in some way incorporate absorption.

As listed in Table 1, all methods assume monochromatic light and point X-ray sources. Some methods, however, incorporate the effects of polychromatic light and extended sources (phase-attenuation duality [22

X. Wu and H. Liu, “X-Ray cone-beam phase tomography formulas based on phase-attenuation duality,” Opt. Express 13, 6000–6014 (2005). [CrossRef] [PubMed]

, 36

X. Wu and H. Liu, “A new theory of phase-contrast X-ray imaging based on Wigner distributions,” Med. Phys. 31, 2378–2384 (2004). [CrossRef] [PubMed]

] and the single-material method [37

G. R. Myers, S. C. Mayo, T. E. Gureyev, D. M. Paganin, and S. W. Wilkins, “Polychromatic cone-beam phase-contrast tomography,” Phys. Rev. A 76, 045804 (2007). [CrossRef]

]). Here, we have chosen a simplified version of these methods, for easy comparison to the others.

For the Bronnikov method and both Fourier methods, the denominator of the filter function sometimes goes to zero. This implies noise amplification at these frequencies. To amend this, Gureyev et al. [25

] suggest a Tikhonov’s regularization term for the Fourier methods. Also known as a Wiener filter [38

J. W. Goodman, Introduction to Fourier Optics , 2nd ed. (McGraw-Hill, 1996).

], the principle is that the filter H _p(w) = 1/h _p(w) is replaced by

H_{wien} (w) = \frac{h_{p}^{*} (w)}{h_{p} (w) h_{p}^{*} (w) + η}

(6)

where the constant η can be adjusted depending on noise level. A large η gives small noise amplification, but also reduces the effect of the phase retrieval. A small η gives a better reconstruction of the object, but more noise amplification. The regularization has been implemented for the Fourier methods [25

], but the Bronnikov method would also benefit from it. In fact, the α introduced in the modified Bronnikov method can be seen as a kind of regularization.

The formulas in Table 1 are given for illumination by a plane wave, and propagation distance d. In most practical cases, the illumination will come from a point source (or at least approximately so) placed a distance R ₁ from the object, and the image will be captured a distance R ₂ behind the object. The effect on the reconstruction formulas is readily included, since the intensity I _{R
₁} captured in the cone-beam case is related to the intensity I _∞ generated by a plane wave by [23

]

I_{R_{1}} (M r_{⊥}, R_{2}) = \frac{1}{M^{2}} I_{\infty} (r_{⊥}, \frac{R_{2}}{M}),

(7)

where M = (R ₁ + R ₂)/R ₁ is the magnification of the system. For incident plane waves, the transverse coordinate r _⊥ was the shared coordinate of object and detector plane both, whereas now r _⊥ is the object-plane transverse coordinate and the detector plane coordinate is M r _⊥. For our case, Eq. (7) implies three changes to the equations in Table 1. First, the intensity measured at distance R ₂ should be multiplied by M ². Second, the intensity I(r _⊥) in the formulas should be replaced by the measured intensity I _{R
₁} (M r _⊥). Third, and most important, the propagation distance d in Table 1 should be replaced by the effective propagation distance R ₂/M. If those changes are applied, all the methods in Table 1 can be used for phase retrieval in cone-beam geometry. We note that the first two changes have rather minor effects on the results. First, the incident intensity I _in is often measured at a distance R ₂, as a background image with no object in it. Then the same factor M ² will be applied to both I(r _⊥) and I _in, and since all methods use I(r _⊥)/I _in, the effect of this change will cancel. Second, the Fourier transform should still be taken with respect to r _⊥ and not with respect to M r _⊥. Consequently, this change has marginal effects on the numerical calculations as long as the scale M r _⊥ of the detector plane is downscaled by M before given as an object coordinate.

3. Choosing your method

Assuming a point source and monochromatic light, you can choose your method of reconstruction following the brief scheme in Fig. 5. First, consider the distance at which your images are registered. If the Fresnel number is much smaller than one you must go with one of the Fourier methods, using either the Born or Rytov approximation. The two methods give fairly similar results, but the Born approximation is slightly easier to implement as the Rytov approximation contains a logarithm in g(I), which causes problems for zero intensity. The Fourier methods tend to amplify noise at specific frequencies, as can be seen from the zeros of the denominator of H _p(w). For short distances (Fresnel numbers around or larger than 1) one of the other five methods will give equivalent or better results. Also, none of the methods apply unless μ ∝ δ or there is no absorption. For samples that do not fulfil one of these requirements there is no suitable reconstruction method, and your best long-distance option is to try one of the Fourier methods anyway.

Fig. 5 Procedure for choosing your phase-retrieval method.

For short distances, i.e., large Fresnel numbers, there are five methods remaining, and three of them are reserved for special cases. First, if the beam energy is high, namely 60 – 500 keV, the phase-attenuation duality method is superior to any of the others. All the other methods place demands on the material, for example that it is homogeneous. For phase-attenuation duality, the high beam energy means μ is automatically proportional to δ for all light materials, and thus the method works for homogeneous or multi-material samples alike. Second, if the sample consists of two known materials and air, and you know the total projected thickness of the sample, the two-material method should be used. Third, if both noise and absorption are very low, the Bronnikov method is the simplest. However, since it is based on the assumption of no absorption at all, it tends to produce artifacts for most objects. It is therefore suggested that the modified Bronnikov method, which is almost as easy to implement, is used instead.

For all remaining cases, i.e., short propagation distance, lower beam energy (< 60 keV), and non-zero absorption, either the single material method or the modified Bronnikov method are the best options. They are very similar, with the main difference in the f (g) function, which implies they give somewhat different reconstructed profiles. The single-material method requires known material constants, but once those are known, the method is deterministic. The modified Bronnikov method involves a constant α which is normally decided from trial and error. Comparing the two methods suggest a good starting point is α = λ/2π · μ/δ. The single-material method is formally limited to homogeneous objects, e.g., samples that consist of one material and air. The derivation of the modified Bronnikov method is not entirely clear on this point, but its similarity to the single-material method implies that the same limit applies. In practice, though, both methods are used for multi-material samples as well, with different kinds of artifacts showing up. The artifacts consist of blurring, if the real μ/δ is higher than the value chosen in the reconstruction, or edge amplification, if the real μ/δ is lower than the value used. The method is used despite those artifacts, simply because there is no method to handle multi-material objects for lower beam energies. Multi-material objects imply a more complicated relationship between absorption and phase, and would require more than one image at each tomographic angle for proper phase retrieval.

Finally, situations where the effects of extended source and/or polychromatic light are significant lie outside the scope of the present survey. As extended versions of the phase-attenuation duality [22

X. Wu and H. Liu, “X-Ray cone-beam phase tomography formulas based on phase-attenuation duality,” Opt. Express 13, 6000–6014 (2005). [CrossRef] [PubMed]

, 36

X. Wu and H. Liu, “A new theory of phase-contrast X-ray imaging based on Wigner distributions,” Med. Phys. 31, 2378–2384 (2004). [CrossRef] [PubMed]

] and single-material [37

G. R. Myers, S. C. Mayo, T. E. Gureyev, D. M. Paganin, and S. W. Wilkins, “Polychromatic cone-beam phase-contrast tomography,” Phys. Rev. A 76, 045804 (2007). [CrossRef]

] methods contain tools to handle this situation, they are a logical starting point.

4. Numerical results

Numerical calculations and simulations have been performed to illustrate the phase-retrieval methods. First, the object in Fig. 2(a) was numerically propagated to the detector plane by numerical evaluation of the Fresnel diffraction integral, to yield an image similar to Fig. 2(b) except it is now free from noise. Figure 6 shows the result when six different phase-retrieval methods were used. The object consists of polystyrene cylinders and spheres in air, with polystyrene parameters taken as β = 3.553 · 10⁻¹⁰ and δ = 1.043 · 10⁻⁶. The diameters are 100, 50, 20, and 10 μm respectively, leading to a maximum absorption of 0.5% and maximum phase of 7.9 radians for the thickest part of the object. The radiation is assumed monochromatic at 15keV and the source Gaussian of full width at half maximum (FWHM) 10 μm. The simulated detector pixels are 9 μm and the detector point-spread function (PSF) is Gaussian of FWHM 25 μm. The distance from source to object is R ₁ = 0.6m and the distance from object to detector is R ₂ = 2.4m, yielding a magnification of M = 5. The resulting phase-retrieved images, displayed in Fig. 6, show that most methods give very similar results in the absence of noise. Two of them stand out: the Bronnikov method in Fig. 6(a) assumes no absorption, and since the object has non-uniform absorption this causes a brighter background in the left half of the image. The phase-attenuation duality image in Fig. 6(c) is not as sharp as the others, and the smallest structures are not reconstructed, since the method is not suited for this photon energy.

Fig. 6 Retrieved phase of the object in Fig. 2(a) from simulated noise-free phase-contrast images, for R ₁ = 0.6 m, R ₂ = 2.4 m, photon energy 15 keV, and simulated detector pixel size 9 μm. The phase is retrieved using (a) the Bronnikov method, (b) the modified Bronnikov method for α = 1.0 · 10⁻³, (c) the phase-attenuation duality, (d) the single-material method, (e) the Fourier-Born method using γ = 5.0 · 10⁻⁴, and (f) the Fourier-Rytov method using γ = 5.0 · 10⁻⁴. All color scales are linear ranging from −6 to 13 radians for (a), −2 to 7 radians for (b) and (d)–(f), and −3 to 10 radians for (c). Part (g) shows line profiles, taken along the white line in figures (a)–(f), for all six methods.

Figure 7 shows the same situation as Fig. 6, except that noise is now included at a pixel signal-to-noise ratio (SNR) of 4. Figure 7(a) now displays the phase-contrast image on the detector, while (b)–(f) show the results of five phase-retrieval methods. The phase-attenuation duality is not included, since there is no point in comparing it to the others in the wrong wavelength region. Here the differences between methods are more obvious. The modified Bronnikov and the single-material methods give the best results, while the Bronnikov method shows the same error as in Fig. 6. The two Fourier methods show amplification of noise at particular frequencies, despite the use of regularization to improve the situation.

Fig. 7 Retrieved phase of the object in Fig. 2(a) from simulated images (a) with a pixel SNR of 4, under the same conditions as Fig. 6. The phase is retrieved using (b) the Bronnikov method, (c) the modified Bronnikov method for α = 1.0 · 10⁻³, (d) the single-material method, (e) the Fourier-Born method using γ = 5.0·10⁻⁴, and (f) the Fourier-Rytov method using γ = 5.0 · 10⁻⁴. For the Fourier methods a Tikhonov’s regularization term η = 10⁻⁶ was used. All color scales are linear ranging from 0.3 to 2 in normalized pixel intensity for (a), −14 to 22 radians in (b), and −7 to 11 radians in (c)–(f). Part (g) shows line profiles, taken along the white line in figures (b)–(f), for all five methods.

Figure 8 shows some special cases. In Fig. 8(a) we see the same retrieved image as in Figs. 6 and 7, but this time using radiation of photon energy 100 keV. Phase retrieval was performed using the phase-attenuation duality method, which at this photon energy gives a correct relationship between μ and δ. The phase is now properly reconstructed. In order to judge the quality of the phase-attenuation duality method, Fig. 8(a) should be used in place of Fig. 6(c).

Fig. 8 Some special cases of phase retrieval. (a) Retrieved phase of the object in Fig. 2 from a simulated noise-free phase-contrast image at 100 keV, using the phase-attenuation duality method. (b) Simulated noise-free phase-contrast image of an object consisting of a 400 μm square rod of PMMA containing 100 μm diameter spheres of water (upper), teflon (middle), and air (lower). (c) Phase retrieved from the image in (b), assuming the encasing material is PMMA and the material of interest is teflon, using the two-material method. (d) Simulated noise-free phase-contrast image of a 20 μm cylinder at R ₁ = 6m and R ₂ = 24m. (e) Phase retrieved from (d) using the Fourier-Rytov method (γ = 5.0 · 10⁻⁴) with regularization term η = 10⁻². (f) Same as (e), except the single-material method is used for phase retrieval.

Figure 8(b) shows a simulated phase-contrast image of an object consisting of more than one material, namely a 400 μm square rod of PMMA (β = 6.459 · 10⁻¹⁰, δ = 1.186 · 10⁻⁶) containing three 100 μm spheres of water (upper, β = 8.968 · 10⁻¹⁰, δ = 1.026 · 10⁻⁶), teflon (middle, β = 2.591 · 10⁻⁹, δ = 1.953 · 10⁻⁶), and air (lower). All parameters are the same as in Fig. 6. Figure 8(c) contains the phase retrieved using the two-material method, assuming the encasing material is PMMA and the embedded material is teflon. The total thickness of the sample is assumed to be the thickness of the rod, 400 μm. The two-material method attempts to reconstruct the thickness of the embedded material, in this case teflon. Ideally, the retrieved image should contain only the sphere in the middle of the image. The water sphere causes artifacts because the method is not adjusted to this material. The air bubble causes similar effects, both in absorption and edge amplification, because it was not included in the total thickness of the sample (it is reasonable to assume that the exact size and location of air bubbles inside a sample are unknown). The edges of the rod, finally, show up because of the approximations in the derivation of the method: the encasing material is assumed to be slowly varying. This implies the method can handle the difference in absorption between the interior and exterior of the rod, so the background color of the retrieved image is uniform. However, it cannot handle the rapid change at the edge of the rod, which shows up as an amplified edge in the retrieved image.

Figure 8(d) shows a simulated noise-free phase-contrast image of a single cylinder of polystyrene, of diameter 20 μm, for a longer propagation distance (R ₁ = 6m, R ₂ = 24m). If the smallest detail is 10 μm this corresponds to a Fresnel number of N _F ≈ 0.2. The radiation is monochromatic at photon energy 15keV, the source has a FWHM of 2 μm, the detector PSF is a Dirac delta function, and the pixel size is 3 μm. Due to the longer propagation distance, the edge amplification has turned into a series of oscillations. Figure 8(e) contains the phase retrieved using the Fourier method in the Rytov approximation. The image displays both amplification of specific frequencies seen in the background, and a fairly sharp image of a cylinder of correct width. Figure 8(f) was retrieved using the single-material method, and shows the cylinder as blurred and of the wrong diameter, but without the background effects.

Figure 9 contains phase retrieval on experimental data. Figure 9(a) shows blood vessels in an extracted rat kidney using CO₂ as contrast medium, obtained using a broad-band Galinstan-based liquid-jet microfocus X-ray source operated at 50 kVp [39

U. Lundström, P. A. C. Takman, L. Scott, H. Brismar, and H. M. Hertz, “Low-dose high-resolution laboratory phase-contrast X-ray imaging,” manuscript in preparation.

]. The source-to-object distance is R ₁ = 0.6 m and the object-to-detector distance R ₂ = 2.4 m, the photon energy is centered around 15 keV, and the detector pixel size is 9 μm. The average dose is 100mGy and the kidney around 7 mm thick. For the reconstruction, the materials are considered as air and soft tissue (β = 9.73 · 10⁻¹⁰, δ = 1.077 · 10⁻⁶). Figure 9(b)-(f) contains the retrieved phase using different phase-retrieval methods. The Bronnikov method shows a clear disadvantage, as it magnifies low-frequency noise to such an extent that the entire image is blurred. The modified Bronnikov method and the single-material method both give good results. As the noise level is relatively low, the two Fourier methods using regularization give nearly the same quality as the other two.

Fig. 9 Phase retrieval on experimental data, in this case blood vessels in a rat kidney using CO₂ as contrast medium. (a) Phase-contrast image of the blood vessels, taken at source-to-object distance R ₁ = 0.6 m and object-to-detector distance R ₂ = 2.4 m at a photon energy centered at around 15 keV, using a detector of pixel size 9 μm. (b) Phase retrieved using the Bronnikov method. (c) Phase retrieved using the modified Bronnikov method for α = 2.1 · 10⁻³. (d) Phase retrieved using the single-material method. (e) Phase retrieved using the Fourier method in the Born approximation, for γ = 1.0 · 10⁻³ and regularization parameter η = 1 · 10⁻⁶. (f) Same as (e), except in the Rytov approximation. All color scales are linear ranging from 0.5 to 1.3 in normalized pixel intensity for (a), −70 to 70 radians in (b), and −7 to 3 radians in (c)–(f).

Table 3 shows the least mean square error for the normalized retrieved phases, for five of the phase-retrieval methods at three different noise levels. The object of Fig. 2(a) is propagated numerically and phase retrieval performed on the result, using the same parameters as in Fig. 6. The result is compared to the known object, generating Table 3. One column includes no noise, the others include Poisson noise at a pixel SNR of 20 and 4, respectively. For the two later cases, the table contains mean and standard deviation calculated from 1000 different realizations of the noisy image. Before calculation of the mean square error, the object has been normalized to a mean value of 1. Then the scale of the retrieved phase has been adjusted for minimal mean square error. The phase-attenuation duality is not included, as it applies to a different wavelength region, but would give results very similar to the single-material method. The two-material method is not included as it is a rather special case. The data in the table support the conclusions drawn from the images: except for the Bronnikov algorithm, all methods give similar results for high SNR. For low SNR, the modified Bronnikov and single-material methods give better results.

Table 3 Least Mean Square Error of Normalized Retrieved Phase Images for Different Noise Levels ¹

Method	No noise	SNR=20	SNR=4
Bronnikov	2.43	2.51 ± 0.11	3.30 ± 0.27
Modified Bronnikov	0.45	0.73 ± 0.06	2.53 ± 0.27
Single material	0.73	0.94 ± 0.09	2.64 ± 0.29
Fourier (Born)	0.45	2.01 ± 0.05	4.01 ± 0.07
Fourier (Rytov)	0.77	1.75 ± 0.05	3.99 ± 0.07

¹ All values should be multiplied by 1 × 10⁻³.

5. Discussion and conclusions

The different phase-retrieval methods suitable for in-line X-ray phase-contrast tomography have been characterized and compared, and most of the information concentrated into Tables 1 and 2. From the expressions in Table 1 and the levels of approximations in Table 2 we can construct a scheme for choosing the most suitable phase-retrieval method, as shown in Fig. 5. Finally, phase retrieval has been performed on simulated and experimental data, illustrating and confirming the scheme in Fig. 5. Comparison is done for the retrieved two-dimensional projections rather than for reconstructed three-dimensional slices, and the tomographic reconstruction is not considered in this paper. Optimally retrieved projections are likely to give improved reconstructed slices, although this has not been explicitly analyzed in this paper.

Figure 8(e)–8(f) contains the phase retrieved from far-field images. Methods suitable only to the near-field (in this case the single-material method) give a blurred reconstruction of the wrong dimensions, while the Fourier methods (in this case in the Rytov approximation) give a better reconstruction of the actual object, but also tend to amplify noise at particular frequencies. This can be seen in Figs. 7 and 9 where noise amplification is also present, from the expressions in Table 1 where the denominator of the spatial-frequency filter goes to zero for specific frequencies, and from the mean square errors in Table 3 which increase significantly in the presence of noise. So as Fig. 5 indicates, the Fourier methods are the only options in the far-field case, but for near-field images other methods give equivalent or better results.

If the photon energy is high and the object made from light materials, the phase-attenuation duality method is the best option as it is not limited to homogeneous objects. For lower beam energies it should not be used, as it gets the values of δ/β wrong as illustrated in Fig. 6. If the sample consists of two materials, there is some benefit to the two-material method, as illustrated in Fig. 8(b)–8(c).

The first phase-retrieval method to be developed, the Bronnikov method, is not included in the scheme in Fig. 5, and the reason for this is illustrated in Figs. 6, 7, and 9. While it gives reasonable results for noise-free data, any low-frequency noise is strongly amplified as the denominator of the spatial-frequency filter goes to zero for low frequencies. This is particularly obvious in experimental data where low-frequency noise is most likely present in the background. In place of the Bronnikov method, the modified Bronnikov method can be used. The methods are very similar, but the modified Bronnikov method avoids the strong amplification at low frequencies and thus gives better results.

The two most widely applicable methods are the modified Bronnikov and the single-material method. Though derived differently, they are very similar, and will give similar results.

In summary, the phase retrieval methods suitable for phase-contrast tomography have been characterized, and found to follow the same mathematical scheme. This information has been used to compare the methods under different circumstances, and to provide a simple strategy for choosing the most suitable phase-retrieval method.

Acknowledgments

Financial support from the Swedish Research Council and from the Swedish Foundation for Strategic Research is gratefully acknowledged.

References and links

1.	See, e.g., P. Suetens Fundamentals of Medical Imaging (Cambridge Univ Press, 2009). [CrossRef]
2.	R. Fitzgerald, “Phase-sensitive X-ray imaging,” Phys. Today 53(7), 23–26 (2000). [CrossRef]
3.	R. A. Lewis, “Medical phase contrast X-ray imaging: current status and future prospects,” Phys. Med. Biol. 49, 3573–3583 (2004). [CrossRef] [PubMed]
4.	U. Bonse and M. Hart, “An X-ray interferometer,” Appl. Phys. Lett. 6, 155–156 (1965). [CrossRef]
5.	T. J. Davis, D. Gao, T. E. Gureyev, A. W. Stevenson, and S. W. Wilkins, “Phase-contrast imaging of weakly absorbing materials using hard X-rays,” Nature 373, 595–598 (1995). [CrossRef]
6.	V. N. Ingal and E. A. Beliaevskaya, “X-ray plane-wave topographyobservation of the phase contrast from a non-crystalline object,” J. Phys. D: Appl. Phys. 28, 2314–2317 (1995). [CrossRef]
7.	A. Snigirev, I. Snigireva, V. Kohn, S. Kuznetsov, and I. Schelokov, “On the possibilities of X-ray phase contrast microimaging by coherent high-energy synchrotron radiation,” Rev. Sci. Instrum. 66, 5486–5492 (1995). [CrossRef]
8.	S. W. Wilkins, T. E. Gureyev, D. Gao, A. Pogany, and A. W. Stevenson, “Phase-constrast imaging using poly-chromatic hard X-rays,” Nature 384, 335–338 (1996). [CrossRef]
9.	K. A. Nugent, T. E. Gureyev, D. F. Cookson, D. Paganin, and Z. Barnea, “Quantitative phase imaging using hard X-rays,” Phys. Rev. Lett. 77, 2961–2964 (1996). [CrossRef] [PubMed]
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.	D. Chapman, W. Thomlinson, R. E. Johnston, D. Washburn, E. Pisano, N. Gmür, Z. Zhong, R. Menk, F. Arfelli, and D. Sayers, “Diffraction enhanced X-ray imaging,” Phys. Med. Biol. 42, 2015–2025 (1997). [CrossRef] [PubMed]
12.	K. A. Nugent, “Coherent methods in the X-ray sciences,” Adv. Phys. 59, 1–99 (2010). [CrossRef]
13.	F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nat. Phys. 2, 258–261 (2006). [CrossRef]
14.	A. Momose, “Recent advances in X-ray phase imaging,” Jpn. J. Appl. Phys. 44, 6355–6367 (2005). [CrossRef]
15.	T. Tuohimaa, M. Otendal, and H. M. Hertz, “Phase-contrast X-ray imaging with a liquid-metal-jet-anode micro-focus source,” App. Phys. Lett. 91, 074104 (2007). [CrossRef]
16.	R. A. Gonsalves, “Phase retrieval by differential intensity measurements,” J. Opt. Soc. Am. A 4, 166–170 (1987). [CrossRef]
17.	K. A. Nugent, “X-ray noninterferometric phase imaging: a unified picture,” J. Opt. Soc. Am. A 24, 536–546 (2007). [CrossRef]
18.	T. E. Gureyev, S. C. Mayo, D. E. Myers, Ya. Nesterets, D. M. Paganin, A. Pogany, A. W. Stevenson, and S. W. Wilkins, “Refracting Röntgen’s rays: propagation-based X-ray phase contrast for biomedical imaging,” J. Appl. Phys. 105, 102005 (2009). [CrossRef]
19.	A. Barty, K. A. Nugent, D. Paganin, and A. Roberts, “Quantitative optical phase microscopy,” Opt. Lett. 23, 817–819 (1998). [CrossRef]
20.	A. V. Bronnikov, “Reconstruction formulas for phase-contrast imaging,” Opt. Commun. 171, 239–244 (1999). [CrossRef]
21.	A. Groso, R. Abela, and M. Stampanoni, “Implementation of a fast method for high resolution phase contrast tomography,” Opt. Express 14, 8103–8110 (2006). [CrossRef] [PubMed]
22.	X. Wu and H. Liu, “X-Ray cone-beam phase tomography formulas based on phase-attenuation duality,” Opt. Express 13, 6000–6014 (2005). [CrossRef] [PubMed]
23.	D. Paganin, S. C. Mayo, T. E. Gureyev, P. R. Wilkins, and S. W. Wilkins, “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” J. Microsc. 206, 33–40 (2002). [CrossRef] [PubMed]
24.	M. A. Beltran, D. M. Paganin, K. Uesugi, and M. J. Kitchen, “2D and 3D X-ray phase retrieval of multi-material objects using a single defocus distance,” Opt. Express 18, 6423–6436 (2010). [CrossRef] [PubMed]
25.	T. E. Gureyev, T. J. Davis, A. Pogany, S. C. Mayo, and S. W. Wilkins, “Optical phase retrieval by use of first Born- and Rytov-type approximations,” Appl. Opt. 43, 2418–2430 (2004). [CrossRef] [PubMed]
26.	S. Zabler, P. Cloetens, J.-P. Guigay, and J. Baruchel, “Optimization of phase contrast imaging using hard x-rays,” Rev. Sci. Intstrum. 76, 073705 (2005). [CrossRef]
27.	J. P. Guigay, “Fourier transform analysis of Fresnel diffraction patterns and in-line holograms,” Optik 49, 121–125 (1977).
28.	L. D. Turner, B. B. Dhal, J. P. Hayes, A. P. Mancuso, K. A. Nugent, D. Paterson, R. E. Scholten, C. Q. Tran, and A. G. Peele, “X-ray phase imaging: Demonstration of extended conditions with homogeneous objects,” Opt. Express 12, 2960–2965 (2004). [CrossRef] [PubMed]
29.	M. R. Teague, “Deterministic phase retrieval: a Green’s function,” J. Opt. Soc. Am. 73, 1434–1441 (1983). [CrossRef]
30.	A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (Siam, 2001). [CrossRef]
31.	R. Grella, “Fresnel propagation and diffraction and paraxial wave equation,” J. Opt. (Paris) 13, 367–364 (1982).
32.	T. E. Gureyev, Y. I. Nesterets, A. W. Stevenson, P. R. Miller, A. Pogany, and S. W. Stevenson, “Some simple rules for contrast, signal-to-noise and resolution in in-line phase-contrast imaging,” Opt. Express 16, 3223–3241 (2008). [CrossRef] [PubMed]
33.	D. M. Paganin, Coherent X-Ray Optics (Oxford Science Publications, 2006). [CrossRef]
34.	J. M. Cowley, Diffraction Physics (Elsevier, 1995).
35.	M. J. Bastiaans, “Application of the Wigner distribution function to partially coherent light,” J. Opt. Soc. Am. A 3, 1227–1237 (1986). [CrossRef]
36.	X. Wu and H. Liu, “A new theory of phase-contrast X-ray imaging based on Wigner distributions,” Med. Phys. 31, 2378–2384 (2004). [CrossRef] [PubMed]
37.	G. R. Myers, S. C. Mayo, T. E. Gureyev, D. M. Paganin, and S. W. Wilkins, “Polychromatic cone-beam phase-contrast tomography,” Phys. Rev. A 76, 045804 (2007). [CrossRef]
38.	J. W. Goodman, Introduction to Fourier Optics , 2nd ed. (McGraw-Hill, 1996).
39.	U. Lundström, P. A. C. Takman, L. Scott, H. Brismar, and H. M. Hertz, “Low-dose high-resolution laboratory phase-contrast X-ray imaging,” manuscript in preparation.

OCIS Codes
(000.1430) General : Biology and medicine
(100.3190) Image processing : Inverse problems
(340.7440) X-ray optics : X-ray imaging

ToC Category:
X-ray Optics

History
Original Manuscript: February 22, 2011
Revised Manuscript: April 28, 2011
Manuscript Accepted: May 4, 2011
Published: May 11, 2011

Virtual Issues
Vol. 6, Iss. 6 Virtual Journal for Biomedical Optics

Citation
Anna Burvall, Ulf Lundström, Per A. C. Takman, Daniel H. Larsson, and Hans M. Hertz, "Phase retrieval in X-ray phase-contrast imaging suitable for tomography," Opt. Express 19, 10359-10376 (2011)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-19-11-10359

Sort: Author | Year | Journal | Reset

References

See, e.g., P. SuetensFundamentals of Medical Imaging (Cambridge Univ Press, 2009). [CrossRef]
R. Fitzgerald, “Phase-sensitive X-ray imaging,” Phys. Today 53(7), 23–26 (2000). [CrossRef]
R. A. Lewis, “Medical phase contrast X-ray imaging: current status and future prospects,” Phys. Med. Biol. 49, 3573–3583 (2004). [CrossRef] [PubMed]
U. Bonse and M. Hart, “An X-ray interferometer,” Appl. Phys. Lett. 6, 155–156 (1965). [CrossRef]
T. J. Davis, D. Gao, T. E. Gureyev, A. W. Stevenson, and S. W. Wilkins, “Phase-contrast imaging of weakly absorbing materials using hard X-rays,” Nature 373, 595–598 (1995). [CrossRef]
V. N. Ingal and E. A. Beliaevskaya, “X-ray plane-wave topographyobservation of the phase contrast from a non-crystalline object,” J. Phys. D: Appl. Phys. 28, 2314–2317 (1995). [CrossRef]
A. Snigirev, I. Snigireva, V. Kohn, S. Kuznetsov, and I. Schelokov, “On the possibilities of X-ray phase contrast microimaging by coherent high-energy synchrotron radiation,” Rev. Sci. Instrum. 66, 5486–5492 (1995). [CrossRef]
S. W. Wilkins, T. E. Gureyev, D. Gao, A. Pogany, and A. W. Stevenson, “Phase-constrast imaging using poly-chromatic hard X-rays,” Nature 384, 335–338 (1996). [CrossRef]
K. A. Nugent, T. E. Gureyev, D. F. Cookson, D. Paganin, and Z. Barnea, “Quantitative phase imaging using hard X-rays,” Phys. Rev. Lett. 77, 2961–2964 (1996). [CrossRef] [PubMed]
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]
D. Chapman, W. Thomlinson, R. E. Johnston, D. Washburn, E. Pisano, N. Gmür, Z. Zhong, R. Menk, F. Arfelli, and D. Sayers, “Diffraction enhanced X-ray imaging,” Phys. Med. Biol. 42, 2015–2025 (1997). [CrossRef] [PubMed]
K. A. Nugent, “Coherent methods in the X-ray sciences,” Adv. Phys. 59, 1–99 (2010). [CrossRef]
F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nat. Phys. 2, 258–261 (2006). [CrossRef]
A. Momose, “Recent advances in X-ray phase imaging,” Jpn. J. Appl. Phys. 44, 6355–6367 (2005). [CrossRef]
T. Tuohimaa, M. Otendal, and H. M. Hertz, “Phase-contrast X-ray imaging with a liquid-metal-jet-anode micro-focus source,” App. Phys. Lett. 91, 074104 (2007). [CrossRef]
R. A. Gonsalves, “Phase retrieval by differential intensity measurements,” J. Opt. Soc. Am. A 4, 166–170 (1987). [CrossRef]
K. A. Nugent, “X-ray noninterferometric phase imaging: a unified picture,” J. Opt. Soc. Am. A 24, 536–546 (2007). [CrossRef]
T. E. Gureyev, S. C. Mayo, D. E. Myers, Ya. Nesterets, D. M. Paganin, A. Pogany, A. W. Stevenson, and S. W. Wilkins, “Refracting Röntgen’s rays: propagation-based X-ray phase contrast for biomedical imaging,” J. Appl. Phys. 105, 102005 (2009). [CrossRef]
A. Barty, K. A. Nugent, D. Paganin, and A. Roberts, “Quantitative optical phase microscopy,” Opt. Lett. 23, 817–819 (1998). [CrossRef]
A. V. Bronnikov, “Reconstruction formulas for phase-contrast imaging,” Opt. Commun. 171, 239–244 (1999). [CrossRef]
A. Groso, R. Abela, and M. Stampanoni, “Implementation of a fast method for high resolution phase contrast tomography,” Opt. Express 14, 8103–8110 (2006). [CrossRef] [PubMed]
X. Wu and H. Liu, “X-Ray cone-beam phase tomography formulas based on phase-attenuation duality,” Opt. Express 13, 6000–6014 (2005). [CrossRef] [PubMed]
D. Paganin, S. C. Mayo, T. E. Gureyev, P. R. Wilkins, and S. W. Wilkins, “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” J. Microsc. 206, 33–40 (2002). [CrossRef] [PubMed]
M. A. Beltran, D. M. Paganin, K. Uesugi, and M. J. Kitchen, “2D and 3D X-ray phase retrieval of multi-material objects using a single defocus distance,” Opt. Express 18, 6423–6436 (2010). [CrossRef] [PubMed]
T. E. Gureyev, T. J. Davis, A. Pogany, S. C. Mayo, and S. W. Wilkins, “Optical phase retrieval by use of first Born- and Rytov-type approximations,” Appl. Opt. 43, 2418–2430 (2004). [CrossRef] [PubMed]
S. Zabler, P. Cloetens, J.-P. Guigay, and J. Baruchel, “Optimization of phase contrast imaging using hard x-rays,” Rev. Sci. Intstrum. 76, 073705 (2005). [CrossRef]
J. P. Guigay, “Fourier transform analysis of Fresnel diffraction patterns and in-line holograms,” Optik 49, 121–125 (1977).
L. D. Turner, B. B. Dhal, J. P. Hayes, A. P. Mancuso, K. A. Nugent, D. Paterson, R. E. Scholten, C. Q. Tran, and A. G. Peele, “X-ray phase imaging: Demonstration of extended conditions with homogeneous objects,” Opt. Express 12, 2960–2965 (2004). [CrossRef] [PubMed]
M. R. Teague, “Deterministic phase retrieval: a Green’s function,” J. Opt. Soc. Am. 73, 1434–1441 (1983). [CrossRef]
A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (Siam, 2001). [CrossRef]
R. Grella, “Fresnel propagation and diffraction and paraxial wave equation,” J. Opt. (Paris) 13, 367–364 (1982).
T. E. Gureyev, Y. I. Nesterets, A. W. Stevenson, P. R. Miller, A. Pogany, and S. W. Stevenson, “Some simple rules for contrast, signal-to-noise and resolution in in-line phase-contrast imaging,” Opt. Express 16, 3223–3241 (2008). [CrossRef] [PubMed]
D. M. Paganin, Coherent X-Ray Optics (Oxford Science Publications, 2006). [CrossRef]
J. M. Cowley, Diffraction Physics (Elsevier, 1995).
M. J. Bastiaans, “Application of the Wigner distribution function to partially coherent light,” J. Opt. Soc. Am. A 3, 1227–1237 (1986). [CrossRef]
X. Wu and H. Liu, “A new theory of phase-contrast X-ray imaging based on Wigner distributions,” Med. Phys. 31, 2378–2384 (2004). [CrossRef] [PubMed]
G. R. Myers, S. C. Mayo, T. E. Gureyev, D. M. Paganin, and S. W. Wilkins, “Polychromatic cone-beam phase-contrast tomography,” Phys. Rev. A 76, 045804 (2007). [CrossRef]
J. W. Goodman, Introduction to Fourier Optics , 2nd ed. (McGraw-Hill, 1996).
U. Lundström, P. A. C. Takman, L. Scott, H. Brismar, and H. M. Hertz, “Low-dose high-resolution laboratory phase-contrast X-ray imaging,” manuscript in preparation.

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper