Using sparsity information for iterative phase retrieval in x-ray propagation imaging

A. Pein; S. Loock; G. Plonka; T. Salditt

doi:10.1364/OE.24.008332

1. Introduction

The first observation of propagation-based phase contrast in the x-ray regime some 20 years ago at the European Synchrotron Radiation Facility (ESRF, Grenoble) [1,2] has started a considerable activity in developing this lens-less full field imaging approach to resolve micro- and nano-structures in the interior of specimens [3–5]. Using the divergent wavefronts emitted from the quasi-point source of an x-ray waveguide channel has enabled a resolution down to 20 nm [6]. The phase retrieval problem in this near field setting can be regarded as better posed compared to the far field counterpart, see also [7] for an illuminating novel view on uniqueness in near field phase retrieval. Nevertheless, a persisting challenge remains in the phase retrieval of images recorded at a single distance and without translations of the object in the beam. Without such added diversity in the data, phase and amplitude of an object’s wave cannot be recovered independently, unless ‘strong’ a priori information, such as compact support, is available. Since many relevant samples are extended, it is generally believed that phase information from single distance measurements can only be recovered for objects where phase and amplitude variations are strongly coupled, or for objects which interact only by either phase or amplitude, i.e. so-called pure phase or pure absorption objects.

As an introduction and for notational clarity, let us briefly recall the imaging process. When an x-ray wave field interacts with a specimen not only its amplitude is modified but its phase distribution as well. Both effects are determined by the specimen’s internal structure, which is reflected in its complex refractive index for x-rays, see e.g. [8],

n (x, y, z) = 1 - δ (x, y, z) + i β (x, y, z),

where β characterizes the absorption and δ the induced phase shift. By allowing the transmitted wave field (the so-called exit wave field) to freely propagate between the specimen and the detector (at distance z), the induced phase variations are transformed into intensity variations in such a way that the detected intensity distribution is a function of both the amplitude and the phase of the exit wave field, see Fig. 1(a) for the experimental setup. By this means, x-ray propagation imaging is able to extract and visualize information about both β and δ of the analyzed specimen. In particular, for thin biological samples of light elements, phase contrast is indispensable, as δ becomes about three orders of magnitude larger than β and absorption contrast is negligibly low [9].

Fig. 1 (a) Setup of an x-ray propagation imaging experiment. The specimen is illuminated by a coherent x-ray probe. Traversing the object, the waves are modified according to the complex refractive index and thereafter the exit wave (object plane) propagates to a detector at a distance z behind the specimen (detector plane). The phase contrast arises from the free space propagation without further optical elements needed. (b) Illustration of the shearlet constraint, incorporated into the ER algorithm. The current approximation is alternately propagated between the object and the detector plane where it is forced to fulfill the shearlet constraint and the modulus constraint, respectively. Here, for the general case of a mixed object, where shearlet thresholding is performed on the phase and the amplitude of the approximation separately. SH and SH⁻¹ denote the discrete shearlet transform and its inverse.

Download Full Size | PDF

More precisely, within the Fresnel regime, the connection between the wave field in the detector plane ψ(x,y,z) and the exit wave field ψ(x,y,0), under the paraxial and the projection approximation and assuming monochromatic plane wave illumination in z-direction, reads [10]

\begin{array}{l} ψ (x, y, z) = D_{z} [ψ (x, y, 0)] \\ : = \exp (i k z) ℱ_{⊥}^{- 1} [\exp (- i z \frac{k_{x}^{2} + k_{y}^{2}}{2 k}) ℱ_{⊥} [ψ (x, y, 0)]] \\ = \exp (i k z) \exp (i k τ) ℱ_{⊥}^{- 1} [\exp (- i z \frac{k_{x}^{2} + k_{y}^{2}}{2 k}) ℱ_{⊥} [\exp (- i k \int_{- τ}^{0} (δ (x, y, z) - i β (x, y, z)) d z)]] . \end{array}

Here, ℱ_⊥ and $ℱ_{⊥}^{- 1}$ denote the Fourier transform in lateral coordinates and its inverse, τ is the specimen’s maximal thickness and k = (k_x,k_y,k_z) is the wave vector with modulus $k = \frac{2 π}{λ}$ . The operator D_z is the so-called Fresnel propagator. However, detectors can only capture the time-averaged squared amplitude, i.e. the intensity, of the propagated wave field

I (x, y, z) = {| ψ (x, y, z) |}^{2} = {| D_{z} [ψ (x, y, 0)] |}^{2} .

The loss of phase information impedes the reconstruction of the exit wave from the measurement by simply applying the inverse Fresnel propagator D₋_z (non - crystallographic phase problem). Many different methods to approach this inverse problem have been developed throughout the years [11–13].

