Three-dimensional propagation in near-field tomographic X-ray phase retrieval

Ruhlandt, A.; Salditt, T.

doi:10.1107/S2053273315022469

research papers

FOUNDATIONS
ADVANCES

ISSN: 2053-2733

Volume 72| Part 2| March 2016| Pages 215-221

https://doi.org/10.1107/S2053273315022469

Open

access

Three-dimensional propagation in near-field tomographic X-ray phase retrieval

Aike Ruhlandt ^a ^* and Tim Salditt ^a

^aInstitut für Röntgenphysik, Georg-August-Universität Göttingen, Friedrich-Hund-Platz 1, Göttingen, Germany
^*Correspondence e-mail: aruhlan@gwdg.de

Edited by J. Miao, University of California, Los Angeles, USA (Received 1 June 2015; accepted 24 November 2015; online 29 January 2016)

This paper presents an extension of phase retrieval algorithms for near-field X-ray (propagation) imaging to three dimensions, enhancing the quality of the reconstruction by exploiting previously unused three-dimensional consistency constraints. The approach is based on a novel three-dimensional propagator and is derived for the case of optically weak objects. It can be easily implemented in current phase retrieval architectures, is computationally efficient and reduces the need for restrictive prior assumptions, resulting in superior reconstruction quality.

Keywords: tomography; phase retrieval; X-ray imaging.

1. Introduction

Over the last two decades the capabilities of X-ray tomography have been significantly improved by phase contrast methods. As opposed to conventional X-ray tomography, which is sensitive only to the amplitude (transmission) of the X-ray wave traversing an object, phase contrast techniques also take the phase shifting properties of the objects into account. This enables the visualization of weakly or non-absorbing soft tissues in biomedical imaging, or nanoscale structures in material science. Phase contrast is easily achieved without additional optical elements by free space propagation. The associated self interference of the object's exit wave over the controllable distance Δ between object and detector (Wilkins et al., 1996 ; Cloetens et al., 1999 ; Paganin & Nugent, 1998 ) encodes phase information into measurable intensities. The standard approach in evaluating phase contrast tomography data is a two-step reconstruction procedure. First, the phase retrieval is carried out, i.e. the complex valued exit wave is retrieved from the intensity measurements, separately for each projection angle. Secondly, all projections are combined to a three-dimensional volume using inverse Radon transform methods, in most cases a filtered back-projection (FBP). This is in sharp contrast to far-field coherent diffractive imaging, where phase retrieval is performed not only on projections, but commonly directly in three dimensions (Miao et al., 2001 , 2005 ; Chapman et al., 2006 ).

The phase retrieval step is considered to be the main challenge. Considerable efforts have concentrated on phase contrast algorithms beyond the simple but flawed holographic reconstruction by numerical back-propagation (Paganin, 2006 ). Deterministic but approximative solutions have been formulated based on the transport-of-intensity equation (TIE) (Paganin et al., 2002 ; Groso et al., 2006 ; Bronnikov, 1999 ), or on the analytic form of the free space contrast transfer functions (CTF) (Cloetens et al., 1999; Turner et al., 2004 ; Gureyev et al., 2004 ; Langer et al., 2012 ; Hofmann et al., 2011 ). The general strategy to overcome reconstruction artifacts such as the well known twin-image problem is to either use data sets with more than a single object-to-detector distance (Cloetens et al., 1999; Allen & Oxley, 2001 ; Latychevskaia & Fink, 2007 ; Hong et al., 2012 ; Krenkel et al., 2013 ), restrictive prior information such as known relationships between phase shift and absorption (Paganin, 2006; Wu et al., 2005 ) or known compact support of the object (Gerchberg & Saxton, 1972 ; Giewekemeyer et al., 2011 ; Bartels et al., 2012 ). The constraints due to prior information are usually implemented by iterative projection algorithms.

