## Abstract

For coherent X-ray imaging of pure phase objects we study the reliability of linear relations in phase-retrieval algorithms based on a single intensity map after free-space propagation. For large phase changes and/or large propagation distances we propose two venues of working beyond linearity: Projection onto an *effective*, linear and *local* model in Fourier space and expansion of intensity contrast in powers of object-detector distance. We apply both algorithms successfully to simulated data.

© 2011 OSA

## 1. Introduction

Third-generation synchrotron light sources produce brilliant photon beams of spatial and temporal coherence properties at the sample which are suitable for routine application of phase-sensitive X-ray imaging methods. Within low-*Z* samples of low atomic number density (polymers, soft biological tissue, etc.), the attenuation of highly energetic X-rays is weak, thus yielding poor absorption contrast. However, such samples introduce sizable phase shifts to X-ray wave fronts creating intensity contrast [1] downstream of the object due to a free-space propagated transmission function [2–5]. Notice that phase contrast can be 10^{3} times larger than absorption contrast [6].

The present work reports on two novel and complementary theoretical approaches addressing the nonlinear, noniterative, and single-distance phase-retrieval problem. Their usefulness is demonstrated by numerical experiments: For pure-phase objects, parallel-beam geometry [7], and in Fresnel diffraction theory we consider the regime where the approximation of a linear and local relation between the phase shift *ϕ _{z}*

_{=0}in the object-exit plane and intensity

*I*at a finite object-detector distance

_{z}*z*is inapplicable. This happens for large relative phase shifts and/or large

*z*values (former relevant to phase-sensitive imaging, e.g., of biological organisms, latter occurring, e.g., in measurements with bulky sample environments). That is, the relation between

*ϕ*

_{z}_{=0}and

*I*can be nonlinear and nonlocal because relative phase shifts in the object plane or the value of

_{z}*z*, where intensity is measured, or both are

*large*.

## 2. Review of linear approaches

Let us briefly review two widely used linear relations between *ϕ _{z}*

_{=0}and

*I*pointing out their strengths and limitations. In the

_{z}*small-z limit*the transport-of-intensity equation (

*TIE*) (imaginary part of the paraxial equation [5, 8]) leads to the approximation

*r⃗*a vector in the plane transverse to the optical axis, ∇

_{⊥}the associated planar nabla operator, $k\hspace{0.17em}\equiv \hspace{0.17em}\frac{2\pi E}{\mathit{\text{hc}}}$ denotes the wave number of the incident monochromatic X-rays, and

*h*,

*c*are Planck’s quantum of action, speed of light in vacuum, respectively. Moreover,

*E*=

*hν*is the energy of a photon associated with a wave of frequency

*ν*. In the following we refer to phase retrieval based on Eq. (1) as

*linearized TIE*because of the truncation at linear order in

*z*of the Taylor expansion of

*I*(

*z*) in the TIE. (This allows to express

*∂*(

_{z}I*z*) at

*z*= 0 linearly in terms of

*g*(

*z*) yielding Eq. (1).) Equation (1) states that out of a given phase map behind the object intensity contrast strengthens with increasing distance

*z*. In the

*small-relative-phase limit*,

*ξ⃗*is the Fourier conjugate to

*r⃗*, an important nonlinear and nonlocal relation between the autocorrelation of the object’s transmission function

*ψ*

_{z}_{=0}(

*r⃗*) and the 2D Fourier transform (

*ℱI*)(

_{z}*ξ⃗*) [9]

*contrast-transfer function*[10–14]. Notice that relation (3) is tied to Fresnel diffraction theory, the Fraunhofer limit of the paraxial approximation would yield a different relation between

*ψ*