For phase retrieval in (far field) coherent diffractive imaging (CDI) iterative algorithms that make use of a priori information of the signal to be reconstructed represent an established strategy [14–17]. Likewise, they have become a powerful tool in (near field) x-ray propagation imaging in the holographic regime [18–21]. The a priori information and/or information from intensity measurements at a single or several defocus distances are incorporated into constraint maps that, in turn, are combined into iterative schemes. For example, in [21], a regularization strategy for phase retrieval is proposed that employs sparsity in an orthogonal wavelet frame. Further remarks and references on application of wavelet sparsity constraints can be found in [27]. Seminal algorithms, which can be applied to both the conventional far field and the optical near field case, are the Gerchberg-Saxton algorithm [22] and Fienups’s Error Reduction (ER) algorithm [23] and Hybrid Input-Output (HIO) algorithm [24]. Aiming at more stable and fast phase retrieval algorithms new relaxation strategies have been developed that led to the Hybrid Projection Reflection (HPR) algorithm [25] and to the Relaxed Averaged Alternating Reflection (RAAR) algorithm [26]. In this paper, we are interested in near field experiments with a single distance measurement and we will use the RAAR method for phase reconstruction that is especially suitable for noisy input data. The standard approach in this setting is to use a priori knowledge about the compact support of the specimen in terms of a corresponding support constraint. However, in case the object’s support cannot be identified or it occupies the whole field of view, alternative constraints of similar strength need to be found.

Recently, Loock and Plonka [27] proposed a new constraint for iterative x-ray phase retrieval in the Fresnel regime, which uses sparsity information of the wavefront, and they presented promising numerical results with this method. In this paper, we will extend this approach to the physical situation of a complex-valued exit wave field and apply the new sparsity constraint separately to its amplitude and phase. Further, we will demonstrate its applicability to experimental data. We start by introducing the idea behind this new constraint (Section 2), followed by the presentation of simulation results (Section 3). In Section 4 the method is applied to experimental data and we end with a short conclusion in Section 5.

2. The shearlet constraint* for iterative phase retrieval*

The a priori information that an object transmission function (exit wave) may be sparsely represented in a certain representation system can be advantageously used for near field phase retrieval. To this end, sparse may mean that only a few expansion coefficients are nonzero, or in a more relaxed sense that only few expansion coefficients contain the ’essential’ information. In [27], it was in particular shown that the so-called shearlet system is very suitable for this purpose, as demonstrated by reconstructing simulated data of real-valued object transmission functions in the presence of noise. Two-dimensional (2D) shearlet systems, introduced in [28] and [29], are constructed by applying scaling, translation and shearing operators to a set of 2D wavelet-type generating functions. By these means, the system’s elements are defined at various scales, locations and orientations, allowing them to efficiently capture anisotropic image elements like discontinuity curves. The associated shearlet transform maps a function f ∈ L²(ℝ²) to the sequence of all shearlet coefficients, defined as the L²-scalar product between the function and each element of the corresponding shearlet system. In [30] detailed information on the theory and applications of shearlets are given. In particular, under certain conditions the shearlet system forms a frame for the L²(ℝ²), providing stable reconstructions of elements f ∈ L²(ℝ²) from their shearlet coefficients (inverse shearlet transform). While the first shearlet systems had been constructed using band-limited generating functions, newer developments also focus on shearlet frames whose elements have compact support in space domain, see [31].

The class of so-called cartoon-like images is mathematically defined as a subspace of L²(ℝ²) that contains all compactly supported functions that are twice continuously differentiable apart from piecewise C²-discontinuity curves, see [32]. Many images encountered in physical applications, like the phase or the amplitude of exit wave fields in x-ray propagation imaging, can be expected to fall into this image class. Let f be a cartoon-like image and let f_N be the approximation that takes the N largest shearlet coefficients in magnitude into account for the reconstruction via the inverse transform. Then, shearlet theory ensures the following approximation rate

{‖ f - f_{N} ‖}_{L^{2}}^{2} \leq C N^{- 2} {(\log N)}^{3}, as N \to \infty, C \in ℝ,