To overcome these limitations and provide a solution in particular for tomography scans at a single detection distance Δ, which is the most relevant case in practice, we and others have recently proposed a coupling of the two steps, phase retrieval and tomographic reconstruction (Kostenko et al., 2013 ; Ruhlandt et al., 2014 ). We could show that a combined phase retrieval and algebraic tomographic reconstruction scheme termed `iterative reprojection phase retrieval' (IRP) enhances the reconstruction quality and allows one to retrieve the phase of objects of mixed composition without the need for additional a priori knowledge (Ruhlandt et al., 2014). In particular, IRP was found to stabilize the reconstruction of low spatial frequencies which have previously hampered single-distance phase retrieval. This differs from combinations of the phase retrieval and FBP published before, which were implemented not to achieve coupling, but enhanced speed (Gureyev et al., 2004; Bronnikov, 2002 ). We attribute the enhancements in quality to the tomographic consistency condition (Helgason, 1965 ; Ludwig, 1966 ), which states that the different projections of an object are not independent from each other. This has been found very useful before in conventional absorption tomography, for example in the reconstruction of incomplete data (missing wedge, Kudo & Saito, 1991 ). Unfortunately, the gain in performance by IRP came at the cost of substantial computational effort, as well known also from previous algebraic tomographic reconstruction (ART) for the case of conventional absorption CT (computed tomography) (Gordon et al., 1970 ; Kak & Slaney, 1988 ). Hence, the combined phase retrieval and ART approach could so far be applied only to small data volumes.

In this work, we show that phase retrieval for optically weak objects can be enhanced by an inversion of the two steps, i.e. by first performing the inverse Radon transform, followed by a computationally efficient three-dimensional phase retrieval. The first step can be implemented for example by a fast FBP, the second by three-dimensional fast Fourier transformations (FFT). This approach has been used before in the context of deterministic one-step phase retrieval (Frank & Penczec, 1995 ; Cloetens et al., 1997 ). In this work we show that in combination with iterative phase retrieval, this three-dimensional propagation scheme exhibits superior reconstruction quality. Furthermore, propagation of entire three-dimensional objects sheds new light on the nature of the three-dimensional phase contrast tomography problem. In this article, first the theoretical concept is introduced, followed by a numerical implementation demonstrating the capabilities of the approach.

2. Theory

In X-ray imaging, the interaction of the radiation with an object is usually described by a three-dimensional index of refraction,

$[n(x,y,z) = 1-\delta(x,y,z)+i\beta(x,y,z), \eqno (1)]$

where $[1 - \delta]$ causes the phase shift and the imaginary part $[i\beta]$ the absorption in the object for a given wavelength λ. Notice that we consider monochromatic radiation. The illumination wave is treated as a plane wave $[\Psi = \Psi_0 \exp(i k_0 z]$ ) propagating along the optical axis z with the wavenumber $[k_0 = 2\pi/\lambda]$ . For a homogeneous slab of thickness t the phase shift induced by the object with respect to the propagation in free space is simply given by $[-k_0 \delta t]$ , and the wave amplitude by $[\exp (- k_0 \beta t)]$ . In most practical applications of hard X-ray coherent imaging, the object is weak enough that the projection approximation holds (Thibault, 2007 $[Thibault, P. (2007). Algorithmic Methods in Diffraction Microscopy. Ithaca: Cornell University.]$ ) and the propagation of the wave within the sample can be neglected. The distribution of the index n(x,y,z) then leads to a spatially modulated exit wave $[\Psi_{\rm e}(x,y,z = 0)]$ in the plane z = 0 directly behind the object, which is determined only by the projection $[{\cal R}(n-1): = \int (n-1) \, {\rm d}z]$ of the index of refraction along the optical axis:

$[\Psi_{\rm e}(x,y,z = 0) = \Psi_0(x,y) \exp[i k_0 {\cal R} (n-1)]. \eqno (2)]$

Further, for optically weak objects $[|k_0 {\cal R} (n-1)| \ll 1]$ the exponential function can be linearized to

$[\Psi_{\rm e}(x,y,z = 0) \simeq \Psi_0 \left[1+ i k_0 {\cal R} (n-1) \right]. \eqno (3)]$

The propagation $[{\cal D}(\Delta)]$ of the exit wave $[\Psi_{\rm e}]$ in free space can be expressed by a multiplication in Fourier space. For a given propagation distance Δ the two-dimensional (2d) Fourier transform $[{\cal F}_{2{\rm d}} \Psi_{\rm e}]$ is multiplied point-wise with the radially symmetric propagator $[P(k_x,k_y,\Delta): =]$ $[\exp[i \Delta (k_0^2-k_x^2 - k_y^2)^{1/2}]]$ (Fourier space coordinates k_x, k_y) followed by a Fourier back-transform $[{\cal F}^{-1}_{2{\rm d}}]$ (Paganin, 2006):

$[\eqalignno { \Psi_\Delta(x,y,z = \Delta)& = {\cal D}(\Delta) \Psi_{\rm e}(x,y,z = 0) &\cr & = {\cal F}^{-1}_{2{\rm d}}\left [P(\Delta) {\cal F}_{2{\rm d}} \Psi_{\rm e}(x,y,z = 0) \right] &\cr & \simeq \Psi_0 \exp(ik_0 \Delta) + \Psi_0 i k_0 {\cal F}^{-1}_{2{\rm d}} [P(\Delta)&\cr &\quad \times {\cal F}_{2{\rm d}} {\cal R} (n-1) ]. &(4)}]$

In the last step, the weak object approximation has been used. According to the Fourier slice theorem, the two-dimensional Fourier transform of a projection is identical to the central slice of the three-dimensional (3d) Fourier transform of the index of refraction n-1, normal to the direction of the projection:

$[\eqalignno{{\cal F}_{2{\rm d}} {\cal R} (n-1)& = \textstyle\int \int\limits_{-\infty}^{\infty} \int (n-1) \, {\rm d}z \exp \left [i(-k_x x -k_y y) \right] \, {\rm d}x \, {\rm d}y &\cr &= \left [{\cal F}_{3{\rm d}} (n-1) \right]_{k_z = 0}. & (5)\cr}]$

Hence, every operation carried out in the two-dimensional Fourier space of the projection, notably the multiplication with the propagator P(k_x,k_y), can equally be applied to the corresponding central slice in the three-dimensional (3d) Fourier space. In the linearized approximation of optically weak objects this allows one to invert the order of the projection and propagation:

$[\eqalignno { &{\cal D}_{2{\rm d}} {\cal R} (n-1) = {\cal F}^{-1}_{2{\rm d}} \left (P_{2{\rm d}} {\cal F}_{2{\rm d}} \left\{{\cal R} [n(x,y,z)-1] \right\} \right) &\cr & = {{1}\over{4\pi^2}} \int \int\limits_{-\infty}^{\infty} \Biggl \{\exp[ i \Delta (k^2-k_x^2-k_y^2)^{1/2}] \int \int\limits_{-\infty}^{\infty} \int (n-1)\, {\rm d}z &\cr &\quad\quad\times \exp[-i(k_x x + k_y y)] \, {\rm d}x \, {\rm d}y \Biggr\}&\cr &\quad\quad\times\exp[ i(k_x x + k_y y) ] \, {\rm d}k_x \, {\rm d}k_y &\cr & = {{1}\over{4\pi^2}} \int \int\limits_{-\infty}^{\infty} \int \Biggl\{ \exp[ i \Delta (k^2-k_x^2-k_y^2-k_z^2)^{1/2} ] \int \int\limits_{-\infty}^{\infty} \int (n-1)&\cr &\quad\quad\times \exp[ -i(k_x x + k_y y + k_z z)] \, {\rm d}x \, {\rm d}y \, {\rm d}z\Biggr\}&\cr &\quad\quad\times\exp[ i(k_x x + k_y y) ] \delta_{\rm Dirac}(k_z) \, {\rm d}k_x \, {\rm d}k_y \, {\rm d}k_z &\cr & = {{1}\over{6\pi^3}} \int\limits_{-\infty}^{\infty} \int \int\limits_{-\infty}^{\infty} \int \Biggr\{ \exp[ i \Delta (k^2-k_x^2-k_y^2-k_z^2)^{1/2} ] \int \int\limits_{-\infty}^{\infty} \int (n-1)&\cr &\quad\quad\times \exp[ -i(k_x x + k_y y + k_z z)] \, {\rm d}x \, {\rm d}y \, {\rm d}z\Biggr\}&\cr &\quad\quad\times\exp[ i(k_x x + k_y y + k_z z) ] \, {\rm d}k_x \, {\rm d}k_y \, {\rm d}k_z \, {\rm d} z &\cr & = {\cal R} {\cal F}_{3{\rm d}}^{-1} \left \{P_{3{\rm d}} {\cal F}_{3{\rm d}} [n(x,y,z)-1] \right\} &\cr &= {\cal R} {\cal D}_{3{\rm d}} [n(x,y,z)-1]. &(6)}]$

To avoid misunderstandings, we note that $[\delta _{\rm Dirac}(k_{z})]$ in the expression above denotes the Dirac delta distribution and not the dispersive decrement of the index of refraction. Propagating the three-dimensional object (volume) first and projecting subsequently leads to the same result as conventional two-dimensional propagation of a projection. Moreover, this holds for every projection angle. For tomography, the object is rotated around the y axis, as shown in Fig. 1(a), and projected under a large number of angles $[\alpha _{i}]$ . In Fourier space, these projections correspond to planes with normal vectors in the xz plane sharing one common axis k_y corresponding to the axis of rotation, as also illustrated in Fig. 1(a). Propagating all the projections simultaneously about a distance Δ is equivalent to multiplying $[{\cal F}_{{3{\rm d}}}(n-1)]$ point-wise with the phase function

$[P_{3{\rm d}}(k_x,k_y,k_z,\Delta) = \exp \left[i \Delta (k_0^2-k_x^2 - k_y^2 -k_z^2)^{1/2} \right], \eqno (7)]$

which appeared in equation (6) as the three-dimensional generalization of the conventional two-dimensional propagation kernel $[P_{2{\rm d}}(\Delta)]$ . Again, k_x, k_y and k_z are the components of the sample's Fourier space and not the components of the wavevector $[k_0 = 2\pi/\lambda]$ . Hence, the three-dimensional propagator $[{\cal D}_{3{\rm d}} = {\cal F}^{-1}_{3{\rm d}}[P_{3{\rm d}} {\cal F}_{3{\rm d}}]]$ describes a method to propagate the entire three-dimensional object in the near field and is the most important concept and tool in this work. We denote the output as `propagated object'. While we only use the above formulation of the propagator in this work, the method is not restricted to this particular choice. Alternative formulations including the paraxial approximation of equation (7) (Paganin, 2006) or of the Rayleigh–Sommerfeld/Fresnel–Kirchhoff diffraction integrals (Voelz & Roggemann, 2009 ) could also be used. The three-dimensional nature of the `propagated object', and the fact that the propagation between object and detection planes is now carried out for the entire three-dimensional volume, has some important consequences and is conceptually different from the conventional propagation of two-dimensional wavefronts. One interesting aspect of the three-dimensional propagator approach is given by the fact that phase retrieval can be performed on central planes through the object with arbitrary orientation, including in particular planes orthogonal to the tomographic rotation axis. In view of computational speed, we found that the single three-dimensional propagation outperformed the conventional N two-dimensional propagations, in particular for data sets with a high number of projections N.