_{z}_{=0}(

*r⃗*) and (

*ℱ I*) (

_{z}*ξ⃗*). Notice also that condition (2) is required to facilitate the truncation at linear order of the exponential of $\text{i}\hspace{0.17em}\left({\varphi}_{z=0}\hspace{0.17em}\left(\overrightarrow{r}\hspace{0.17em}-\hspace{0.17em}\frac{\pi z}{k}\overrightarrow{\xi}\right)\hspace{0.17em}-\hspace{0.17em}{\varphi}_{z=0}\hspace{0.17em}\left(\overrightarrow{r}\hspace{0.17em}+\hspace{0.17em}\frac{\pi z}{k}\overrightarrow{\xi}\right)\right)$. This exponential arises upon substituting ${\psi}_{z=0}\hspace{0.17em}=\hspace{0.17em}\sqrt{{I}_{z=0}}{e}^{\text{i}{\varphi}_{z=0}}$ into Eq. (3).

The position-space version of the contrast-transfer function reads

*r⃗*and

*ξ⃗*to guarantee the validity of Eq. (4).

In Fourier space relation (4) reads

*CTF*. Truncating the power-series expansion of the sine at linear order, yields the Fourier-space version of Eq. (1). Thus Eq. (4) reduces to Eq. (1) in the limit of small values of $\frac{z}{k}{\overrightarrow{\xi}}_{\text{max}}^{2}$ where ${\left|\overrightarrow{\xi}\right|}_{\text{max}}\hspace{0.17em}=\hspace{0.17em}\frac{1}{2\Delta x}$ is the upper cutoff for transverse spatial frequency set by the effective linear pixel size Δ

*x*of the detector.

Figures 1(b), 1(c), and 1(d) depict function log |*ℱ g _{z}*|(

*ξ⃗*) based on the standard Lena test pattern

*ϕ*

_{z}_{=0}(

*r⃗*) of Fig. 1(a), scaled to a maximal phase shift

*ϕ*

_{max}= 0.01, 1, 6, respectively. The frame-like structure in log |

*ℱ g*|(

_{z}*ξ⃗*) is due to the Gaussian-blurred transition from the zero-padded region to the actual phase map, see Fig. 1(a), and truncation rods occur (crosses centered at origin). Notice the decrease in visibility of the sinusoidal, radial modulation (rings of minima) from Figs. 1(b) to 1(d). This indicates the failure of CTF at large

*ϕ*

_{max}. Phase retrieval then is affected by quasiperiodic texture superimposed onto the exact phase map, see Figs. 2(b), 2(c): While Eq. (5) expects a sinusoidal modulation of (

*ℱ g*)(

_{z}*ξ⃗*) this modulation does

*not*occur at large

*ϕ*

_{max}. As a result, the zeros of the sine in Eq. (5)

*artificially*enhance the spectrum of

*ϕ*

_{z}_{=0}by simple poles in

*ξ⃗*

^{2}. Since |

*ℱ g*|(

_{z}*ξ⃗*) decays towards large |

*ξ⃗*| the dominating effect comes from the singularities associated with the most central ring of minima. Upon inverse Fourier transformation this yields in transverse position space a 2D quasiperiodic artifact of wavelength $\sqrt{\frac{2\pi z}{k}}$. The amplitude of these quasioscillations rapidly grows with increasing

*ϕ*

_{max}. In contrast to that, quasiperiodic artifacts are absent in linearized TIE, compare Figs. 2(c) and 2(d). For small relative phase shifts Fig. 2(a) demonstrates the superiority of CTF over linearized TIE. When increasing

*ϕ*

_{max}this is no longer true, see Figs. 2(b) and 2(c).

Condition (2) can also be satisfied for large *ϕ*
_{max} if
$z\hspace{0.17em}\ll \hspace{0.17em}\frac{k\Delta x}{\sqrt{2}\pi \hspace{0.17em}|{\nabla}_{\perp}\hspace{0.17em}{\varphi}_{z=0}{|}_{\text{max}}}$ where the maximum of |∇_{⊥}
*ϕ _{z}*

_{=0}| is taken over the entire field of view. Thus Eq. (4) is valid at large phase shifts when using sufficiently small values of

*z*. In practice this may or may not be an option considering constraints imposed by a limited exposure time (imaging of processes, dose restriction, detector dynamics) and/or a bulky sample environment. Small exposure times imply large statistical noise beating the weakly developed intensity contrast at small

