1. Introduction

Lensless coherent diffractive x-ray imaging (CDI) has opened up a new field of high resolution structure analysis beyond the ensemble averaging of conventional x-ray diffraction [1–3]. Typically, lensless x-ray imaging setups are closer to diffraction experiments than to a classical microscope setup with lenses, except that they require a sufficiently coherent probing wavefront or beam as well as sufficient sampling of the diffraction pattern. In other words, the imaging systems is essentially based on free space propagation between object and detector. Depending on whether the data is recorded in the optical near-field or far-field, the propagation is modeled by the Fresnel propagator or a Fourier transform (Fraunhofer far-field regime), respectively. As in conventional diffraction, each detector pixel carries information about all object pixels. Therefore, if the data is modeled in the detector plane, high spatial resolution can only be achieved in terms of the average structure. Contrarily, if the data is inverted to reconstruct the object, the individual real space configuration is depicted, beyond an ensemble average (or more precisely an average over the entire illuminated volume). Importantly, the image formation is still very similar to a plain diffraction experiment, while the data analysis is not. Data Modeling for example by least-square fitting is replaced by image reconstruction. Model formulation is replaced by the formulation of a priori knowledge (constraints), which are required to compensate for the missing information on the phase of the diffraction field, and hence to achieve a unique solution. Accordingly, the difference is brought about by the inversion of the diffraction process: by solving the inverse problem of (non-crystallographic) diffraction, we obtain access to the individual configuration instead of the average sizes and correlation lengths in the object. It is for this reason that iterative algorithms [4–8], being to date the engine of CDI, have received so much attention. Iteratively cycling between the detector and object planes, they feed in both measured data and additional a priori information on the solution.

CDI uses a priori knowledge, for example related to the object support or its optical constants (positivity, pure phase contrast). Diversity in the data generated by illuminating the same object pixels by different wavefronts may be exploited by ptychographic algorithms [9–13]. In general, however, one is interested in finding experimental settings allowing for robust phase retrieval with least invasive constraints, dose and accumulation times. In many applications for example, one is interested in reconstructing from single recordings, without scanning or multiple exposures. At the same time, the object may be composed of several materials, both with phase and amplitude contrast. It is therefore often of advantage to perform the imaging in the optical near-field rather than far-field. Here, phase information is much more directly encoded in the diffraction pattern, by ways of interference between diffracted and reference waves. Recently, improved uniqueness results have been presented for near-field imaging using the theory of entire functions [14]. In fact, contrary to a common belief, see e.g. [15, 16], measurements at only one distance are sufficient to determine both the real and the imaginary part of any compactly supported refractive index distribution. These results are much stronger than those for the Fraunhofer regime, as obtained e.g. by [17]. Loosly speaking, the near-field case can thus be expected to require less information to complement the measurements than the far-field case. In practice, however, many iterative algorithms which have been adapted to the near-field case [18–20] require similar data diversity or constraints.

Notwithstanding the merits of many highly performing iterative algorithms, it is therefore necessary to broaden the perspective of the phase retrieval field. Here, we present iteratively regularized Gauss-Newton methods (IRGNM) [21] as an alternative approach to phase retrieval and other imaging problems. In this method, each iterate is computed to provide an optimal compromise between agreement with the measured data and additional constraints on the basis of a local linearization of contrast formation, as we discuss further below. The approach is related to the regularized gradient descent methods for phase reconstruction proposed in [22], but promises improved convergence owing to the Newton-like solution of linearized subproblems. We apply the general IRGNM framework to near-field phase contrast with x-rays, for the reasons mentioned above, and since recent progress in propagation imaging has narrowed the gap in resolution compared to CDI. In fact, using highly divergent and coherent quasi-spherical wavefronts, x-ray imaging in the optical near-field has been recently demonstrated down to 20 nm resolution [23].

The scope of the present work is two-fold. Firstly, we give a concise overview of iteratively regularized Gauss-Newton methods in view of x-ray imaging, summarizing and explaining recent mathematical literature in this field for an applied audience. Secondly, we demonstrate the performance of this approach in solving the phase problem on the level of real state-of-the-art experimental data. In particular, we apply an IRGNM approach to three-dimensional (3d) imaging, i.e. we show how the method can be used to perform phase retrieval and tomographic reconstruction simultaneously. This strategy has been argued to enable improved accuracy compared to splitting the reconstruction into phase retrieval problems for each angle to recover the fields in the object plane and a subsequent inversion of the Radon transform [24–27]. We show that IRGNM approaches offer significant flexibility in treating various experimental setups and different a priori information. In the long run, we also expect advantages owing to the fact that the noise characteristics of the data can be suitably accounted for in this framework. For the most important example of Poisson data, Newton-type regularization methods with Kullback-Leibler-type data fidelity terms have already been proposed [28].