Figure 1
(a) Schematic of holographic phase contrast imaging. A plane wave illuminates the object, leading to an exit wave in the xy plane as given by the projected optical indices, followed by free space propagation resulting in holographic phase contrast formation as recorded in the detection plane at distance Δ. The projection $[{\cal R}]$ and propagation $[{\cal D}]$ process can be described theoretically in Fourier space by multiplying the central slices of the object's Fourier transform with a radially symmetric phase factor. (b) For weakly interacting objects, the projection operation $[{\cal R}]$ and the propagation operation $[{\cal D}]$ can be permuted, allowing for the propagation of the three-dimensional index of refraction.

Next, we show how the three-dimensional propagator can be used advantageously in iterative phase retrieval. To this end, we do not restrict ourselves to small propagation distances, i.e. validity of the TIE regime (Paganin, 2006), but consider the problem for general Fresnel numbers including the holographic regime. We remain, however, within the weak object approximation, as stated above. Note that, strictly speaking, the linearization with respect to the object needed here is valid beyond the strict weak phase object [see for example the slowly varying phase condition introduced in Guigay (1977 ) and Gureyev & Nesterets (2015 )]. Since in an experiment only the intensity $[I \propto |\Psi_\Delta|^2 = \Psi_\Delta \Psi_\Delta^*]$ of the propagated exit field can be measured, we do not have direct access to the propagated three-dimensional object. However, we can apply an inverse Radon transform $[{\cal R}^{-1}]$ on the intensities I_i, combining all $[i = 1 \ldots N]$ projections to a three-dimensional `intensity field' $[I_{3{\rm d}}: = {\cal R}^{-1} (I_i / |\Psi_0|^2)]$ . Using equation (6), the usual decomposition of the measured intensities into the different components of the hologram is given as

$[\eqalignno{&I_{3{\rm d}} - {\cal R}^{-1}(1) =&\cr & i k_0 \left [\underbrace{\exp(-i k_0 \Delta) {\cal D}_{3{\rm d}} (n-1)}_{{\rm object}} - \underbrace{\exp(i k_0 \Delta) {\cal D}_{3{\rm d}}^* (n-1)^*}_{{\rm twin\hbox{-}object}} \right]&\cr &\quad + \underbrace{k_0^2 {\cal R}^{-1} |{\cal R} (n-1)|^2}_{\ll 1}, & (8)\cr}]$

formulating the well known twin-image problem of holography in three dimensions. The expression shows that the three-dimensional `intensity field' I_3d is consistent, since it can be written as the sum of the propagated object and the complex conjugated `twin object' (n-1)^* propagated about the distance $[-\Delta]$ , since $[{\cal D}(\Delta)^{*} = {\cal D}(-\Delta)]$ . Tomographic consistency, described in detail by Helgason and Ludwig (Helgason, 1965; Ludwig, 1966; Kudo & Saito, 1991), states that projections from the same object are not independent from each other. The Fourier slice theorem already shows that all projections share a common intersection line in Fourier space. For objects of finite size, the central slices in Fourier space can be seen as broadened by a convolution with the object's support. This introduces a dependence between all slices. Explicitly, one cannot change the value at any given point in Fourier space without violating consistency for an object of finite size. From equation (8) we can infer that the twin-object/phase problem does not introduce such tomographical inconsistencies. In return, this property cannot be exploited to improve phase retrieval directly. In particular, we expect only limited improvements for three-dimensional phase retrieval by deterministic one-step algorithms like a direct CTF inversion, where the three-dimensional inversion has been used previously (Frank & Penczec, 1995; Cloetens et al., 1997) if carried out in three dimensions rather than two dimensions, even if the geometry of the three-dimensional problem is different from two-dimensional imaging. With the definition $[\exp(i\chi): = \exp(-ik_{0}\Delta)P_{{3{\rm d}}}]$ and $[(n-1) = -\delta+i\beta]$ , equation (8) is rewritten in Fourier space as