*z*. To investigate this we introduce Poisson noise to the computed intensity at each pixel before phase retrieval. We use a mean value of

*N*times the computed, noise-free intensity ratio

_{c}*I*/

_{z}*I*

_{z}_{=0}where

*N*is the expected number of photons detected by a pixel. At

_{c}*I*/

_{z}*I*

_{z}_{=0}∼ 1 and

*N*= 8000 this yields ∼ 1% noise which we employ from now on. In Figs. 3(a) and 3(b) results of CTF for two distinct distances and energies are shown. Compared to Fig. 3(b) no large texture artifact occurs in Fig. 3(a) because CTF acts like a linearized TIE retrieval: $\text{sin}\hspace{0.17em}\left(2\frac{{\pi}^{2}z}{k}{\overrightarrow{\xi}}^{2}\right)\hspace{0.17em}\sim \hspace{0.17em}2\frac{{\pi}^{2}x}{k}{\overrightarrow{\xi}}^{2}$ for all

_{c}*ξ⃗*. Notice, however, the importance of noise artifacts due to the small propagation distance. Figure 3(b) depicts CTF at

*z*= 50cm and

*E*= 10keV where noise effects are marginal but a large texture artifact occurs. On the other hand, linearized TIE at

*z*= 50cm and

*E*= 10keV yields reasonable results but suffers from a limited resolution, compare Fig. 3(c) with Fig. 3(a).

## 3. Nonlinear, noniterative approaches

The above discussion indicates the need for nonlinear approaches to the retrieval of large relative phases.

#### 3.1. Quasiparticle approach

One way to work beyond *fundamental* linearity is to appeal to an *effective* linear and local model in Fourier space. Namely, we propose a modified CTF retrieval, dubbed *projected CTF*, in the following sense: Artificial peaks in the spectrum of *ϕ _{z}*

_{=0}introduced by CTF are removed by considering an intensity contrast which is filtered in Fourier space, $(\tilde{\mathcal{F}\hspace{0.17em}{g}_{z}})\hspace{0.17em}(\overrightarrow{\xi})\hspace{0.17em}\equiv \hspace{0.17em}\Theta \hspace{0.17em}\left(\left|\text{sin}\hspace{0.17em}\left(\frac{2{\pi}^{2}z}{k}{\overrightarrow{\xi}}^{2}\right)\right|\hspace{0.17em}-\hspace{0.17em}\varepsilon \right)\hspace{0.17em}\cdot \hspace{0.17em}(\mathcal{F}\hspace{0.17em}{g}_{z})(\overrightarrow{\xi})$. Here $\frac{2{\pi}^{2}z}{k}{\overrightarrow{\xi}}^{2}\hspace{0.17em}>\hspace{0.17em}\frac{\pi}{2}$ Θ denotes the Heaviside step function,

*ε*is the threshold for this binary filter (0 ≤

*ε*< 1), and on the left-hand side of Eq. (5) (

*ℱ g*)(

_{z}*ξ⃗*) is replaced by $(\tilde{\mathcal{F}\hspace{0.17em}{g}_{z}})\hspace{0.17em}(\overrightarrow{\xi})$. Let us discuss two physical aspects of projected CTF. First, to assure the absence of singularities in the effective phase

*ℱ ϕ*

_{z}_{=0}the above binary filtering is used. Second, the use of a

*local*relation between $(\tilde{\mathcal{F}\hspace{0.17em}{g}_{z}})\hspace{0.17em}(\overrightarrow{\xi})$ and (

*ℱ ϕ*

_{z}_{=0})(

*ξ⃗*) renders the latter an effective quantity. Namely, nonlinear effects in the forward propagation of the exact phase fundamentally imply a

*nonlocal*relation between the Fourier transform of the latter and

*ℱ g*. Exemplarily, this can be checked by truncating at quadratic order the expansion of the exponential of phase after substitution of ${\psi}_{z=0}\hspace{0.17em}=\hspace{0.17em}\sqrt{{I}_{z=0}}\hspace{0.17em}\text{exp}\hspace{0.17em}(\text{i}{\varphi}_{z=0})$ into Eq. (3). The resulting relation between