The manuscript is organized as follows: §2 introduces IRGNM in view of image reconstructions problems from the principal idea to practical implementation. In §3, the framework is applied to (2d) near-field phase contrast with hard x-rays, imaging the phase shifts and absorption induced by a nano-structured test pattern. A Kaczmarz-type IRGNM suitable for tomographic imaging is presented in §4 and applied to resolve the structure of a colloidal micro-crystal.

2. Regularized Newton methods for imaging

2.1. Basic approach

We consider an abstract imaging system of the form

Here I^obs ∈ Y denotes some observable intensity data, given by the image of the unknown object f^† ∈ X (e.g. a spatially varying refractive index or scatterer positions) under a known forward operator F and superimposed measurement errors ε. Note that ε may depend on F(f^†) as is the case for Poisson data. For many models, an explicit inversion formula for F is not known. Moreover, even if the inverse F⁻¹ is available, it is often not continuous and would thus amplify errors ε by large magnitudes if applied directly to the data I^obs. Such ill-posedness of the reconstruction despite uniqueness is well-established e.g. for computed tomography [29]. We seek a method to stably recover f^† from an ill-posed problem of the form (1).

Notably, the operator F : X → Y modeling the imaging system is nonlinear in general, for example whenever a phase retrieval problem described by a squared modulus operation is involved. Nevertheless, reasonable results can often be achieved using a linearization of contrast formation, as given for instance by the contrast transfer function in electron microscopy or coherent x-ray imaging [30–33].

Mathematically, such first order approximations are justified by the Fréchet differentiability of the forward operator F, i.e. for any f there exists a bounded linear map F′[f] such that lim_h→0(F(f + h) − F(f) − F′[f]h)/∥h∥ = 0. A natural approach to solve (1) is then given by Newton-type iterations

As opposed to methods based on a static linearization of contrast formation, the linearizations in (2) are computed about the current iterate f_k and thereby account for moderate nonlinearity.

However, iterations of the form (2) are often neither feasible nor desirable for imaging because the linearized problems - just like the nonlinear equation (1) - are typically ill-posed. Hence, the solution of (2) is unstable and may not even exist. Again, one may think of ambiguities in phase retrieval problems. A remedy is given by iteratively regularized Gauss-Newton methods (IRGNM) as first proposed by [21], corresponding to Tikhonov regularization of the Newton steps:

Here, ∥·∥_Y and ∥·∥_X denote the norms in Hilbert spaces X and Y, f₀ ∈ X is the initial guess and α_k > 0 is a regularization parameter. In this setting, it can be shown that (3) always has a unique solution given by [34]

F′[f_k]^* denotes the adjoint of the linear map F′[f_k] : X → Y. In (4), only the inverse of the selfadjoint positive-definite operator F′[f_k]^*F′[f_k] + α_k has to be computed, which is bounded according to the estimate

Accordingly, the iterate f_k+1 depends continuously on I^obs, i.e. the impact of data errors on the reconstruction is regularized.

As the IRGNM is based on linearizations of the imaging operator, yet iteratively updated, the approach is best suited for weakly or moderately nonlinear problems. Formally, convergence of the method to f^† for ε → 0 can indeed be shown given bounds on the nonlinearity of F along with suitably chosen α_k and f₀ [35, 36].

2.2. Parameter choice and constraints

The first term on the right hand side of (3) measures the agreement of the object with the observed data I^obs based on the current linearization. To achieve competitive results, the choice of the norm ∥·∥_Y should reflect the statistical properties of the data errors ε, e.g. by taking the negative log-likelihood of the measured signal I^obs. For additive Gaussian white noise, this consideration leads to the choice of the standard L²-norm, i.e.

For Poisson noise, the resulting data fidelity term is the Kullback-Leibler-divergence. This distance measure can be implemented in the framework of generalized Newton methods as demonstrated by [28]. Within the IRGNM, a quadratic approximation about its minimum may be used as a norm

On the right hand side of (3), the data residual is balanced with the regularization term ${\alpha }_k{\Vert f-f_0\Vert }_X^2$ bounding the deviation from the initial guess f₀. The regularization parameter weights the different contributions: if α_k is very small, there is essentially no regularization and the norm bound (5) diverges, allowing for large amplifications of the data error. If α_k is chosen too large on the other hand, the Newton iterate computed via (4) need not have much to do with the actual image to be reconstructed. A good choice of α_k must thus balance data- and approximation errors. One possible strategy is the following:

The stopping criterion, known as Morozov’s discrepancy principle [37], requires good knowledge of the magnitude of data errors. When this is not available, noise level-free stopping rules need to be applied, see e.g. [34, 38, 39]. Here, the principal idea is to make s_k as small as possible while limiting heuristic measures for the impact of data errors such as the object norm ∥f_k − f₀∥_X or the inverse regularization parameter $\frac{1}{{\alpha }_k}$.

By the choice of the norm ∥·∥_X, we may define desirable properties of the object f to be reconstructed. Choosing the standard L²-norm here prevents isolated spikes, promoting more evenly distributed values. Typical images are expected to be of higher regularity, for example being composed of smoothly varying areas bounded by sharp edges. This may be exploited to obtain higher robustness to high-frequency errors by regularizing with Sobolev norms

Beyond smoothness it is often desirable to impose additional constraints, which corresponds to restricting the set of admissible solutions f_k ∈ C ⊂ X. Prominent examples are real-valuedness, support constraints or positivity. Geometrically, the former two types are represented by linear subspaces C ⊂ X. Imposing these constraints in the IRGNM simply amounts to substituting F′[f_k]^* with $\mathrm{P}F\prime {[f_k]}^{\ast }$ in (4), where P : X → C is the orthogonal projection onto C. Positivity, on the other hand, is a nonlinear convex constraint. It may be included within a generalized Newton framework via a nonsmooth regularization term. In practice, this amounts to solving the minimization problem (3) restricted to f ∈ C, as can be done using semismooth Newton methods [40]. One approach to approximate sign constraints within the IRGNM framework of the present work lies in supplementing (3) with the penalty term

For numerical implementation of the IRGNM, all that needs to be done is to exchange the imaging operator F and its derivative F′ in (3) by suitable discrete approximations. The norms ∥·∥_Xdis and ∥·∥_Ydis in the discretized object- and image spaces $X_{\mathrm{dis}}={ℝ}^{N_X}$, $Y_{\mathrm{dis}}={ℝ}^{N_Y}$are characterized by their Gramians G_X,G_Y with respect to the Euclidean norm, i.e.

The Hermitean positive-definite linear problem in the Newton step (3) can be solved efficiently by a conjugate gradients (CG) method. A major advantage of this IRGNM-CG scheme is that only the forward maps F_dis, ${F^{\prime }}_{\mathrm{dis}}[f]$ and ${F^{\prime }}_{\mathrm{dis}}{[f]}^{\mathrm{T}}$ need to be evaluated, which are often easy to implement. Moreover, note that the matrices ${F^{\prime }}_{\mathrm{dis}}[f]$ usually need not (and should not) be set up explicitly to compute the matrix-vector products $h\to {F^{\prime }}_{\mathrm{dis}}{[f]}^{(}{}^{\mathrm{T})}h$.

2.3. Newton-Kaczmarz methods

The Newton-type method presented so far is an all-at-once approach to image reconstruction: linearizations of the problem (1) are solved in each iteration, incorporating all constraints and the complete measured data I^obs. In view of 3d or even 4d tomographic or time-resolved imaging data, this results in computations with huge arrays, posing numerical challenges. For linear problems, a well-known remedy is given by Kazcmarz-type methods [41] such as the Algebraic Reconstruction technique in computed tomography [29, 42], which cyclically solve small under-determined subproblems. Indeed, many nonlinear imaging problems also lend themselves to a natural separation into different subproblems so that (1) becomes

For problems of this form, [43] proposed regularized Newton-Kaczmarz methods, where in each iteration only one of the linearized subproblems $F_{jk}\left(f_k\right)+{F^{\prime }}_{jk}\left[f_k\right]\left(f_k{}_{+1}-f_k\right)=I_{jk}^{\mathrm{obs}}$ is considered. A Kaczmarz-type equivalent of the IRGNM-update (3) takes the form

The additional regularization term bounds the deviations from the preceding iterate f_k according to the weights β_k ∈ [0;1], ensuring that the reconstruction does not change too much within one Newton iteration. The sequence (j_k) determines the processing order of the different sub-problems which may be adjusted to the requirements of the particular imaging problem.

3. Application to propagation-based phase contrast

Next, we apply the regularized Newton framework presented in section 2 to propagation-based phase contrast x-ray imaging. An exemplary experimental setup is sketched in Fig. 1: quasi-monochromatic undulator radiation is focused by a pair of elliptical mirrors onto an x-ray waveguide. The exit of the waveguide serves as a quasi point source illuminating an object, which is placed at a distance of several mm behind the waveguide. The resulting diffraction patterns (holograms) are recorded by a detector at about 5m distance behind the focal plane, capturing the entire cone beam emanating from the waveguide with the sample induced interference pattern. The imaging data presented in this work has been recorded at the GINIX endstation at P10 beamline, DESY, Hamburg, described in [44, 45]. For details on the waveguide system and comparable high resolution propagation imaging results, see also [23].

 