for band-limited shearlets [33] as well as for compactly supported shearlets [34], which is nearly optimal (apart from the log-factor). Hence, shearlet frames provide sparse representations of cartoon-like functions. In other words, in many relevant experimental situations encountered in x-ray propagation imaging, we may assume the a priori knowledge that the amplitude/phase of the exit wave field is sparsely represented in a suitable shearlet frame. As shown in [27], this information can be used as a constraint for iterative phase retrieval algorithms.
Let $ψ \in ℂ^{N_{x} \times N_{y}}$ be the discrete, complex-valued exit wave field and let $D_{z} ψ \in ℂ^{N_{x} \times N_{y}}$ be the discrete, propagated field in the detector plane, i.e., D_z denotes the discrete Fresnel transform. Further, SH denotes the discrete shearlet transform that transfers an image of size N_x×N_y into a vector of shearlet coefficients of length N_x · N_y · S, where S > 1 indicates the redundancy of the transform that depends on the number of different scales and shearings that are taken into account by the discrete shearlet frame. Further, SH⁻¹ denotes the discrete inverse shearlet transform, which is a pseudo inverse that can be efficiently computed by conjugate gradient methods [35]. In this work, the transforms are computed with the ShearLab 3D toolbox (freely accessible from shearlab.org) that employs a system of compactly supported shearlets [36].
The basic idea of the shearlet (sparsity) constraint for complex-valued exit waves is to enforce a sparse shearlet representation of the amplitude and phase of the exit wavefront by applying a shearlet shrinkage at each iteration step. Let us explain the approach in more detail within the framework of the Error Reduction algorithm, see Fig. 1(b). Starting with an initial guess ψ₀, the current approximation ψ_n is alternately propagated between the object and the detector plane, where it is forced to fulfill the object plane constraint (here: the shearlet constraint) and the modulus constraint, respectively. The modulus constraint enforces the approximation to be in agreement with the measurement $m = \sqrt{I} \in ℝ^{N_{x} \times N_{y}}$ , see Eq. (3), on a grid in the detector plane. The corresponding projection onto the set of digital images that agree with the measurement reads, see e.g. [27],
$P_{M} ψ = D_{- z} g with g (i, j) : = {\begin{array}{l} m (i, j) \frac{D_{z} ψ (i, j)}{| D_{z} ψ (i, j) |} & , | D_{z} ψ (i, j) | \neq 0 \\ m (i, j) & , | D_{z} ψ (i, j) | = 0 \end{array} for (i, j) \in Ω,$ with Ω = {1,…,N_x} × {1,…,N_y}.
In the object plane we employ the shearlet constraint. For that purpose, we apply first the shearlet transform separately to the amplitude and to the phase of ψ, then cut off small shearlet coefficients by a threshold operator before applying the inverse shearlet transform. In practice, instead of such a strict cut off, we use a soft threshold operator for the shearlet coefficients $c = {(c_{j})}_{j = 1 \dots N_{x} \cdot N_{y} \cdot S}$ , defined as
${(T_{^{θ}} c)}_{j} : = {\begin{array}{l} c_{j} - θ & , c_{j} > θ \\ c_{j} + θ & , c_{j} < - θ \\ 0 & , | c_{j} | < θ \end{array},$ with threshold parameter θ. The choice of the soft threshold operator is motivated in [27], since it coincides with the so-called proximity operator corresponding to the l₁-norm of the sequence of shearlet coefficients. At the same time, it causes a smoothing effect, similarly as for wavelet denoising, see e.g. [37]. The shearlet constraint operator can be written as $P_{S ℋ}^{θ_{a}, θ_{p}} ψ : = S ℋ^{- 1} T_{θ_{a}} S ℋ | ψ | \cdot \exp (i \cdot S ℋ^{- 1} T_{θ_{p}} S ℋ φ (ψ)),$ where φ(ψ) denotes the phase of the complex-valued wave field ψ. Different threshold parameters θ_a,θ_p for the amplitude and phase can be taken into account.
For the reconstructions presented in this work we use the more complex Relaxed Averaged Alternating Reflection (RAAR) algorithm [26]. Formulated as a fixed-point iteration with the shearlet constraint as object plane constraint, this algorithm reads
$ψ_{n + 1} = (\frac{α}{2} (R_{S ℋ}^{θ_{a}, θ_{p}} R_{M} + Id) + (1 - α) P_{M}) ψ_{n},$ with relaxation parameter α and reflectors R_M := 2P_M − Id, $R_{S ℋ}^{θ_{a}, θ_{p}} : = 2 P_{S ℋ}^{θ_{a}, θ_{p}} - Id$ (Id being the identity operator). A mathematical treatment regarding convergence algorithm within this setting may be obtained by extending the approach in [27]. Here, we focus on the evaluation of the numerical results.
Equation (7) can be applied to the general case of a mixed object that causes both phase shifts and absorption of the illuminating wave field, the so-called probe p. As we are dealing with monochromatic plane wave illumination we have p = 1. For so-called pure phase objects, the absorption of the probe can be neglected, i.e. the amplitude of the exit wave can be set to one and Eq. (7) becomes $P_{S ℋ}^{θ_{p}} ψ : = \exp (i \cdot S ℋ^{- 1} T_{θ_{p}} S ℋ φ (ψ))$ . Likewise, for pure absorption objects, where the phase shift can be neglected, it reads $P_{S ℋ}^{θ_{a}} ψ : = S ℋ^{- 1} T_{θ_{a}} S ℋ | ψ |$ , i.e. the phase is set to zero. Another constraint based on a priori information, which is frequently used for phase retrieval in x-ray propagation imaging experiments, concerns the range of values that the phase and amplitude of the exit wave can obtain. Within the projection approximation we know that the amplitude of the exit wave field can only be diminished due to absorption inside the object, i.e. we have the a priori information |ψ| ≤ 1. Furthermore, we know that with x-ray energies the real part of the refractive index is smaller than one, i.e. δ ≥ 0 (cf. Eq. (1)) and we only observe negative phase shifts φ(ψ) ≤ 0 (cf. Eq. (2)). This knowledge can define an object plane constraint by itself (the so-called range constraint) or it can be additionally combined with other constraints. Table 1 summarizes the effects of all constraints used within this paper for different types of objects.
Table 1. Effects of different object plane constraints for mixed objects, pure phase objects and pure absorption objects. |ψ|, $φ (ψ) \in ℝ^{N_{x} \times N_{Y}}$ denote the amplitude and phase of the discrete exit wave field. For the support & range constraint, it is known that the compact support of the exit wave field is contained in some region D ⊂ Ω.
View Table | View all tables in this article
3. Simulation
Figure 2 shows a simulated x-ray propagation imaging experiment for an object which causes both phase shifts and absorption of the probe (denoted here as mixed object). Figures 2(a) and 2(b) depict the phase shift and amplitude distribution of the exit wave field behind the object, which are structurally completely uncorrelated. The simulation’s setup parameters are: δ = 1.6 · 10⁻⁶, β = 8 · 10⁻⁷, τ = 10 µm, wavelength λ = 10⁻¹⁰ m and Fresnel number F = 4 · 10⁻³. Figure 2(c) shows the simulated intensity measurement, degraded by artificial Poisson noise with 50 photons per pixel. For reconstructing the exit wave field we use the RAAR algorithm (relaxation parameter α = 0.55) with the holographic reconstruction (simple back-propagation D₋_z of m) as initial guess. Figures 2(d)–2(f) show the reconstructions obtained after 100 iteration steps with different object plane constraints, namely the range constraint, the support & range constraint and the shearlet & range constraint. The support considered for the support & range constraint is indicated by the dashed black frames in Figs. 2(a) and 2(b). For the shearlet & range constraint we use a shearlet system of 5 scales (81 shearlets in total not counting translations) and a fixed threshold parameter of θ_a = θ_p = 0.005 for both the thresholding of the phase and the amplitude.

Fig. 2 Simulated x-ray propagation imaging experiment. (a), (b) Phase and amplitude of the complex-valued exit wave field in the object plane (512×512 pixel). The dashed black frames indicate the support considered for the support & range constraint (341×256 pixel). (c) Simulated intensity measurement with Fresnel number F = 4 10⁻³ and artificial Poisson noise (50 photons per pixel). The coarse outlines of the fish ·and the elephant are still visible. (d), (e), (f) Exit waves reconstructed with the RAAR algorithm (after 100 iterations) and different object plane constraints. The phase is shown in the upper row and the amplitude in the lower row. The amplitudes obtained with the range constraint and the support & range constraint are displayed after only 3 iterations since the corresponding error increases with the number of iteration steps, see (g). (g) RMS error decay of the phase and amplitude from the three different reconstructions, calculated according to Eq. (9) and (10), while taking the region inside the box support into account.
Download Full Size | PDF
The reconstructed phase and amplitude obtained with the shearlet & range constraint in Fig. 2(f) show the strongest resemblance to the original exit wave (cf. Figs. 2(a) and 2(b)), compared to the results obtained with the range constraint in Fig. 2(d) and the support & range constraint in Fig. 2(e). For instance, unlike in Figs. 2(d) and 2(e), the small lines on the fish’s pectoral fin and on the elephant’s trunk are clearly visible in Fig. 2(f). Furthermore, Fig. 2(f) shows a smaller noise level and a better agreement in the range of pixel values with the original exit wave field (cf. Figs. 2(a) and 2(b)) in comparison to Figs. 2(d) and 2(e). The reconstructions in Figs. 2(d) and 2(e) resemble each other strongly, i.e. the additional support information does not seem to be of particular importance in this setting. This may be attributed to the fact that the restriction of the exit wave to the box support is not a very strong constraint as the diffraction fringes in the measurement merely exceed the region of the box.
In Fig. 2(g) we have plotted the root mean square (RMS) error decay of the phase and the amplitude for the three different reconstructions, calculated by
${RMS}_{phase}^{n} = {[\frac{1}{N_{x} N_{y}} \sum_{n_{x} = 1}^{N_{x}} \sum_{n_{y} = 1}^{N_{y}} {(φ (ψ (n_{x}, n_{y})) - φ (ψ_{n} (n_{x}, n_{y})))}^{2}]}^{1 / 2},$ ${RMS}_{amplitude}^{n} = {[\frac{1}{N_{x} N_{y}} \sum_{n_{x} = 1}^{N_{x}} \sum_{n_{y} = 1}^{N_{y}} {(| ψ (n_{x}, n_{y}) | - | ψ_{n} (n_{x}, n_{y}) |)}^{2}]}^{1 / 2} .$
For a fair comparison, only the region inside the box support is taken into account for the calculations. The shearlet & range constraint provides the smallest error for both the phase and the amplitude after 100 iterations. For the range constraint and the support & range constraint, the RMS error of the amplitude does not decrease but increases with the number of iterations, probably due to the high noise level leading to an inconsistent feasibility problem. With the shearlet & range constraint, the algorithm seems to be more noise tolerant. As unstructured noise is associated with small shearlet coefficients, the shearlet thresholding always implies a denoising of the reconstruction, as visible in Fig. 2(f).
The simulation demonstrates that the constraint proposed in [27] can be successfully extended to the physical situation with a complex-valued exit wave field behind a mixed object. Moreover, in this setting the shearlet & range constraint distinctly outperforms the range constraint and the support & range constraint.
4. Experimental data
Figure 3 shows the analysis of an experimental data set. X-ray propagation images were recorded for a Siemens star on a high-resolution test pattern (ATN/XRESO-50HC, NTT Advanced Technology) consisting of radial stripes in the form of isosceles trapezoids of different size, arranged in rings around a center. The trapezoidal small bases of each ring have a defined length as indicated by numbers in the pattern. They reach from 4 µm down to 50 nm in the center. The pattern is defined in a 200 nm thick layer of tantalum on a Ru/SiC/SiN membrane. Due to the small thickness, the specimen’s absorption is negligible at the given photon energy of E = 8 keV, i.e. we can treat the pattern as a pure phase shifting object. As the star is etched into the tantalum layer, the phase shifts of the structure are positive relative to the surroundings (i.e. here the range constraint reads |ψ| = 1, φ(ψ) ≥ 0).

Fig. 3 Experimental data. (a) Reconstructed phases from experiment 1 obtained with 100 iterations of the RAAR algorithm and different object plane constraints. Reconstructed is only a 1061 × 1061 pixel central detail of the 2048 × 2048 pixel near field hologram. The inlays show the details of the reconstructions as indicated by the dashed white boxes. (b) PSDs of the three reconstructions in (a). The rings denote the location of spatial frequency components corresponding to a half-period resolution of 250 nm and 350 nm as indicated by the legend. (c) Histograms of the phase values from the three reconstructions in (a). For comparison we plotted the same range of phase values for all three reconstructions, however, some runaway values are not captured by this region. Note the different y-scales. The indicated maxima were determined by Gaussian fitting. (d) Reconstructed phases from experiment 2 obtained with 100 iterations of the RAAR algorithm and different object plane constraints. Depicted are the central 700 × 700 pixel details of the reconstructions. The larger magnification factor in experiment 2 compared to experiment 1 provides a higher resolution.
Download Full Size | PDF
The Siemens star was imaged at beamline | P10 of the PETRAIII storage ring (DESY, Hamburg) at instrument settings given in [38]. For detection, a 2048 × 2048 pixel SCMOS camera with a gadolinium oxysulfide scintillator of 15 µm thickness was used. A divergent x-ray waveguide probe was used in order to reach nanoscale resolution. The effective near field parameters of the equivalent parallel beam geometry are computed from the variable transformation given by the Fresnel Scaling theorem [10]. Two measurements with different waveguide-specimen and specimen-detector distances, i.e. with different magnification factors, were performed. Table 2 summarizes the relevant setup parameters for both experiments. A detail (1061 × 1061 pixel) of the recorded near field hologram from experiment 1 after applying the dark field and flat field correction is shown in the experimental setup in Fig. 1(a). For both experiments the initial guess for the phase retrieval iteration is composed of a wave field with amplitude equal to one and the phase of the holographic reconstruction.
Table 2. Relevant setup parameters of the two near field diffraction experiments with the Siemens star (Section 4), including the effective near field parameters of the parallel beam geometry.
View Table | View all tables in this article
Figure 3(a) shows the reconstruction results (1061 × 1061 pixel detail) for experiment 1 after 100 iterations of the RAAR algorithm (relaxation parameter α = 0.55) with different object plane constraints. Since the test pattern covers the whole field of view, the support constraint is not applicable. For constraints based on shearlet thresholding a shearlet system of 7 scales and a threshold parameter of θ_p = 0.001 is used. Compared to the reconstruction with the range constraint, the reconstructions based on shearlet thresholding show less noise and appear to exhibit a higher resolution. The rings corresponding to 4 µm down to 500 nm spatial half-period at the innermost radius are distinguishable in all three reconstructions. However, the stripes of the next smaller ring (200 nm half-period at innermost radius) are merely identifiable in the range constraint reconstruction due to the high noise level, see inlays in Fig. 3(a). The effective pixel size of about 120 nm impedes to resolve the stripes of the two smallest rings (with spatial half-periods of 100 nm and 50 nm) in all reconstructions. The higher performance in terms of resolution of the constraints based on shearlet thresholding compared to the simple range constraint is further reflected in the power spectral densities (PSDs) of the reconstructed phases, see Fig. 3(b). For the range constraint, structures can be distinguished up to spatial frequencies corresponding to a 350 nm half-period resolution before information is erased by noise. Contrarily, for the constraints based on shearlet thresholding the structures reach up to the 250 nm half-period resolution ring and, in parts, even beyond.
Figure 3(c) shows the histograms of phase values corresponding to the three reconstructions in Fig. 3(a). As the Siemens star consists of only two different thickness levels along the direction of illumination, we expect a histogram containing two sharp peaks. In the histogram of the range constraint we observe one very high peak at zero rad and two broader peaks corresponding to the two different object thicknesses. These peaks are smeared out due to noise and they would reach into the negative phase region. However, as the range constraint forces the phase to be nonnegative, the high peak at zero rad is formed. In the histogram of the shearlet & range constraint we have the same phenomenon, but the two relevant peaks are sharper and separated more clearly compared to those from the range constraint. The histogram corresponding to the shearlet constraint exhibits two very broad peaks that melt together rather strongly. Altogether, the shearlet & range constraint yields a histogram in closest resemblance to a binary object. This also reflects the high image contrast present in the corresponding reconstruction. Regarding the quantitativeness of the reconstructions, we observe the following phase shift differences between the relevant maxima (determined by Gaussian fitting) in the histograms of the range constraint and the shearlet & range constraint: Δφ_range ≈ 0.31 rad and Δφ_{shearlet&range} ≈ 0.32 rad. These values agree very well with the value Δφ ≈ 0.33 rad expected from x-ray optical constants (www.cxro.lbl.gov), i.e. quantitative phase retrieval is provided in both cases.
The RAAR algorithm is also used for the phase retrieval in experiment 2. Here, a shearlet system with 8 scales (the whole 2048 × 2048 pixel camera field is reconstructed) and a threshold parameter of θ_p = 0.001 are employed. Figure 3(d) shows central details (700 × 700 pixel) of the reconstructed phases obtained with different object plane constraints. As in experiment 1, we observe very clean reconstructions, but at higher resolution (stronger magnification). Now, we can even distinguish the trapezoids of the ring with 100 nm spatial half-period. By-eye-inspection the shearlet & range constraint provides again the best result.
These experiments demonstrate that the constraints based on shearlet thresholding are robust against noise and imperfections of real data. Due to their denoising property they can improve the resolution of the reconstruction and provide results superior to those from the simple range constraint.
5. Conclusion
In this paper we have analyzed the performance of the shearlet constraint for iterative phase retrieval in x-ray propagation imaging, both for simulations and for experimental data, and have compared reconstruction quality to other commonly used constraints. The shearlet constraint as proposed in [27] uses the knowledge that the exit wave field can be sparsely represented in a suitable shearlet frame. We have shown that the original formulation can be extended to the physical situation of a specimen that causes both phase shifts and absorption of the illuminating probe (complex-valued exit wave), providing reasonable results for the challenging setting of reconstructing from a single distance exposure (Section 3). Next, the approach was tested for experimental data (Section 4). The reconstructions were found to be consistent with conventional constraints, and partly outperform these. In particular, the shearlet thresholding implies a denoising effect which can lead to better image quality and resolution of the reconstruction.
The shearlet constraint can be applied without further information about the specimen needed. However, the used threshold parameter crucially influences the reconstruction. The optimal value depends not only on the object to be reconstructed but also on the measurement, i.e. the level of noise and the Fresnel number. Within this paper, the selected parameters were found empirically.
In order to obtain a deeper understanding regarding the optimum choice of the threshold and the constraint’s performance in general, it is necessary to further analyze the shearlet constraint within the RAAR algorithm. Further future work would be aimed at extending the shearlet constraint to phase retrieval in three-dimensional x-ray tomography and to far field CDI. The results presented in this paper suggest that constraints based on sparsity information can be advantageously used for phase retrieval tasks.
Acknowledgments
We thank Martin Krenkel and Anna-Lena Robisch for numerical tools and related fruitful collaborations, and acknowledge financial support by SFB 755 Nanoscale Photonic Imaging.
References and links
1. 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(12), 5486–5492 (1995). [CrossRef]
2. 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(1), 133–146 (1996). [CrossRef]
3. 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]
4. 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(14), 2961–2964 (1996). [CrossRef]
5. A. Pogany, D. Gao, and S. W. Wilkins, “Contrast and resolution in imaging with a microfocus x-ray source,” Rev. Sci. Instrum. 68(7), 2774–2782 (1997). [CrossRef]
6. M. Bartels, M. Krenkel, J. Haber, R. N. Wilke, and T. Salditt, “X-ray holographic imaging of hydrated biological cells in solution,” Phys. Rev. Lett. 114, 048103 (2015). [CrossRef]
7. S. Maretzke, “A uniqueness result for propagation-based phase contrast imaging from a single measurement,” Inverse Probl. 31, 065003 (2015). [CrossRef]
8. J. Als-Nielsen and D. McMorrow, Elements of Modern X-ray Physics (John Wiley & Sons, 2011). [CrossRef]
9. R. A. Lewis, “Medical phase contrast x-ray imaging: current status and future prospects,” Phys. Med. Biol. 49(16), 3573–3583 (2004). [CrossRef]
10. D. M. Paganin, Coherent X-Ray Optics (Oxford University, 2006). [CrossRef]
11. M. R. Teague, “Deterministic phase retrieval: a Green’s function solution,” J. Opt. Soc. Am. 73(11), 1434–1441 (1983). [CrossRef]
12. P. Cloetens, W. Ludwig, J. Baruchel, D. Van Dyck, J. Van Landuyt, J. P. Guigay, and M. Schlenker, “Holoto-mography: quantitative phase tomography with micrometer resolution using hard synchrotron radiation x rays,” Appl. Phys. Lett. 75(19), 2912–2914 (1999). [CrossRef]
13. J. P. Guigay, M. Langer, R. Boistel, and P. Cloetens, “Mixed transfer function and transport of intensity approach for phase retrieval in the Fresnel region,” Opt. Lett. 32(12), 1617–1619 (2007). [CrossRef] [PubMed]
14. H. N. Chapman, A. Barty, S. Marchesini, A. Noy, S. P. Hau-Riege, C. Cui, M. R. Howells, R. Rosen, H. He, J. C. H. Spence, U. Weierstall, T. Beetz, C. Jacobsen, and D. Shapiro, “High-resolution ab initio three-dimensional x-ray diffraction microscopy,” J. Opt. Soc. Am. A 23(5), 1179–1200 (2006). [CrossRef]
15. J. Miao, P. Charalambous, J. Kirz, and D. Sayre, “Extending the methodology of x-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens,” Nature 400, 342–344 (1999). [CrossRef]
16. S. Marchesini, “Invited article: a unified evaluation of iterative projection algorithms for phase retrieval,” Rev. Sci. Instrum. 78, 011301 (2007). [CrossRef]
17. V. Elser, “Phase retrieval by iterated projections,” J. Opt. Soc. Am. A 20(1), 40–55 (2003). [CrossRef]
18. T. E. Gureyev, “Composite techniques for phase retrieval in the Fresnel region,” Opt. Commun. 220(1–3), 49–58 (2003). [CrossRef]
19. K. Giewekemeyer, S. P. Krüger, S. Kalbfleisch, M. Bartels, C. Beta, and T. Salditt, “X-ray propagation microscopy of biological cells using waveguides as a quasipoint source,” Phys. Rev. A 83, 023804 (2011). [CrossRef]
20. A.-L. Robisch, K. Kröger, A. Rack, and T. Salditt, “Near-field ptychography using lateral and longitudinal shifts,” New J. Phys. 17, 073033 (2015). [CrossRef]
21. V. Davidoiu, B. Sixou, M. Langer, and F. Peyrin, “Nonlinear approaches for the single-distance phase retrieval problem involving regularizations with sparsity constraints,” Appl. Opt. 52(17), 3977–3986 (2013). [CrossRef] [PubMed]
22. R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik 35(2), 237–246 (1972).
23. J. R. Fienup, “Reconstruction of an object from the modulus of its Fourier transform,” Opt. Lett. 3(1), 27–29 (1978). [CrossRef] [PubMed]
24. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21(15), 2758–2769 (1982). [CrossRef] [PubMed]
25. H. H. Bauschke, P. L. Combettes, and D. R. Luke, “Hybrid projection-reflection method for phase retrieval,” J. Opt. Soc. Am. A 20(6), 1025–1034 (2003). [CrossRef]
26. D. R. Luke, “Relaxed averaged alternating reflections for diffraction imaging,” Inverse Probl. 21(1), 37–50 (2005). [CrossRef]
27. S. Loock and G. Plonka, “Phase retrieval for Fresnel measurements using a shearlet sparsity constraint,” Inverse Probl. 30, 055005 (2014). [CrossRef]
28. D. Labate, W.-Q. Lim, G. Kutyniok, and G. Weiss, “Sparse multidimensional representation using shearlets,” Proc. SPIE 5914, 254–262 (2005).
29. K. Guo, G. Kutyniok, and D. Labate, “Sparse multidimensional representations using anisotropic dilation and shear operators,” in Proceedings of the International Conference on the Interactions between Wavelets and Splines, G. Chen and M. Lai, eds. (Nashboro, 2006), pp. 189–201.
30. G. Kutyniok and D. Labate, Shearlets: Multiscale Analysis for Multivariate Data (Birkhäuser, 2012). [CrossRef]
31. P. Kittipoom, G. Kutyniok, and W.-Q. Lim, “Construction of compactly supported shearlet frames,” Constr. Approx. 35, (1)21–72 (2012). [CrossRef]
32. D. L. Donoho, “Sparse components of images and optimal atomic decomposition,” Constr. Approx. 17(3), 353–382 (2001). [CrossRef]
33. K. Guo and D. Labate, “Optimally sparse multidimensional representation using shearlets,” SIAM J. Math. Anal. 39(1), 298–318 (2007). [CrossRef]
34. G. Kutyniok and W.-Q. Lim, “Compactly supported shearlets are optimally sparse,” J. Approx. Theory 163(11), 1564–1589 (2011). [CrossRef]
35. W.-Q. Lim, “Nonseparable shearlet transform,” IEEE Trans. Image Process. 22(5), 2056–2065 (2013). [CrossRef] [PubMed]
36. G. Kutyniok, W.-Q. Lim, and R. Reisenhofer, “Shearlab 3D: faithful digital shearlet transforms based on compactly supported shearlets,” ACM Trans. Math. Software 42(1) 100 (2015).
37. D. L. Donoho, “De-noising by soft-thresholding,” IEEE Trans. Inform. Theory 41(3), 613–627 (1995). [CrossRef]
38. T. Salditt, M. Osterhoff, M. Krenkel, R. N. Wilke, M. Priebe, M. Bartels, S. Kalbfleisch, and M. Sprung, “Compound focusing mirror and x-ray waveguide optics for coherent imaging and nano-diffraction,” J. Synchrotron Radiat. 22(4), 867–878 (2015). [CrossRef] [PubMed]

constraint	effect for a… mixed object	pure phase object	pure absorption object
range constraint	\|ψ\| ≤ 1, φ(ψ) ≤ 0	\|ψ\| = 1, φ(ψ) ≤ 0	\|ψ\| ≤ 1, φ(ψ) = 0
support & range constraint	\|ψ(i, j)\| = 1, φ(ψ(i, j)) = 0 for (i, j) ∉ D \|ψ(i, j)\| ≤ 1, φ(ψ(i, j)) ≤ 0 for (i, j) D	\|ψ(i, j)\| = 1 for (i, j) ∈ Ω φ(ψ(i, j)) = 0 for (i, j) ∉ D φ(ψ(i, j)) ≤ 0 for (i, j) ∈ D	\|ψ(i, j)\| = 1 for (i, j) ∉ D \|ψ(i, j)\| ≤ 1 for (i, j) ∈ D φ(ψ(i, j)) = 0 for (i, j) ∈ Ω
shearlet constraint	soft thresholding for \|ψ\| and φ(ψ) with threshold θ_a and θ_p	soft thresholding for φ(ψ) with threshold θ_p and \|ψ\| = 1	soft thresholding for \|ψ\| with threshold θ_a and φ(ψ) = 0
shearlet & range constraint	soft thresholding for \|ψ\| and φ(ψ) with threshold θ_a and θ_p and \|ψ\| < 1, φ(ψ) ≤ 0	soft thresholding for φ(ψ) with threshold θ_p and \|ψ\| = 1, φ(ψ) ≤ 0	soft thresholding for \|ψ\| with threshold θ_a and \|ψ\| ≤ 1, φ(ψ) = 0

setup parameter	experiment 1	experiment 2
photon energy E	8 keV
detector pixel size Δ	6.54 µm
exposure time	1 s
waveguide-specimen distance	0.094 m	0.024 m
specimen-detector distance	5.034 m	5.104 m
magnification factor M	54.61	214.48
effective distance specimen-detector z_eff	0.092 m	0.024 m
effective pixel size Δ_eff	0.120 µm	0.030 µm
Fresnel number F	1.004 · 10⁻³	2.521 · 10⁻⁴

constraint	effect for a… mixed object	pure phase object	pure absorption object
range constraint	\|ψ\| ≤ 1, φ(ψ) ≤ 0	\|ψ\| = 1, φ(ψ) ≤ 0	\|ψ\| ≤ 1, φ(ψ) = 0
support & range constraint	\|ψ(i, j)\| = 1, φ(ψ(i, j)) = 0 for (i, j) ∉ D \|ψ(i, j)\| ≤ 1, φ(ψ(i, j)) ≤ 0 for (i, j) D	\|ψ(i, j)\| = 1 for (i, j) ∈ Ω φ(ψ(i, j)) = 0 for (i, j) ∉ D φ(ψ(i, j)) ≤ 0 for (i, j) ∈ D	\|ψ(i, j)\| = 1 for (i, j) ∉ D \|ψ(i, j)\| ≤ 1 for (i, j) ∈ D φ(ψ(i, j)) = 0 for (i, j) ∈ Ω
shearlet constraint	soft thresholding for \|ψ\| and φ(ψ) with threshold θ_a and θ_p	soft thresholding for φ(ψ) with threshold θ_p and \|ψ\| = 1	soft thresholding for \|ψ\| with threshold θ_a and φ(ψ) = 0
shearlet & range constraint	soft thresholding for \|ψ\| and φ(ψ) with threshold θ_a and θ_p and \|ψ\| < 1, φ(ψ) ≤ 0	soft thresholding for φ(ψ) with threshold θ_p and \|ψ\| = 1, φ(ψ) ≤ 0	soft thresholding for \|ψ\| with threshold θ_a and \|ψ\| ≤ 1, φ(ψ) = 0

setup parameter	experiment 1	experiment 2
photon energy E	8 keV
detector pixel size Δ	6.54 µm
exposure time	1 s
waveguide-specimen distance	0.094 m	0.024 m
specimen-detector distance	5.034 m	5.104 m
magnification factor M	54.61	214.48
effective distance specimen-detector z_eff	0.092 m	0.024 m
effective pixel size Δ_eff	0.120 µm	0.030 µm
Fresnel number F	1.004 · 10⁻³	2.521 · 10⁻⁴

Using sparsity information for iterative phase retrieval in x-ray propagation imaging

Abstract

1. Introduction

2. The shearlet constraint* for iterative phase retrieval*

3. Simulation

4. Experimental data

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (3)

Tables (2)

Equations (10)

Optics Express

Abstract

1. Introduction

2. The shearlet constraint for iterative phase retrieval

3. Simulation

4. Experimental data

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (3)

Tables (2)

Equations (10)

Optics Express

2. The shearlet constraint* for iterative phase retrieval*