$[\eqalignno{{\cal F}_{3{\rm d}} \left [{{I_{3{\rm d}} - {\cal R}^{-1}(1)}\over{k_0}} \right] &\simeq \underbrace{ i \left [\exp(i \chi) - \exp(-i \chi) \right] }_{-2 \sin(\chi)} {\cal F}_{3{\rm d}} (-\delta)&\cr &\quad- \underbrace{ \left [\exp(i\chi) + \exp(-i \chi)\right]}_{2 \cos(\chi)} {\cal F}_{3{\rm d}}(\beta), &(9)\cr}]$

illustrating the structure of the CTF in three dimensions. The CTF_3d with zero values on spherical shells does not contain all information about the object. Therefore, one has to resort to further constraints in phase retrieval, i.e. a priori information, such as positivity, range restrictions of the object functions, or known support of the object. However, we have shown previously that the performance of iterative phase retrieval with constraints such as positivity can be dramatically improved if tomographic consistency is enforced (Ruhlandt et al., 2014). If the phase retrieval is carried out in three dimensions according to the concept presented, consistency is guaranteed automatically.

3. Simulation and numerical results

To validate and to illustrate the concept introduced above, we show an exemplary simulation. A sample was designed consisting of 30 spheres with a diameter between 20 and 50 px (px stands for pixel) distributed randomly on a 256³ voxel grid. For each sphere, a combination of $[0.001 \leq k_0\delta \leq 0.01]$ and $[0 \leq k_0 \beta \leq 0.0007]$ was selected randomly, resulting in a maximum phase shift of 1.5 rad and a maximum absorption of 6% in the projections.¹ The object was projected to 402 equidistant angles in the range 0° $[\leq \alpha \le]$ 180° to satisfy the angular sampling criterion. To simulate the intensity data I_meas in the detection plane, all projections were propagated about a distance $[\Delta = 2000]$ px with a wavelength of $[\lambda = 0.01]$ px, resulting in a Fresnel number of F₁ = 0.05 for one pixel. This choice corresponds for example to a relevant experimental setting in X-ray phase contrast tomography with photon energy of 12.4 keV, an effective pixel size of 10 nm and an effective sample-to-detector distance of 20 µm. The phase retrieval algorithm cycles between the object plane and detection plane, always propagating the wavefield Φ. In the detection plane, the modulus constraint