Fig. 1 Schematic setup for propagation-based phase contrast imaging with hard x-rays: GINIX at P10 beamline, DESY [45, 48]. Quasi-monocromatic x-rays are focused onto a waveguide, illuminating a downstream object by a cone beam emanating from this coherent quasi-point source. Rotation of the sample allows for tomographic measurements.

Download Full Size | PDF

For the imaging model, a fully coherent illumination of the sample by a plane wave is assumed. Note that we consider an effective parallel beam geometry equivalent to the experimental cone-beam setup. As in previous studies of cone beam x-ray propagation imaging [46, 47], we make use of the simple variable transformation, accounting for the geometric magnification, mapping the spherical beam illumination to an effective plane wave case. For the present imaging system, the data treatment is detailed by [23]. Moreover, we assume the object to be sufficiently thin and weak for the wave field to remain spherical and for the transmission to be well described by geometrical optics (projection approximation). This is typically well-satisfied for hard x-ray imaging [7, 47]. Under these assumptions, the intensity at the detector can be modeled as

From (14), we see that phase contrast imaging is nonlinear in general. Nevertheless, the model can be linearized for weak phase shifts and attenuation ϕ, µ, yielding the contrast transfer function (CTF) [32, 49]

Assuming vanishing attenuation µ ≈ 0 or single-material objects µ ∝ ϕ, this formula can be directly inverted for image reconstruction. This approach yields satisfactory reconstructions especially if data from multiple defocus distances is available [33, 48].

From the CTF solution (16), it can be seen that regularized Newton methods are well-suited for phase contrast imaging in two respects: firstly, the zeros of the oscillatory prefactors make the image reconstruction problem ill-posed even in the simplest case of a weak non-absorbing object so that regularization is required. Secondly, the relative success of linear CTF-based methods indicate that the nonlinearity of the imaging method is sufficiently weak in a large regime of experimental setups.

For suitable spaces X,Y, the imaging operator F : X → Y defined by (14) can be shown to be Fréchet differentiable. Using linearity of $\mathcal{D}$ and the expansions $\exp (x+h)=\exp (x)(1+h)+\mathcal{O}(h^2)$ and ${\left|x+h\right|}^2={\left|x\right|}^2+2\Re \left(\overline{x}h\right)+\mathcal{O}\left(h^2\right)$ we obtain the first-order Taylor approximation

This allows to read off the Fréchet derivative given by the linear term in h:

Here the overbar denotes complex conjugation and ℑ the pointwise imaginary part. The special case f = 0 in (18) reproduces the CTF. Implementing (14) and (18) for discrete images in the spaces $X:={ℂ}^{M_X\times N_X}$ and $Y:={ℝ}^{M_Y\times N_Y}$, using suitably padded fast Fourier transforms for the discrete Fresnel propagator, the reconstruction problem of propagation-based phase contrast imaging can be solved by the IRGNM of section 2.

As a proof of concept, we consider the diffraction pattern of an IRP logo, engraved into a thin gold film. The data has been recorded on the aforementioned GINIX setup at an energy of E = 7.9keV is shown in Fig. 2(a). The plotted intensity data has been normalized by the corresponding flat-field, i.e. the image of the empty beam, in order to meet the assumption of plane wave illumination as justified by the analysis of [50]. The Fresnel number is N_F = 1.77·10⁻⁴ at an effective pixel size of a = 21.7nm of the 1080 × 1920-sized images and a maximum flux per pixel of ≈ 3400 photons.

 

Fig. 2 X-ray phase contrast imaging of a nano-structured object using the IRGNM reconstruction algorithm. The test pattern (institute logo) was defined in a thin gold film by focused ion beam milling. (a) Diffraction pattern recorded at GINIX endstation, P10 beam-line, DESY [45] (flat-field corrected, maximum photon flux per pixel ≈3400). (b) Reconstructed phase image of the entire field of view, assuming a fixed ratio of 0.21 between absorption µ and phase shifts ϕ. (c) Magnification of (b) in the framed region around the logo. For comparison, the dashed inset shows the corresponding part of the phase map reconstructed by direct inversion of the CTF (16) via the methods of [48]. Scale bars: 2 µm in the effective geometry. Fringe artifacts on the upside of the logo are due to diffraction fringes leaving the field of view in (a). The phase shifts ϕ ≈ 0.20 induced by the object correspond to a thickness of the gold layer of about 100nm.

Download Full Size | PDF