_{z}*ℱ I*and

_{z}*ℱ ϕ*

_{z}_{=0}is (by construction) nonlinear and (by the Fourier convolution theorem) nonlocal. As a result, a

*ξ⃗*dependent rescaling between the Fourier transform of the exact phase and

*ℱ ϕ*

_{z}_{=0}as retrieved from our modified, local version of Eq. (5) occurs. But for phase retrieval in transverse position space, an

*inverse*Fourier transformation needs to be applied which

*averages out*spurious spectral information contained in

*ℱ ϕ*

_{z}_{=0}. According to our simulations, the result indeed is close to the exact phase. An analogy is the successful quasiparticle concept in the quantum theory of condensed-matter physics and quantum field theory where a free dispersion law is altered by strong interactions yielding a free quasiparticle. In our case, the concept of a quasiparticle translates into the linear and local relation between effective phase and filtered intensity while the average over particle momenta in some normalized derivative of the partition function is analogous to inverse Fourier transformation. Namely, an

*r⃗*dependent ‘action’

*S*≡ 2

_{r⃗}*πξ⃗*·

*r⃗*is used to define the partition function

*Z*as

_{r⃗}*Z*= tr exp(i

_{r⃗}*S*) where the trace symbol is understood as a sum over all wave-vector states. We stress that projected CTF is nonperturbative. Figure 4(a) depicts $\text{log}\hspace{0.17em}\left|\tilde{\mathcal{F}\hspace{0.17em}{g}_{z}}\right|$ at

_{r⃗}*ε*= 0.01. Comparison of Figs. 4(b) and 4(c) reveals that an increase of resolution is traded for an increase of the texture artifact when decreasing the threshold

*ε*. Projected CTF essentially improves CTF, compare Fig. 3(b) with Figs. 4(b), 4(c). Notice that in the limit

*ε*→ 1 projected CTF yields results that are similar to those of linearized TIE, compare Fig. 4(c) with Fig. 3(c).

#### 3.2. Perturbation theory

Let us now consider another nonlinear phase-retrieval algorithm which, in contrast to CTF, is perturbative. Namely, expanding the intensity contrast as
${g}_{z}\hspace{0.17em}(\overrightarrow{r})\hspace{0.17em}\equiv \hspace{0.17em}{\sum}_{i=1}^{N}{g}^{(i)}\hspace{0.17em}(\overrightarrow{r}){z}^{i}$, the idea is to determine coefficient *g*
^{(1)}(*r⃗*) as precisely as possible to solve
${g}^{(1)}\hspace{0.17em}(\overrightarrow{r})\hspace{0.17em}=\hspace{0.17em}-\hspace{0.17em}\frac{1}{2k}\hspace{0.17em}{\nabla}_{\perp}^{2}\hspace{0.17em}{\varphi}_{z=0}\hspace{0.17em}(\overrightarrow{r})$, see Eq. (1). In [15] *g*
^{(1)}(*r⃗*) is extracted from a fit to a stack of intensity measurements {*Iz _{i}*|

*i*= 1, ⋯ ,

*N*} the better the larger

*N*is. In [16] a theoretical alternative to that approach, which only uses a single-distance intensity map, was proposed: All coefficients

*g*

^{(}

^{i}^{)}(

*r⃗*) are, by virtue of the full paraxial equation, expressed by 2

*i*transverse derivatives, appearing in terms that are linear and nonlinear in

*ϕ*

_{z}_{=0}. Truncating the expansion at

*N*≥ 2, one obtains a nonlinear partial differential equation (PDE) for

*ϕ*

_{z}_{=0}subject to the single-distance source term

*g*(

_{z}*r⃗*) [16]. At

*N*= 2 or next-to-leading order (

*NLO*) the explicit expansion of

*g*(

_{z}*r⃗*) reads [16]

*ϕ*

_{z}_{=0}to the nonlinear PDE (6) perturbatively is to solve it at leading order