$[\tilde{\Phi}^\prime = \Phi_{\rm meas} {{\tilde{\Phi}}\over{|\tilde{\Phi}|}} \eqno (10)]$

is applied, setting the modulus of the wavefield $[\tilde{\Phi}^\prime]$ to the measured values $[\Phi_{\rm meas} = (I_{\rm meas})^{1/2}]$ while keeping the phase information. In the object plane we demand positivity for both the δ and β parts of the index of refraction which follows directly from its definition and the projection approximation:

$[\Phi^\prime = 1- |{\rm Re}(\Phi_g^\prime)-1| - i|{\rm Im}(\Phi_g^\prime)|. \eqno (11)]$

Additionally, we enforced $[\beta \leq 0.1 \times \delta]$ , which is justified for a large class of hard X-ray experiments and samples, including in particular biological tissues. This constraint is implemented in the following way: we first set $[\delta ^\prime = \beta]$ in pixels where $[\beta \,\gt\, \delta]$ and only afterwards set $[\beta ^\prime = 0.1 \times \delta ^\prime]$ where $[\beta\ \gt]$ $[ 0.1 \times \delta]$ . 1000 iterations of this `soft coupling' scheme were carried out for (i) the conventional method with the phase information of each detector image retrieved individually followed by a filtered back-projection (`two-dimensional reconstruction') and (ii) for the three-dimensional volume starting with a filtered back-projection of the $[\Phi_{\rm meas}]$ . Typical results are depicted in Fig. 2, showing clearly the better quality of the three-dimensional approach. In the supporting information for this paper, additional simulation results are presented for noisy intensity data (Poissonian noise corresponding to 10 000 ph px⁻¹; ph stands for photon) as well as results without the `soft coupling' and thus without any assumption on the object. While the overall quality of these reconstructions is not as good as in Fig. 2, both cases still show a major improvement of the three-dimensional approach compared to the two-dimensional reconstruction and illustrate the robustness of the method. The entire three-dimensional phase retrieval process was carried out within about 2 min using MATLAB and an NVIDIA GTX Titan GPU. Example code in MATLAB/Octave for the simulation and the explained three-dimensional reconstruction is available as supporting information.