For the Newton reconstruction, we choose the Sobolev norm in (8) with an exponent s = 0.5 as a regularization term and an L²-norm for the data residual. The regularization parameters α_k are chosen according to the procedure outlined in section 2.2. The algorithm is stopped as soon as the reduction of the residual s_k − s_k−1 achieved in the k-th iteration falls below 1% of the maximum decrease max_j<k s_j − s_j−1. For the initial guess, we simply take f₀ = 0. We exploit that the specimen is made of gold by prescribing its characteristic between absorption µ and phase shifts ϕ of c_µ/ϕ ≈ 0.21 [51]. The constraint is imposed by substituting $f=\left(1\frac{\mathrm{i}}{2}{c}_{\mu /\varphi }\right)\varphi$ and reconstructing the real-valued ϕ. Apart from this, neither positivity nor a support constraint is assumed. The latter means that the phase shifts ϕ are computed within the entire field of view of 1080×1920 pixels without any oversampling

Convergence of the IRGNM is reached after 12 Newton steps, corresponding to a total of 242 CG-iterations. The resulting phase map is plotted in Fig. 2(b). Despite the lack of oversampling, it can be seen that both the object and the background come out quite clean except for some stripes near the boundary due to the applied smooth replicate padding. Occasional white spots and the smooth background variations, on the other hand, may be attributed to some dirt in the imaging optics and physical variations of the thickness of the gold layer, respectively. The magnification of the reconstruction in a region around the IRP-logo depicted in Fig. 2(c) confirms the high uniformity of the recovered phase shifts in the regions within and without the logo, reproducing the binary test pattern. Slight fringe artifacts can be identified in the vertical direction, especially on the upper edge of the logo. These may be explained by incompleteness of the data. Indeed, it can be seen from Fig. 2(a) that parts of the diffraction pattern lie outside of the recorded field of view. The reconstruction of the missing object information is based on a priori constraints only, which naturally gives rise to artifacts.

It should be noted that the imaging setting considered so far still lies well within the regime of applicability of reconstruction method based on the CTF (16): for once, the phase shifts within the logo of ≈ 0.20rad, induced by the 100nm-thick gold layer of refractivity kδ ≈ 1.96 µm⁻¹ [51], are sufficiently weak for a global linearization of contrast formation to be reasonably accurate. Secondly, the assumption of a single material, i.e. of proportional fields ϕ and µ, allows for a direct inversion of (16). In the considered example, the latter approach indeed yields results of comparable quality as shown by the dashed inset in Fig. 2(c). If ϕ and µ are considered as independent variables and only a single diffraction pattern is available, directly inverting (16) is impossible. This suggests that the imaging problem is non-unique in this case. On the contrary, theoretical investigation has recently shown that a unique recovery is indeed feasible if µ and ϕ are compactly supported [14]. The mathematical reason is that compactness of the support translates into correlations in Fourier space, which can be used to disentangle ϕ and µ in (16). By its ability to incorporate a priori constraints, the IRGNM approach may thus enable reconstructions beyond the limitations of direct inversion methods.

For a first experimental verification of a simultaneous recovery of both phase and absorption from a single hologram, we repeat the IRGNM reconstruction using the same parameters but without imposing a coupling of ϕ and µ. Following the uniqueness analysis, we impose a loose circular support constraint around the logo. Note that this induces a strong oversampling in the data by a factor of ≈ 17, owing to the considerably reduced number of object pixels on which ϕ and µ have to be reconstructed. Moreover, negativity of µ and ϕ (the test pattern represents missing material!) is imposed in the IRGNM via the penalty term approach outlined in section 2.2. The recovered phase and absorption obtained after 10 Newton steps (199 CG-iterations) is plotted in Fig. 3. For comparison, dashed insets show the corresponding image parts obtained from a second reconstruction without negativity constraints.

 

Fig. 3 Simultaneous IRGNM reconstruction of phase shifts ϕ and absorption µ from the data in Fig. 2(a) without assuming a fixed ratio µ/ϕ. Negative values of µ and ϕ in-dicate missing material in the gold film. The circular support visible in the phase- and absorption maps and negativity of µ and ϕ has been imposed as a constraint. All other parameters are retained as in the computation of Fig. 2(b). The dashed inset in the ab-sorption image shows the reconstruction without negativity constraint in the IRGNM and demonstrates that simultaneous recovery tends to introduce low-frequency artifacts. These are effectively suppressed by exploiting physical a priori knowledge on the sign of phase shifts and absorption. Scale bar: 2 µm.

Download Full Size | PDF