*N*= 1, to use this solution in estimating the next-to-leading order (we refer to the result as

*PNLO*), and to finally invert the Laplacian on $-\frac{k}{z}({g}_{z}(\overrightarrow{r})\hspace{0.17em}-\hspace{0.17em}\text{PNLO})$[16]. This type of nonlinear phase retrieval is dubbed

*linearized TIE + PNLO*. As Fig. 5 indicates linearized TIE + PNLO is useful at very large relative phase shifts. Obviously, CTF yields unacceptable results. Comparing Figs. 5(d) and 5(c), linearized TIE + PNLO improves linearized TIE within regions where the phase varies sufficiently slowly, compare also with [16]. This lowers the mean retrieval error per pixel by about 40%. Comparing linearized TIE + PNLO with projected CTF, see Fig. 5(d) and Fig. 5(b), the mean retrieval error per pixel in the former is reduced by about 42% w.r.t. the latter. Linearized TIE + PNLO, however, may artificially hollow regions which are bounded by strong phase slopes (transition across Lena’s eyes). Finally, let us remark that PNLO can be modified to

*PNLO*′ by using CTF to estimate the nonlinear terms in Eq. (6) [17]. In simulations we see that CTF + PNLO′ works whenever CTF already is excellent (typically for

*ϕ*

_{max}< 0.01). Practically, this is irrelevant because of the smallness of the correction PNLO′. For

*ϕ*

_{max}≥ 0.1 PNLO′ exhibits large texture artifacts.

## 4. Summary

To summarize, we have proposed two nonlinear, noniterative single-distance phase-retrieval algorithms. In analogy to the quasiparticle approach in strongly coupled systems [18], an effective phase in Fourier space, satisfying certain minimal requirements, assures at moderately large phase shifts an essential improvement over fundamentally linear methods. Alternatively, we have considered a perturbative approach [16] to the power-series expansion of intensity contrast [15] which works well for very large phase shifts. Our results could be useful for *static* phase-contrast imaging, but they also are applicable to the imaging of structure *dynamics* in nonabsorbing materials for which single-distance retrieval is essential. This extends to the case of weak and quasistatic absorption: An absorption image *I _{z}*

_{=0}would be taken to define intensity contrasts {

*g*} for a sequence of images recorded at a single distance

_{z}*z*> 0, think of the quick motion of nonabsorbing air bubbles inside a quasistatic, absorbing frame containing a nonabsorbing fluid. Another way of including absorption effects in single-distance phase retrieval is to make the assumption that

*ϕ*

_{z}_{=0}=

*αB*

_{z}_{=0}in ${\psi}_{z=0}\hspace{0.17em}=\hspace{0.17em}\sqrt{{I}_{z=0}}\text{exp}(\text{i}{\varphi}_{z=0}\hspace{0.17em}-\hspace{0.17em}{B}_{z=0})$, compare with Eq. (3). Here

*α*is a real constant that is known a priori, see for example [19].

Let us briefly compare our noniterative, single-distance, nonlinear phase-retrieval algorithms with iterative nonlinear methods, see [20]. First, iterative methods employ a sequence of forward-backward propagations. Backward propagation is simply implemented in case of the full Helmholtz theory (elementary wave fronts are spherical) and the Fraunhofer far-field approximation of the paraxial theory (elementary wave fronts are scaled 2D planes). In Fresnel theory the backward-propagation part is problematic since the theory is not symmetric under *z* → –*z* (elementary wave fronts are paraboloids). We insist on Fresnel theory because in Fraunhofer theory the Fourier transform of intensity at *z*, *ℱ I _{z}*, is related to the integral (and not the Fourier transform) of the autocorrelation of the object function

*ψ*

_{0}over the object plane. This implies that, unlike in Fresnel theory, the small-relative-phase condition (2) does not yield a linear relation between

*ℱ I*and

_{z}*ℱ ϕ*

_{z}_{=0}. Namely, the lowest-order relation between

*ϕ*

_{z}_{=0}and

*ℱ I*is