Figure 2
Comparison of the phase shifts reconstructed using the conventional two-dimensional propagation and the three-dimensional propagation. The upper row shows a typical projection of the phantom/reconstructed volume, the lower row shows an x–z slice, normal to the axis of rotation of the volume. In the left column the original phantom is depicted, the central column contains the results of the three-dimensional propagation method presented here and the right column shows the conventional results. The scale bar is the same for all images; the colour bars apply to the corresponding rows.

Based on the simulations, we can expect that the approach also performs well on noisy experimental data. To this end, we have tested the algorithm on the data presented previously in Bartels et al. (2012) and Ruhlandt et al. (2014), corresponding to holographic projection images of freeze-dried Deinococcus radiodurans bacteria, dispersed on ultra-thin Si₃N₄ membranes and recorded using 13.8 keV radiation exiting from a waveguide [see Bartels et al. (2012) for experimental details]. The tomographic scan (single defocus distance data set) comprised 83 projection angles distributed over 162°. The reconstructed phase of the bacteria as resulting from the three-dimensional propagator approach in combination with ART is displayed for a representative slice through the object (Fig. 3). For comparison, reconstructions are shown as obtained by the previously presented algorithms, the modified HIO (mHIO) (Giewekemeyer et al., 2011; Bartels et al., 2012) and the IRP (Ruhlandt et al., 2014). A major difference between mHIO on the one side, and IRP and the scheme reported here (three-dimensional propagation) on the other is that mHIO needs additional support constraints for phase retrieval, while the latter do not. Note that, despite this additional constraint, mHIO leads to artifacts such as the increase of density towards the top and bottom corners, while the reconstruction of IRP and the method presented here yield a much more plausible density distribution. This advantage of IRP over mHIO has been stressed before (Ruhlandt et al., 2014), but came at the expense of significant numerical complexity, while the three-dimensional propagation reaches a quality comparable to IRP in a fraction of the computation time.