Both the phase shifts ϕ in 3(a) and the attenuation coefficient µ plotted in 3(b) clearly represent the object shape, where the magnitudes roughly reproduce the material-specific ratio of µ/ϕ ≈ 0.21. The only visible artifacts are those attributed to the missing fringes, which have been observed previously. In the IRGNM-reconstruction without negativity constraints, additional low-frequency errors appear especially in the recovered attenuation, as can be seen from the central dark spot surrounded by a bright halo in the inset of Fig. 3(b). Simulations as well as preliminary analytical studies indicate that the susceptibility to these halo artifacts is characteristic of simultaneous retrieval of phase and absorption. As observed in the present reconstruction, however, they seem to be suppressed very effectively by imposing suitable sign constraints. Owing to the flexibility of exploiting physical a priori knowledge in the IRGNM, we thus indeed achieve an almost perfect simultaneous recovery of µ and ϕ from a single hologram - which has to date been considered impossible [16, 52].

Although proportionality of µ and ϕ is not imposed in the reconstruction in Fig. 3, the surprising quality might still be owing to some implicit preference of the IRGNM for homogeneous objects. We test this by repeating the reconstruction for a simulated object that is composed of two different materials. To this end, a binarized version of the recovered IRP logo with constant ϕ = 0.2 and µ = 0.04 is embedded into a purely phase-shifting disc with ϕ = 0.2 and µ = 0. All setup- and reconstruction parameters are chosen exactly as in Fig. 3. A single synthetic hologram (not shown) is simulated by mapping the exact object with the forward operator and superimposing Gaussian white noise with standard deviation σ = 0.02. The IRGNM reconstruction of ϕ and µ from this data is shown in Fig. 4 along with the exact solutions.

 

Fig. 4 Supplementary IRGNM reconstruction of a non-homogeneous object from a simulated near-field hologram with 2% Gaussian white noise. Constraints, reconstruction- and setup parameters are chosen as for the experimental data in Fig. 3. (a) and (b): Exact object composed of an absorbing logo structure with ϕ = 0.2, µ = 0.04 (gold) embedded into a non-absorbing disc that induces the same phase shifts ϕ = 0.2. (c) and (d): Recovered phase and absorption images. Except for a slight halo, features in ϕ and µ are correctly identified, revealing the hidden logo structure by the enabled absorption contrast.

Download Full Size | PDF

Good agreement with the exact object is found both in phase and absorption. The different object features are correctly attributed to the ϕ- and µ-components in Figs. 4(c) and (d) and no artifacts except for the known slight halo and fringe structures are visible. In particular, the absorption contrast enabled by the simultaneous recovery clearly reveals the logo structure despite the low signal-to-noise-ratio in µ. The logo is not represented in the phase image and thus likely invisible to any reconstruction method that assumes proportionality of µ and ϕ.

4. Phase contrast tomography of a colloidal crystal

By rotating the specimen within the incident beam in the imaging setup of Fig. 1, propagation-based phase contrast can be extended to a tomographic imaging method, capable of resolving three-dimensional variations of the complex refractive index n = 1−δ +iβ. Within the geometrical optics approximation of section 3, the imprinted phase- and absorption images ϕ and µ are proportional to the projection of δ and β, respectively, along the optical axis. For a tomographic incident angle θ, we have

Since ℛ is linear and bounded, F_tomo is Fréchet differentiable with a derivative of similar form as (18).

Implementing (20) in a regularized Newton-type framework corresponds to an all-at-once approach to phase contrast tomography, as phase- and tomographic reconstruction are not carried out in subsequent steps but simultaneously. In particular, this implies that tomographic correlations between images under different incident angles are incorporated already in the phase retrieval step. Similarly as in coherent diffractive imaging [24, 25], this has been shown to improve stability and accuracy of the reconstruction [26, 27]. Unfortunately, a numerical inversion of F_tomo by an IRGNM is not viable for high resolution data sets due to memory constraints and high computational costs of the evaluation of the Radon transform.

As a remedy, we propose to solve the reconstruction problem of phase contrast tomography by a regularized Newton-Kaczmarz methods as presented in section 2.3. The tomographic setting suggests a decomposition of the forward operator F_tomo = (F₁,…,F_d), corresponding to small sub-data sets of diffraction patterns from only very few incident angles. Each Newton step (13) then only requires evaluations of comparably small Radon transforms and processing of feasible chunks of data. At the same time, tomographic consistency is exploited by directly reconstructing a three-dimensional object instead of single projections. Furthermore, the F_j may be chosen such that the most strongly correlated images for only slightly differing incident angles are processed simultaneously.

We implement the Newton-Kaczmarz step (13) for the imaging operator (20), using the simple parameter choice α_k = α₀ and β_k = β₀ ≪ 1, where α₀ is estimated according to section 2.2. Firstly, this regularization imposes the current reconstruction as a strong prior for the next iteration. Secondly, it also yields a favorable condition number for the linear Newton steps so that CG-methods typically converge after 3–5 iterations.