_{z}The present work should be relevant for soft-matter and life sciences. In fact, preliminary experimental results on the propagation phase contrast induced by Xenopus (African Clawed Frog) embryos show that CTF is not applicable and indicate that projected CTF indeed yields a substantial improvement over linearized TIE [21] in terms of spatial resolution.

## Acknowledgments

The authors acknowledge useful conversations with V. Altapova, A. Bronnikov, D. Hänschke, H. Suhonen, and L. Waller.

## References and links

**1. **B. Henke, E. Gullikson, and J. Davis, “X-ray interactions: photoabsorption, scattering, transmission, and reflection at E=50-30000 eV, z=1–92,” At. Data. Nucl. Data Tables **54**, 181–342 (1993). [CrossRef]

**2. **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]

**3. **P. Cloetens, “Contribution to phase contrast imaging, reconstruction and tomography with hard synchrotron radiation,” PhD dissertation, Vrije Universiteit Brussel (1999).

**4. **S. W. Wilkins, T. E. Gureyev, D. Gao, A. Pogany, and A. W. Stevenson, “Phase-contrast imaging using polychromatic hard X-rays,” Nature **384**, 335–338 (1996). [CrossRef]

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

**6. **A. Momose and J. Fukuda, “Phase-contrast radiographs of nonstained rat cerebellar specimen,” Med. Phys. **22**, 375–379 (1995). [CrossRef] [PubMed]

**7. **The cone-beam case relates to the parallel-beam case via a simple rescaling operation, see D. M. Paganin, *Coherent X-Ray Optics* (Oxford University Press, 2006). [CrossRef]

**8. **M. R. Teague, “Deterministic phase retrieval: a Greens function solution,” J. Opt. Soc. Am. **73**, 1434–1441 (1983). [CrossRef]

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

**10. **P. Cloetens, W. Ludwig, J. Baruchel, J.-P. Guigay, P. Rejmankova-Pernot, M. Salome, M. Schlenker, J. Y. Buffiere, E. Maire, and G. Peix, “Hard X-ray phase imaging using simple propagation of a coherent synchrotron radiation beam,” J. Phys. D: Appl. Phys. **32**, A145–A151 (1999). [CrossRef]

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

**12. **T. E. Gureyev, D. M. Paganin, G. R. Myers, Y. I. Nesterets, and S. W. Wilkins, “Phase-and-amplitude computer tomography,” Appl. Phys. Lett. **89**, 034102 (2006). [CrossRef]

**13. **T. E. Gureyev, Y. Nesterets, D. Paganin, A. Pogany, and S. Wilkins, “Linear algorithms for phase retrieval in the Fresnel region. 2. partially coherent illumination,” Opt. Commun. **259**, 569–580 (2006). [CrossRef]

**14. **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**, 073705 (2005). [CrossRef]

**15. **L. Waller, L. Tian, and G. Barbastathis, “Transport of intensity phase-amplitude imaging with higher order intensity derivatives,” Opt. Express **18**, 12552–12561 (2010). [CrossRef] [PubMed]

**16. **J. Moosmann, R. Hofmann, A. V. Bronnikov, and T. Baumbach, “Nonlinear phase retrieval from single-distance radiograph,” Opt. Express **18**, 25771–25785 (2010). [CrossRef] [PubMed]

**17. **J. Moosmann, R. Hofmann, and T. Baumbach, “Nonlinear, single-distance phase retrieval and Schwinger regularization,” To appear in Phys. Status. Solidi A.

**18. **L. D. Landau, “The theory of a Fermi liquid,” Sov. Phys. JETP **3**, 920–925 (1957).

**19. **M. Langer, P. Cloetens, and F. Peyrin, “Regularization of phase retrieval with phase-attenuation duality prior for 3-D holotomography,” IEEE Transactions on Image Processing **19**, 2429–2436 (2010). [CrossRef]

**20. **J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. **21**, 2758–2769 (1982). [CrossRef] [PubMed]

**21. **V. Altapova, D. Hänschke, J. Moosmann, R. Hofmann, and T. Baumbach, “Imaging in evolutionary biology: Single-distance phase contrast, nonlinear-noniterative retrieval, and tomographic reconstruction,” In preparation (2011).