Figure 3
Comparison of different reconstructions from the same experimental data. All images show the same slice through the reconstructed three-dimensional volumes of a Deinococcus radiodurans bacteria. (a) illustrates the `conventional' mHIO reconstruction reported before (Bartels et al., 2012

; Ruhlandt et al., 2014

) where additional support constraints have been used. (b) The IRP result of the same data shows a much more homogeneous signal distribution without the use of support information. The method described here leads to the result depicted in (c) showing a comparable reconstruction quality as IRP but calculated in a fraction of the time.

4. Summary and conclusions

While three-dimensional reconstruction in far-field diffraction has long been known to aid phase retrieval (Marchesini et al., 2003 ), the generalization to the near-field case was less obvious since the measurements (near-field patterns) cannot in general be considered to be sub-manifolds of a three-dimensional reciprocal space of the sample. Therefore, an algebraic tomography scheme was previously proposed for the purpose of three-dimensional reconstruction of arbitrary objects (Kostenko et al., 2013; Ruhlandt et al., 2014), however at the cost of numerous iterations of performing Radon and inverse Radon transforms.

In this work, we have introduced a new approach to near-field tomographic phase retrieval in the limit of an optically weak object. The usual sequence of first retrieving the phase information of all projections individually followed by an inverse Radon transform is inverted. A three-dimensional volume is computed from the measurements, followed by iterative propagations and the application of constraints in three dimensions. Since the inverse Radon transform is performed only once, a tremendous enhancement in reconstruction speed is obtained with respect to previous combined schemes of iterative phase retrieval and tomographic reconstructions. At the same time, the three-dimensional propagation method preserves the essential advantages of tomographic consistency, which is intrinsically enforced by the three-dimensional scheme and was found to stabilize phase retrieval with otherwise under-determined data.

Even though the three-dimensional propagation does not lead to major improvements in the quality of direct CTF phase retrieval, it can be useful to investigate the influence of regularization parameters directly in three dimensions.

Another immediate advantage of the three-dimensional scheme not yet exploited here is directly obvious: in three dimensions, the formulation of entirely new and more powerful constraints is possible, since they can be applied directly on the object rather than its projections. This applies for example to a much more accurate and constraining support determination, to positivity as well as to sparsity or to regularization procedures. In general, constraints can not only be formulated in physically correct and direct terms in three dimensions, but can possibly also be applied to a higher fraction of voxels. Finally, the measurement scheme could easily be generalized to several detection planes, and – somewhat less straightforward – also to more complex illumination wavefields than plane waves.

The theory presented here is based on the assumption of an optically weak object. However, better phase retrieval results were obtained in three dimensions even if the maximum phase shift in some projections was as high as $[\simeq]$ 1.5 rad, as shown in the example of Fig. 2.

In future work, extensions to optically thick samples could be investigated. The corresponding exit fields of a sample may be expected to fulfil generalized consistency criteria which might be exploitable, and problems of phase wrapping may also be treated much better if the object is represented in three dimensions throughout the phase retrieval process.

More straightforward extensions of the present work will extend the numerical comparison between two-dimensional and three-dimensional phase retrieval and establish how the reconstruction quality gain depends on all experimental parameters, including number of projections, Fresnel number, added noise or even systematically perturbed data. To this end, the MATLAB example code (available as supporting information) may be helpful.

Supporting information

Supplementary figures. DOI: https://doi.org/10.1107/S2053273315022469/mq5036sup1.pdf

Example code in MATLAB/Octave. DOI: https://doi.org/10.1107/S2053273315022469/mq5036sup2.txt

Footnotes

¹This value clearly violates the assumption of a weakly interacting object, but illustrates that the advantage of the three dimensions persists even to larger phase shifts.

Acknowledgements

We thank Johannes Hageman for helpful discussions and in particular for his advice on numerical issues, and acknowledge support by the Collaborative Research Center 755 `Nanoscale Photonic Imaging' of the German Science Foundation (DFG).

References