For a proof of concept, we consider phase contrast images of a colloidal crystal composed of 415nm-sized polystyrene beads, imaged at an energy E = 7.9keV in the experimental setup of Fig. 1. The tomographic data set is composed of 249 diffraction patterns of size 1024 × 1024 acquired at a maximum photon flux of 770 per pixel under incident angles θ between0 and 173°. The raw data correction included an alignment of the tomographic projections to the common center of mass and iterative reprojections, following the approach of [53]. The photon flux corresponds to an upper bound of 110kGy for the absorbed dose over the total tomographic data acquisition. The flat-field corrected intensities are visualized in Fig. 5(a). In the effective geometry, the Fresnel number is N_F = 2.41·10⁻⁴ at an effective pixel size of 29.3nm.

 

Fig. 5 X-ray phase contrast tomography of a colloidal micro-crystal composed of 415nm-size polystyrene beads, reconstructed by a regularized Newton-Kaczmarz method. The maximum photon flux is 770 per pixel corresponding to an absorbed dose ≤ 110kGy.(a) Slice visualization of the tomographic data given by 1024 × 1024 flat-field corrected holograms measured under 249 incident angles between 0 and 173° at GINIX endstation, P10 beamline, DESY [45]. (b) Central slice of the reconstructed 256³ voxel volume, plotting the increment δ of the refractive index n = 1−δ +iβ which is proportional to electron density. A priori constraints were positivity δ ≥ 0 and vanishing absorption β = 0. The scale bar is 1 µm. (c) Fourier shell correlation (FSC) computed for reconstructions from complementary data sets of 125 and 124 incident angles (green curve). For comparison, the blue curve shows the result without positivity constraint and the red dashed one plots the 1/2-bit threshold curve [54]. The intersection between the green- and the red curve indicates an achieved resolution of 95nm. (d) Binary representation of the slice in (b), determined by deconvolving the reconstructed δ with the form factor of a 415nm-sized homogeneous sphere after Gaussian filtering.

Download Full Size | PDF

From the fringes in the data set in Fig. 5(a) as well as preliminary reconstructions, one can infer the object location in the center of the field of view. This is exploited by restricting the reconstruction to a central 256³ voxel cube, imposing a loose 3d-support constraint. Moreover, the hydrocarbon composition of the polystyrene spheres and the photon energy allow using β = 0 as a constraint (non-absorbing object), i.e. to reconstruct only the refractive increment δ. Unfortunately, a Sobolev-type regularization term (8) would bottleneck the proposed Newton-Kaczmarz-scheme as the required FFTs in the Gramian would constitute the only ℴ (N³ logN) operations in the Newton updates. We therefore recur to simple L²-regularization and L²-data fidelity terms, corresponding to identity Gramians G_X,G_Y in (10). However, positivity of δ is imposed by supplementing the Kaczmarz-Newton step (13) with the penalty term in (9). The weights of the penalty terms are set to γ = α₀ and β_k = 0.001. For the initial guess, we simply choose f₀ = 0. The forward operator F_tomo is decomposed such that angular “wedges” of six adjacent diffraction patterns are processed simultaneously in each Newton step. The number of iterations is determined such that each hologram is visited twice where the processing order of the wedges is random to minimize directional bias in the reconstructed object.

The reconstruction runs on a simple laptop with Intel i7 CPU and 8 Gigabytes of RAM, computing for about 20 minutes. Each of the Newton steps takes four CG-iterations, which means that the total solution - processing every incident angle twice - requires only 10 forward and adjoint operations of the full imaging operator F_tomo. Figure 5(b) shows the central slice of the reconstructed δ perpendicular to the tomographic axis. It can be seen that the binary character of the specimen and the spherical shapes of the beads are well-resolved, being clearly distinguishable both from one another and from the background. Notably, no artifacts caused by angular undersampling or due to the missing wedge of 6° can be identified.

We estimate the resolution by computing the Fourier shell correlation (FSC) [55] of two auxillary reconstructions from complementary sets of 125 and 124 incident angles, respectively. In Fig. 5(c), the resulting correlation is compared to the 1/2-bit threshold curve as proposed by [54]. For comparison, the FSC is also computed for a reconstruction without positivity constraint. The beneficial effect of the latter can be seen from the gap between the blue and green curves in Fig. 5(c), showing a significantly higher correlation especially in the resolution-critical part around the intersections with the threshold curve. The intersections indicate a resolution of about 95nm and 105nm (half of the corresponding Fourier wavelengths) for the reconstructions with and without positivity constraint, respectively. Notably, the local minima of the FSC curves at |ξ| ≈ 0.022nm⁻¹ are no artifacts but neatly coincide with the first order zero of the form factor of 415nm-sized spheres. This emphasizes the sensitivity of the imaging method to structural features.