Allen, L. J. & Oxley, M. P. (2001). Opt. Commun. 199, 65–75. Web of Science CrossRef CAS Google Scholar
Bartels, M., Priebe, M., Wilke, R. N., Krüger, S., Giewekemeyer, K., Kalbfleisch, S., Olendrowitz, C. C., Sprung, M. & Salditt, T. (2012). Opt. Nanoscopy, 1, 10. CrossRef Google Scholar
Bronnikov, A. V. (1999). Opt. Commun. 171, 239–244. Web of Science CrossRef CAS Google Scholar
Bronnikov, A. V. (2002). J. Opt. Soc. Am. A, 19, 472–480. Web of Science CrossRef Google Scholar
Chapman, H. N., Barty, A., Marchesini, S., Noy, A., Hau-Riege, S. P., Cui, C., Howells, M. R., Rosen, R., He, H., Spence, J. C. H., Weierstall, U., Beetz, T., Jacobsen, C. & Shapiro, D. (2006). J. Opt. Soc. Am. A, 23, 1179–1200. Web of Science CrossRef Google Scholar
Cloetens, P., Ludwig, W., Baruchel, J., Van Dyck, D., Van Landuyt, J., Guigay, J. P. & Schlenker, M. (1999). Appl. Phys. Lett. 75, 2912–2914. Web of Science CrossRef CAS Google Scholar
Cloetens, P., Pateyron-Salom, M., Buffire, J. Y., Peix, G., Baruchel, J., Peyrin, F. & Schlenker, M. (1997). J. Appl. Phys. 81, 5878–5886. CrossRef CAS Web of Science Google Scholar
Frank, J. & Penczec, P. (1995). Optik, 98, 125–129. CAS Google Scholar
Gerchberg, R. W. & Saxton, W. O. (1972). Optik (Jena), 35, 237–246. Google Scholar
Giewekemeyer, K., Krüger, S. P., Kalbfleisch, S., Bartels, M., Beta, C. & Salditt, T. (2011). Phys. Rev. A, 83, 023804. Web of Science CrossRef Google Scholar
Gordon, R., Bender, R. & Herman, G. T. (1970). J. Theor. Biol. 29, 471–476. CrossRef CAS PubMed Web of Science Google Scholar
Groso, A., Abela, R. & Stampanoni, M. (2006). Opt. Express, 14, 8103–8110. Web of Science CrossRef PubMed CAS Google Scholar
Guigay, J. P. (1977). Optik, 49, 121–125. Google Scholar
Gureyev, T. E., Davis, T. J., Pogany, A., Mayo, S. C. & Wilkins, S. W. (2004). Appl. Opt. 43, 2418–2430. Web of Science CrossRef PubMed Google Scholar
Gureyev, T. E. & Nesterets, Ya. I. (2015). Linear algorithms for phase retrieval in the Fresnel region: validity conditions. https://arxiv.org/abs/1503.06949. Google Scholar
Helgason, S. (1965). Acta Math. 113, 153–180. CrossRef Google Scholar
Hofmann, R., Moosmann, J. & Baumbach, T. (2011). Opt. Express, 19, 25881–25890. Web of Science CrossRef CAS PubMed Google Scholar
Hong, Y., Zhang, K., Wang, Z., Zhu, Z., Zhao, X., Huang, W., Yuan, Q., Zhu, P. & Wu, Z. (2012). Chin. Phys. B, 21, 104202. Web of Science CrossRef Google Scholar
Kak, A. C. & Slaney, M. (1988). Principles of Computerized Tomographic Imaging. New York: IEEE Press. Google Scholar
Kostenko, A., Batenburg, K. J., King, A., Offerman, S. E. & van Vliet, L. J. (2013). Opt. Express, 21, 12185–12196. Web of Science CrossRef PubMed Google Scholar
Krenkel, M., Bartels, M. & Salditt, T. (2013). Opt. Express, 21, 2220–2235. Web of Science CrossRef CAS PubMed Google Scholar
Kudo, H. & Saito, T. (1991). J. Opt. Soc. Am. A, 8, 1148–1160. CrossRef Web of Science Google Scholar
Langer, M., Cloetens, P., Pacureanu, A. & Peyrin, F. (2012). Opt. Lett. 37, 2151–2153. CrossRef PubMed Google Scholar
Latychevskaia, T. & Fink, H.-W. (2007). Phys. Rev. Lett. 98, 233901. Web of Science CrossRef PubMed Google Scholar
Ludwig, D. (1966). Commun. Pure Appl. Math. 19, 49–81. CrossRef Google Scholar
Marchesini, S., Chapman, H., Hau-Riege, S., London, R., Szoke, A., He, H., Howells, M., Padmore, H., Rosen, R., Spence, J. & Weierstall, U. (2003). Opt. Express, 11, 2344–2353. Web of Science CrossRef PubMed CAS Google Scholar
Miao, J., Förster, F. & Levi, O. (2005). Phys. Rev. B, 72, 052103. Web of Science CrossRef Google Scholar
Miao, J., Hodgson, K. O. & Sayre, D. (2001). Proc. Natl Acad. Sci. USA, 98, 6641–6645. Web of Science CrossRef PubMed CAS Google Scholar
Paganin, D., Mayo, S. C., Gureyev, T. E., Miller, P. R. & Wilkins, S. W. (2002). J. Microsc. 206, 33–40. Web of Science CrossRef PubMed CAS Google Scholar
Paganin, D. & Nugent, K. A. (1998). Phys. Rev. Lett. 80, 2586–2589. Web of Science CrossRef CAS Google Scholar
Paganin, D. M. (2006). Coherent X-Ray Optics. New York: Oxford University Press. Google Scholar
Ruhlandt, A., Krenkel, M., Bartels, M. & Salditt, T. (2014). Phys. Rev. A, 89, 033847. Web of Science CrossRef Google Scholar
Thibault, P. (2007). Algorithmic Methods in Diffraction Microscopy. Ithaca: Cornell University. Google Scholar
Turner, L., Dhal, B., Hayes, J., Mancuso, A., Nugent, K., Paterson, D., Scholten, R., Tran, C. & Peele, A. (2004). Opt. Express, 12, 2960–2965. Web of Science CrossRef PubMed CAS Google Scholar
Voelz, D. G. & Roggemann, M. C. (2009). Appl. Opt. 48, 6132–6142. Web of Science CrossRef PubMed Google Scholar
Wilkins, S. W., Gureyev, T. E., Gao, D., Pogany, A. & Stevenson, A. W. (1996). Nature (London), 384, 335–338. CrossRef CAS Web of Science Google Scholar
Wu, X., Liu, H. & Yan, A. (2005). Opt. Lett. 30, 379–381. Web of Science CrossRef PubMed CAS Google Scholar