According to the slice plotted in Fig. 5(b), the positions of the polystyrene spheres and thus the crystal structure could be determined in principle by visual inspection of the reconstructed δ. In order to determine the colloid locations numerically and independent of an operator, we smoothen the recovered object by a Gaussian of 95nm FWHM to reduce the impact of noise before deconvolving with the form factor of a homogeneous 415nm-sized sphere (with regularization around the zeros in Fourier space). This procedure results in Gaussian peaks centered at the positions of polystyrene beads. By computing the maxima of this field using quadratic interpolation about the peaks, we thus obtain a representation of the crystal structure.

By the simple deconvolution procedure outlined above, 448 colloid positions are identified. Convolving the obtained Dirac delta-array of bead locations with the ideal 415nm-sized sphere yields a binary representation of the colloidal crystal. A slice of this binarization corresponding to Fig. 5(b) is plotted in 5(d). Figure 6 shows a 3d-rendering of the determined colloidal crystal. The corresponding site coordinates of the colloids, i.e. the centers of the spherical beads, are provided in Data File 1. The imaged micro-crystal contains regions of (approximately) hexagonal close-packing as well as cubic and amorphous regions, induced by an interplay of bulk-and surface effects.

 

Fig. 6 3d-rendering of a colloidal crystal of polystyrene nano-beads reconstructed by Newton-based phase contrast tomography. The underlying binary object, a slice of which is shown in Fig. 5(d), has been determined by deconvolving the smoothened reconstruction of the refraction parameter δ with the form factor of a homogeneous sphere of diameter 415nm. The coordinates of the colloid sites are provided in Data File 1. Scale bar: 1 µm.

Download Full Size | PDF

5. Conclusions

In this work, we have presented iteratively regularized Gauss-Newton methods (IRGNM) as a generic approach to solve nonlinear ill-posed image reconstruction problems. The principal idea is to reconstruct an unknown object by iteratively inverting linearizations of a known imaging model, which describes contrast formation from object to observable data. In order to compensate for missing information, the iterates are computed to provide an optimal compromise between the measurements and additional a priori information on the unknown object. The IRGNM approach differs from well-known alternating-projection-type algorithms typically used in CDI in that it exploits differentiability and simultaneously processes constraints and observed data. This promises improved convergence.

By applying regularized Newton methods to near-field phase constrast imaging, both in 2d and in a tomographic 3d setting, we have demonstrated their flexibility in treating different experimental setups. The reconstruction of a nano-structured test pattern shows that IRGNM constitute a reasonable generalization of direct inversion methods based on a global linearization of contrast formation. Owing to its potential to incorporate moderate nonlinearity and a priori knowledge, e.g. on support or positivity, the approach permits faithful reconstructions beyond the scope of such direct methods. This has been demonstrated by the simultaneous recovery of magnitude and phase of the test object from a single diffraction pattern without assuming proportionality of phase shifts and absorption. The validity of the reconstruction is supported by results for a simulated object composed of two materials, which are correctly identified by the algorithm. However, numerical simulations also indicate that the stability of the simultaneous recovery deteriorates significantly with increasing Fresnel number. Thus, it is possibly only feasible for deeply holographic near field data as considered in this work.

A further benefit of the presented regularized Newton framework is that it requires only a forward model for the imaging setup and no explicit knowledge of (approximate) inverses. This enables IRGNM reconstructions also for complicated multi-staged imaging setups, possibly including the influence of various experimental parameters. In the considered example of phase contrast tomography this flexibility permits to directly recover a probed 3d-object from the complete tomogram of diffraction patterns instead of subsequently performing (2d) phase retrieval and tomographic backprojection. We have demonstrated the potential of this approach by imaging a colloidal crystal composed of 415nm-sized polystyrene beads at 95nm resolution, reconstructed by an efficient regularized Newton-Kaczmarz method.

X-ray imaging techniques currently reach out for novel frontiers, ranging from new contrast modalities in material science and dose reduction in biomedical analysis to in-operando studies and investigations of highest spatial and temporal resolution by FEL pulses. To enable these goals, the mathematical modeling underlying image reconstruction has to become increasingly accurate in order not to bottleneck the achievable contrast and resolution: apart from nonlinearity of image formation and Poissonian noise statistics, future imaging models may need to better account for partial coherence, non-uniform illumination, mechanical vibrations and movements, misalignment or geometrical aberrations, just to name a few. Owing to their flexibility, we are convinced that regularized Newton methods, as presented in this work, may greatly foster such developments and thus contribute to a bring lensless x-ray techniques as well as other cutting-edge imaging methods to the level of quantitative structural measurements.