Maximum likelihood reconstruction for grating-based X-ray microscopy

Andreas Wolf; Max Schuster; Veronika Ludwig; Gisela Anton; Stefan Funk

doi:10.1364/OE.380940

1. Introduction

The Talbot effect, first observed in 1836 [1], is a diffraction phenomenon which can lead to repeated self-images of periodic structures. Besides having found a wide applicability in the optical regime [2], the Talbot effect also provides a basis for grating-based interferometry and imaging with X-rays [3,4], electrons [5,6], neutrons [7,8], and molecules [9]. In X-ray Talbot interferometry (XTI), the Talbot self-image of an X-ray diffraction grating serves as a spatial reference pattern. Distortions of the X-ray wave front, such as those caused by a sample in the beam, lead to deformations of the said pattern [4]. In particular, the pattern is modified by refraction in the sample, thus enabling X-ray phase-contrast imaging. In comparison with conventional X-ray absorption imaging, phase-sensitive techniques offer a higher contrast for samples consisting of low-Z materials [10].

Due to its compatibility with spherical-wave illumination, XTI can be combined with full-field transmission X-ray microscopes [11–15]. It adds the possibility of quantitative phase imaging with spatial resolution at the micro- and nanoscale to these instruments. First, Takeda et al. [11] incorporated an X-ray Talbot interferometer to an imaging microscope. The interferometer consisted of two gratings and was sensitive to the gradient of the wave front. With increasing magnification of the microscope, this leads to a reduction of the sensitivity [16] which consequently has to be traded off against the achievable spatial resolution [15]. The problem of reduced sensitivity was solved by Yashiro et al. [12] through the usage of a single, highly magnified grating. It enables the retrieval of phase difference images instead of differential phase images, as it is outlined in Section 2 of this paper. The same group published detailed theoretical descriptions of the wave field propagation in such phase difference microscopes [13] and transferred them to laboratory X-ray tubes [14], following up on the initial work on X-ray Talbot-Lau interferometry by Pfeiffer et al. [17]. More recently, Takano et al. [18] installed a Talbot-Lau interferometer into a commercially available X-ray microscope. They demonstrated that this system allows quantitative phase measurements while providing a superior phase sensitivity in comparison with the Zernike phase-contrast mode available on the same instrument [19].

In conventional XTI, the phase shift of the respective sample can be determined by simply integrating the differential phase images [17]. More sophisticated algorithms are needed in order to complete the same task in phase difference microscopy. This is particularly true when the sample is larger than the characteristic length scale on which the interferometer is sensitive to phase differences. In Ref. [12] a discrete shift-and-sum algorithm has been formulated. As an alternative, a deconvolution method based on the Fourier transform was proposed in Ref. [13]. However, these techniques are prone to noise and suffer from deconvolution artefacts [19]. To remedy these issues, Takano et al. [20] developed an iterative deconvolution process. It is based upon repeatedly updating those regions of the spatial frequency domain for which the transfer function of the phase difference process, i.e. the denominator of Eq. (3) in Ref. [20], is zero. The updates are driven by assumptions about the sample, in this case by a support that was generated via filtering and thresholding of an initial phase map obtained through a simple deconvolution.

In this work, we seek to demonstrate an alternative image reconstruction method for X-ray microscopy with single-grating XTI. It is mainly motivated by the observation that the aforementioned phase difference images can only be obtained under a small set of experimental conditions. This results in demands on the interferometer, the sample as well as on the layout and the illumination function of the microscope, hence posing limitations to the applicability of XTI. In the general case, local transverse shifts of the interference fringes – they are usually quantified through the argument of the first order Fourier term of the interference pattern – still depend on phase differences across the wave front, but unlike the special situation of phase difference imaging, this dependency is nonlinear. In order to fully take this into account, the proposed scheme is based on a complete forward model of XTI which is then used to set up a regularized maximum likelihood type cost function. Minimizing the cost function by the application of iterative nonlinear optimization techniques enables the simultaneous reconstruction of both the phase and the transmission image of the sample.

The paper is structured as follows: In Section 2 we review the principles of XTI, with a focus on setups for X-ray microscopy employing a single grating setup. We highlight possible limitations to the applicability of this phase-contrast imaging technique. Section 3 details the proposed maximum likelihood reconstruction. The results are presented in Section 4. Here, the proposed reconstruction method is tested on simulated data. We close the paper with a brief conclusion and outlook in Section 5.

2. Theoretical background

To review the principles of XTI, we consider the X-ray imaging microscope shown in Fig. 1. The focusing device, e.g. a Fresnel zone plate or a compound refractive lens, is illuminated by a quasi-monochromatic and spatially coherent wave. The sample and the X-ray detector are placed in the object and image plane of the lens respectively. Hence, the distance $z_1$ from the sample to the lens and the distance $z_2$ from the lens to the detector satisfy the classical lens law $1/z_1 + 1/z_2 = 1/f$ of geometrical optics [21], with $f$ the focal length of the focusing system. Without the grating, the setup resembles a conventional full-field X-ray microscope. Hence, it would allow the acquisition of magnified absorption images. The magnification $M$ of the microscope is given by the ratio $z_2/z_1$.

Fig. 1. Schematic of an X-ray imaging microscope with XTI.

Download Full Size | PDF

The Talbot effect and, by extension, the image formation process of XTI is governed by Fresnel diffraction [22]. The latter can be regarded as a paraxial approximation of the general Rayleigh-Sommerfeld diffraction formula for scalar waves [21]. Given the values of a complex scalar wave field $E_{\mathrm {in}}(x,y)$ in the plane $z=z_{\mathrm {in}}$, the Fresnel diffracted field $E_{\mathrm {out}}$ in the plane $z=z_{\mathrm {out}}$ is given by

(1)$$E_{\mathrm{out}}(x,y) = \frac{\exp(ik\Delta)}{i\lambda \Delta} \int E_{\mathrm{in}}(x',y') \exp\left[\frac{ik}{2\Delta}\left(\left(x-x'\right)^2 + \left(y-y'\right)^2 \right)\right] \mathrm{d}x' \mathrm{d}y' \;,$$

with $\Delta = z_{\mathrm {out}} - z_{\mathrm {in}}$ the propagation distance, $k=2\pi / \lambda$ the wavenumber, and $\lambda$ the wavelength of the X-rays. Formally, we can express the propagation of a wave field over a distance $\Delta$ in the $z$-direction, as in Eq. (1), through the action of an operator, namely

(2)$$E_{\mathrm{out}} = \mathcal{P}_{\Delta}E_{\mathrm{in}}\;.$$

We now assume that a diffraction grating with a period of length $d$ is located at a distance $R_1$ downstream of the focus point, cf. Figure 1. For simplicity and by virtue of the Fresnel scaling theorem [22], we work out the Fresnel propagation from the grating to the detector plane in an effective geometry, where the illumination by the focus, which is assumed to be equivalent to a point source, is replaced by a plane-wave illumination. This mapping leads to an effective propagation distance $z_{\mathrm {eff}} = z_{\mathrm {gd}} / M_{\mathrm {g}}$ and to demagnified coordinates $x = x_{\mathrm {d}} / M_{\mathrm {g}}$ in the detector plane. Here, the magnification of the grating $M_{\mathrm {g}}$ is given by $R_2/R_1$. The resulting effective wave field $E_{\mathrm {d}}^{\mathrm {eff}}$ in the detector plane can be transformed by means of $E_{\mathrm {d}}^{\mathrm {eff}}(x)/M_{\mathrm {g}} = E_{\mathrm {d}}(M_{\mathrm {g}}x)$ into the actual wave field $E_{\mathrm {d}}$. The former one can be written as

(3)$$E_{\mathrm{d}}^{\mathrm{eff}} = \mathcal{P}_{z_{\mathrm{eff}}}\left[E_{\mathrm{g}} T \right]\;,$$

where $E_{\mathrm {g}}$ is the wave field just in front of the grating and $T$ is the complex transmission function of the grating. Assuming an infinitely large grating that is uniform in the $y$-direction, the complex transmission function can be expressed in form of a Fourier series, i.e.

(4)$$T(x,y) = \sum_n a_n \exp\left(i2\pi\frac{n}{d}x\right) \;.$$

After inserting Eq. (4) into Eq. (3), the resulting diffraction integral can be simplified by first interchanging the order of integration and summation and then completing the square in the exponent of each term of the Fourier series. This yields

(5)$$E_{\mathrm{d}}^{\mathrm{eff}}(x) = \sum_n a_n \exp\left({-}i\pi\frac{n^2\lambda}{d^2} z_{\mathrm{eff}}\right) \exp\left(i2\pi\frac{n}{d}x\right) \mathcal{P}_{z_{\mathrm{eff}}}E_{\mathrm{g}}\left(x - n\frac{\lambda z_{\mathrm{eff}}}{d}\right) \;,$$

where the y-dependence has been omitted for clarity. Also note that $\mathcal {P}_{z_{\mathrm {eff}}}E_{\mathrm {g}}$ is the wave field which would arrive at the detector in the absence of the grating.

For the further analysis of Eq. (5), it is beneficial to write the effective propagation distance as

(6)$$z_{\mathrm{eff}} = m_{\mathrm{t}} \frac{d^2}{\lambda}$$

with $m_{\mathrm {t}}$ the so-called Talbot order. Assuming that $E_{\mathrm {g}}$ is a plane wave, the wave field at a distance $z_{\mathrm {eff}}$ behind the grating is the same as the one directly behind the grating, if $m_{\mathrm {t}}$ is an even integer [22]. In the case that $m_{\mathrm {t}}$ is an odd integer, the field behind the grating is reproduced with a transverse shift of half the period. When using phase gratings, there are no intensity modulations in these distances. However, under certain circumstances, intensity patterns that resemble the shape of the phase grating can be observed in fractional Talbot distances. For a binary $\pi /2$ or a $\pi$-shifting phase grating with a duty-cycle of $0.5$ such self-images manifest if $m_{\mathrm {t}}$ is an odd multiple of $1/2$ or $1/8$ respectively [23].

Moving from the effective plane-wave illumination back to the real situation of spherical-wave illumination, there can be two possible positions $R_1$ for a given grating and a given setup length $R_2$, which will result in the same effective propagation distance $z_{\mathrm {eff}}$. The two positions – they are the solutions to the quadratic equation obtained by combining the conditions $R_2 - R_1 = z_{\mathrm {gd}}$, $z_{\mathrm {eff}} = z_{\mathrm {gd}}/M_{\mathrm {g}}$, and $M_{\mathrm {g}} = R_2/R_1$ – are given by

(7)$$R_1 = \frac{R_2}{2} \left( 1 \pm \sqrt{1 - \frac{4 z_{\mathrm{eff}}}{R_2}} \right)\;.$$

The solution with the minus sign leads to a larger period of the intensity fringe pattern and consequently facilitates its direct resolvability.

The intensity measured at the detector is obtained by taking the absolute square of Eq. (5). This results in

(8)$$\begin{aligned}I(x) = \sum_{n, m} &a_n a_m^{{\ast}} \exp\left({-}i\pi m_{\mathrm{t}} \left(n^2 - m^2\right) \right) \exp\left(i2\pi\frac{n-m}{d}x\right) \\ &\mathcal{P}_{z_{\mathrm{eff}}}E_{\mathrm{g}} \left(x - n \frac{\lambda z_{\mathrm{eff}}}{d}\right) \mathcal{P}_{z_{\mathrm{eff}}}E_{\mathrm{g}} \mkern1mu^{{\ast}}\left(x - m \frac{\lambda z_{\mathrm{eff}}}{d}\right)\;. \end{aligned}$$

We see from Eq. (8) that the single grating Talbot setup acts as a multiple beam interferometer. In other words, the intensity pattern on the detector is formed by pairwise interference of the grating diffraction orders, where the separation distance between neighboring diffraction orders is given by $d_{\mathrm {s}} = \lambda z_{\mathrm {eff}} / d$. This is the fundamental reason why XTI is a phase-sensitive technique. Generally speaking, nearly any deformation to the phase of the wave field will result in alterations of the intensity pattern. A key question for phase-contrast imaging is then how to infer the phase from the intensity measurements. In XTI, one usually tries to extract the first order Fourier term of the intensity fringe pattern and therewith the interference terms between adjacent diffraction orders. This can be done by the fringe scanning method, where a series of images is taken for different transverse positions of the grating [4], or by a suitable single-shot scheme, e.g. the Fourier transform method proposed by Takeda et al. [24]. From Eq. (8), a general expression for the $p$-th order Fourier term of $I(x)$ is obtained by only keeping those terms of the double sum that fulfil $n-m = p$. One gets

(9)$$I^{(p)}(x) = \exp\left(i2\pi \frac{p}{d} x\right) \sum_n C_{n,p,m_{\mathrm{t}}} \mathcal{P}_{z_{\mathrm{eff}}}E_{\mathrm{g}} {\bigg(}x - n d_{\mathrm{s}}{\bigg)} \mathcal{P}_{z_{\mathrm{eff}}}E_{\mathrm{g}} \mkern1mu^{{\ast}}{\bigg(}x - (n-p) d_{\mathrm{s}}{\bigg)}$$

with $C_{n,p,m_{\mathrm {t}}} = a_n a_{n-p}^{\ast } \exp \left (-i\pi m_{\mathrm {t}}p(2n-p)\right )$.

In an imaging microscope, as the one depicted in Fig. 1, assuming an ideal lens, the wave field $\mathcal {P}_{z_{\mathrm {eff}}}E_{\mathrm {g}}$ is a flipped and magnified version of the wave field $E_{\mathrm {s}}$ in the sample plane, i.e.

(10)$$\mathcal{P}_{z_{\mathrm{eff}}}E_{\mathrm{g}}(x) = \frac{M_{\mathrm{g}}}{M} E_{\mathrm{s}}\left(-\frac{M_{\mathrm{g}}}{M}x \right)\;.$$

Here, the inclusion of the factor $M_{\mathrm {g}}$ is a result of the effective geometry. At this point, it should be noted, that the calculations presented thus far and in the following can be generalized to the case of realistic focusing elements. As shown in Ref. [13], the deviations of the focusing element from an ideal lens result in an amplitude spread function for the X-ray microscope, with which the wave field $\mathcal {P}_{z_{\mathrm {eff}}}E_{\mathrm {g}}$ would have to be convolved in Eq. (9). In the projection approximation [22], the modification of the wave field $E_{\mathrm {s}}$ through the insertion of a sample can be obtained by the multiplication with the complex transmission function

(11)$$S(x,y) = \left|S(x,y)\right| \exp\left(i\Phi(x,y)\right)$$

of the sample. The amplitude transmission coefficient $|S|$ and the phase shift $\Phi$ of the sample are determined by the integral of the imaginary part $\beta$ and the decrement $\delta$ of the complex refractive index $n=1-\delta +i\beta$ over the respective path length inside the sample [22], i.e.

(12)$$-\ln{\left|S(x,y)\right|} = k \int_{\mathrm{sample}} \beta(x,y,z) \;\mathrm{d}z$$

and

(13)$$\Phi(x,y) ={-}k \int_{\mathrm{sample}} \delta(x,y,z) \;\mathrm{d}z\;.$$

In total, the wave field $E_{\mathrm {s}}$ directly behind the sample can then be expressed via its amplitude $A_{\mathrm {s}}$ and phase $\varphi _{\mathrm {s}}$ as

(14)$$A_{\mathrm{s}}(x,y) \exp\left(i\varphi_{\mathrm{s}}(x,y)\right) = A_{\mathrm{ref}}(x,y) \left|S(x,y)\right| \exp\left(i\varphi_{\mathrm{ref}}(x,y) + i\Phi(x,y) \right)$$

where $A_{\mathrm {ref}}$ and $\varphi _{\mathrm {ref}}$ denote the amplitude and the phase if no sample is in the beam. Writing the wave field $\mathcal {P}_{z_{\mathrm {eff}}}E_{\mathrm {g}}$ as $A\exp (i\varphi )$ in Eq. (9) – $A$ and $\varphi$ can be linked to $A_{\mathrm {s}}$ and $\varphi _{\mathrm {s}}$ in Eq. (14) by means of Eq. (10) – results in

(15)$$\begin{aligned}I^{(p)}(x) = \exp\left(i2\pi \frac{p}{d} x\right) \sum_n &C_{n,p,m_{\mathrm{t}}} A{\bigg(}x - n d_{\mathrm{s}}{\bigg)} A{\bigg(}x - (n-p) d_{\mathrm{s}}{\bigg)} \\ &\exp\Big( i\varphi{\bigg(}x - n d_{\mathrm{s}}{\bigg)} - i\varphi{\bigg(}x - (n-p) d_{\mathrm{s}}{\bigg)}\Big)\;. \end{aligned}$$

With the terminology introduced above, we can now distinguish two distinct operation regimes of XTI. For the first case, we assume that both the phase $\Phi$ can be approximated by a linear function and that the amplitude transmission coefficient can be assumed to be constant in an interval determined by the separation distance $d_{\mathrm {s}}$ and the highest relevant diffraction order $n_{\mathrm {max}}$, i.e.

(16)$$\begin{aligned} & \Phi(x) - \Phi(x+pd_{\mathrm{s}}) \approx{-}pd_{\mathrm{s}} \frac{\partial \Phi}{\partial x} (x_{\mathrm{c}}) \\ & |S(x)| \approx |S(x_{\mathrm{c}})| \end{aligned} \quad \mathrm{for}\; x \in [x_c - n_{\mathrm{max}}d_{\mathrm{s}}, x_c + n_{\mathrm{max}}d_{\mathrm{s}}]\;.$$

These requirements are usually only fulfilled if $n_{\mathrm {max}}d_{\mathrm {s}}$ is small compared to the spatial resolution. After applying the approximations of Eq. (16) to Eq. (15), a differential phase image can be extracted from the complex arguments of suitable Fourier orders of the intensity measurements $I_{\mathrm {sample}}$ and $I_{\mathrm {ref}}$ with and without the sample in the beam. One finds that

(17)$$\arg\left(I^{(p)}_{\mathrm{sample}}\right)-\arg\left(I^{(p)}_{\mathrm{ref}}\right) \approx{-}pd_{\mathrm{s}} \frac{\partial \Phi}{\partial x}$$

yields an image that is proportional to the component of the gradient of $\Phi$ in the direction perpendicular to the grating bars. In setups for X-ray microscopy, separation distances $d_{\mathrm {s}}$ much larger than the resolution of the detection system can be realized. Consequently, significant changes in the complex transmission function of the sample are to be expected over this length scale and the approximations of Eq. (16) do not hold anymore. In the general case, the complex argument of $I^{(p)}$ is a nonlinear function of phase differences $\varphi (x) - \varphi (x+pd_{\mathrm {s}})$. A particular case arises for the choice of a $\pi /2$-shifting grating with a duty-cycle of $0.5$. Here, all even Fourier coefficients except for $n=0$ vanish, so that only the two terms with $C_{0,1,m_{\mathrm {t}}}$ and $C_{1,1,m_{\mathrm {t}}}$ remain from the series for $I^{(1)}$ in Eq. (15). Additionally, we have $C_{0,1,m_{\mathrm {t}}} = C_{1,1,m_{\mathrm {t}}}$. Thus, sum-to-product identities for the sine and the cosine can be applied to the two remaining terms. Under the condition

(18)$$A(x - d_{\mathrm{s}}) = A(x + d_{\mathrm{s}})\;,$$

which effectively amounts to requiring a constant amplitude $A$, a phase difference image of the sample can then be obtained by combining the complex arguments of sample and reference measurements, analogous to the case in Eq. (17), i.e.

(19)$$\arg\left(I^{(1)}_{\mathrm{sample}}\right)-\arg\left(I^{(1)}_{\mathrm{ref}}\right) = \frac{1}{2} \Big( \Phi(x - d_{\mathrm{s}}) - \Phi(x + d_{\mathrm{s}}) \Big)\;.$$

The reconstruction schemes proposed so far in Refs. [12,13,20] are all based upon Eq. (19) and the fact, that the image obtained thereout linearly depends on the the sample’s phase shift $\Phi$. This allows the direct application of deconvolution techniques. The premise of Eq. (18), which is a necessary condition for the linear dependency to hold, poses limitations to the general applicability of XTI. It does not only imply a negligible absorption due to the sample, i.e. $|S| \approx 1$, but also a uniform amplitude of the X-ray illumination. Wave front measurements at nanofocusing instruments by either ptychography [25] or grating interferometry [26] reveal significant variations in the amplitude of the focused beam. This could lead to erroneous values if the phase shift is retrieved from Eq. (19). Beyond posing limitations to the probing X-ray beam, Eq. (18) hinders the combination of XTI with X-ray projection microscopes, where the sample is placed in the divergent beam downstream of the back focal plane. For such a setup, $\mathcal {P}_{z_{\mathrm {eff}}}E_{\mathrm {g}}$ would not be an image of $E_{\mathrm {s}}$, but the corresponding Fresnel diffraction image. Here, even a purely phase-shifting sample would lead to characteristic interference fringes and hence to variations in $A(x)$. Similarly, unwanted defocusing in an imaging microscope can also result in the nullity of Eq. (19) and hence in image artefacts. Lastly, the usage of other gratings than a $\pi /2$-shifting one may prove beneficial. A $\pi$-shifting grating for example offers a greater flexibility in practice, as its Talbot distances lie closer together. With $I^{(2)}$ being its dominant Fourier order, it also provides a higher frequency Talbot pattern than an interferometer based on a $\pi /2$-shifting grating with a similar separation distance $d_{\mathrm {s}}$. For the single-shot reconstruction method of Ref. [24], this amounts to a larger bandwidth available to the phase and thus to an enhanced spatial resolution. However, an expression similar to Eq. (19) does not exist for $\pi$-shifting gratings in the first place.

3. Reconstruction method

A statistical image reconstruction (SIR) method determines the image that fits best to the input data by minimizing a cost function. Following the work of Fessler [27], SIR methods can be divided into five main components: a discrete representation of the image or the images to be reconstructed, a forward model linking the unknown quantities to the measurements, assumptions about the statistical fluctuations of the measurement around its expectation value, a cost function whose minimization provides an image estimate, and an algorithm to perform the said minimization. In the following, we will outline a SIR method for X-ray microscopy with XTI and shortly address each of these constituents.

For XTI, the quantities to be reconstructed are images of both the amplitude $A$ and the phase $\varphi$ of the wave field $\mathcal {P}_{z_{\mathrm {eff}}}E_{\mathrm {g}}$. A sort of natural representation is obtained by parameterizing both quantities on a finite pixel grid. In this work, the pixel pitch of this grid corresponds to the virtual detector pixel size in the grating plane. Hereinafter, it is denoted by $\Delta x$. Consequently, $A$ and $\varphi$ are represented through the values $A_{i,j}$ and $\varphi _{i,j}$ respectively, where the indices $i$ and $j$ indicate the row and the column of the pixel matrix.

Inserting the discrete representation of $A_{i,j}\exp (i\varphi _{i,j})$ into Eq. (8), we obtain the following forward model for the intensity in the detector plane:

(20)$$\begin{aligned}I_{i,j}(A,\varphi) = \sum_{n,m} &a_n a_m^{{\ast}} \exp\left({-}i\pi m_{\mathrm{t}} \left(n^2 - m^2\right) \right) \exp\left(i2\pi\frac{n-m}{d}x_j\right) \\ &\mathrm{sinc}\left(\frac{n-m}{d}\Delta x\right) A_{i,j_n}\exp\left(i\varphi_{i,j_n}\right) A_{i,j_m}\exp\left({-}i\varphi_{i,j_m}\right)\;. \end{aligned}$$

Here $x_j$ denotes the x-coordinate of the pixel centers in the $j$-th column. In order to model the sampling of the continuous intensity $I(x)$ by the detector pixels, the term with the sinc function had to be added compared to Eq. (8) [28]. The indices $j_n$ are given by $j - \left \lfloor n d_{\mathrm {s}} / \Delta x + 0.5 \right \rfloor$, where $\lfloor . \rfloor$ denotes the floor function. This rounding method effectively corresponds to a nearest-neighbor interpolation of $A(x-nd_{\mathrm {s}})$ and $\varphi (x-nd_{\mathrm {s}})$ on the pixel grid. It was chosen for the first proof-of-principle studies in order to keep the forward model and thus the cost function and above all its minimization as simple as possible.

The central component of the cost function is a data term. Its purpose is to evaluate the likeliness that a pair of $A$ and $\varphi$ belongs to the measured intensity $\tilde {I}$. Under the assumption that the statistical fluctuations in the measured data are dominated by the Poissonian noise of the detected number of photons, a suitable data term for the cost function is given by the negative log-likelihood [29]

(21)$$L\left(A,\varphi | \tilde{I}\right) = \sum_{i,j} -\tilde{I}_{i,j} \log{\bigg(}I_{i,j}(A,\varphi){\bigg)} + I_{i,j}(A,\varphi)\;.$$

It is to be noted, that all terms which do not depend on the intensities $I_{i,j}$ estimated by Eq. (20) and hence do not influence the later minimization are omitted in Eq. (21). We would like to point out that a weighted least squares type cost function might be more appropriate as soon as noise sources other than photon shot noise become relevant.

One downside of the proposed forward model is that the cost function is not convex. Rather, it does not possess a strict global minimum. One cause for the last point is the possible occurrence of phase wrapping. Furthermore, a closer inspection of the continuous model in Eq. (8) reveals that spatial frequencies which are an integer multiple of $1/d_{\mathrm {s}}$ can be added to the phase $\varphi (x)$ without affecting the expected intensity. Consequently, the data term requires regularization. This amounts to imposing a certain degree of a priori information onto the parameters of the forward model. We decided to include regularization terms in form of an $\epsilon$-smoothed total variation (TV) norm [30], i.e.

(22)$$R_{\mathrm{TV}}(f,\beta_{\mathrm{x}},\beta_{\mathrm{y}},\epsilon) = \sum_{i,j} \beta_{\mathrm{y}} \sqrt{|f_{i+1,j}-f_{i,j}|^2+\epsilon} + \beta_{\mathrm{x}} \sqrt{|f_{i,j+1} - f_{i,j}|^2+\epsilon}\;.$$

The purpose of such a regularization term is to penalize differences between neighboring pixel values of an image $f$. This will mainly be in conflict with the data term and the parameters $\beta$ can be used to adjust the trade-off between both terms. Since the directionality of the grating introduces an asymmetry between the x- and y-direction, we opted for an anisotropic version of the TV norm. The complete cost function of our SIR method reads

(23)$$C = L\Big(A,\varphi | \tilde{I}\Big) + R_{\mathrm{TV}}\Big(A,\beta_{\mathrm{x}}^{\mathrm{A}},\beta_{\mathrm{y}}^{\mathrm{A}},\epsilon^{\mathrm{A}}\Big) + R_{\mathrm{TV}}\Big(\varphi,\beta_{\mathrm{x}}^{\mathrm{\varphi}},\beta_{\mathrm{y}}^{\mathrm{\varphi}},\epsilon^{\mathrm{\varphi}}\Big)\;.$$

An estimate for the amplitude $A$ and the phase $\varphi$ of the wave front is then obtained by minimizing the cost function of Eq. (23) with respect to the entries $A_{i,j}$ and $\varphi _{i,j}$ of the respective pixel matrix. The minimization is carried out with an implementation of the L-BFGS-B algorithm [31]. L-BFGS-B is a limited memory version of the Broyden-Fletcher-Goldfarb-Shanno algorithm (BFGS), that can also handle bound constraints on its variables. As a quasi-Newton method, L-BFGS-B uses an approximation of the inverse Hessian matrix in order to determine a descent direction in every iteration. To limit the computer memory required, the approximation of the inverse Hessian is stored implicitly by a certain number of vectors [32]. This makes L-BFGS-B applicable to large-scale optimization problems. At this point, it should be emphasized, that the main focus of this work is the derivation of a forward model based on physically modelling the imaging system and to illustrate its suitability for a cost function which takes the nonlinear relationship between the interference fringe shift and the sample phase shift into account. An investigation of the performance of the applied optimization algorithm is beyond the scope of this paper. We are aware that more suitable optimization techniques might exist. Yet the L-BFGS-B algorithm minimized Eq. (23) reliably enough for the presented investigations.

In order to isolate the phase shift $\Phi$ and the amplitude transmission coefficient $|S|$ of a sample, the recovered phases $\varphi _{\mathrm {sample}}$ and $\varphi _{\mathrm {ref}}$ and amplitudes $A_{\mathrm {sample}}$ and $A_{\mathrm {ref}}$ of measurements with and without the sample in the beam have to be combined. It can be inferred from Eq. (14) that $\Phi = \varphi _{\mathrm {sample}} - \varphi _{\mathrm {ref}}$ and $|S| = A_{\mathrm {sample}}/A_{\mathrm {ref}}$. The forward model, the cost function, and its gradient were implemented using the Python package JAX [33].

4. Results

To validate the proposed SIR method, numerical simulations of XTI in an X-ray microscope were conducted. The simulations are based on a numerical wave field propagation based on scalar diffraction theory [34]. The transmission of the wave field through the sample and the optical elements is treated in the projection approximation. We assume a perfectly coherent wave field and an ideal detector without any crosstalk. We only consider the Poissonian photon shot noise. The X-ray energy was fixed at $11\,\mathrm {keV}$ in all simulations. Two different phantoms were used to test the algorithm. The first phantom consists of polystyrene spheres with diameters of $8\,\mathrm {\mu m}$, $12\,\mathrm {\mu m}$, and $18\,\mathrm {\mu m}$. With a minimal amplitude transmission coefficient of $|S|\approx 0.998$, they can be treated in good approximation as pure phase objects with phase values $\varphi _{i,j} \in [-1.86,0]\,\mathrm {rad}$. The second phantom is a $1\,\mathrm {\mu m}$ thick Siemens star made out of tantalum. It has a diameter of $100\,\mathrm {\mu m}$ and $20$ spokes. With $|S|\approx 0.86$, the Siemens star cannot be considered a pure phase object. In the simulations, an X-ray detector with $2048 \times 2048\,\mathrm {pixel}^2$ with a pixel size of $6.5\,\mathrm {\mu m}$ was used. The magnification of the X-ray microscope was set to $M=88.75$ which results in an effective pixel size of roughly $73\,\mathrm {nm}$ in the sample plane. There, the phantoms were generated with an initial size of $2048 \times 2048\,\mathrm {pixel}^2$. To avoid artefacts in the numerical wave propagation, both phantoms were first padded to a size of $4096 \times 4096\,\mathrm {pixel}^2$. Afterwards, they were upsampled by a factor of six to guarantee a proper sampling of the grating transmission function. The resulting intensity values in the detector plane were then rebinned to the pixel size before Poissonian noise was added for different mean fluences.

Figs. 2(a) and 2(d) shows the phase shift of the two samples for a central region of $1748 \times 1748\,\mathrm {pixel}^2$ and $1448 \times 1448\,\mathrm {pixel}^2$ respectively. The reconstructions were later carried out for these central regions. Two exemplary realizations of simulated intensity measurements for a Talbot interferometer based on a $\pi /2$-shifting grating are depicted in Figs. 2(b) and (e). Therein, the mean fluence was $N\approx 110$ photons per pixel. For a better visualization, only a small detector region of $200\times 200\,\mathrm {pixel}^2$ is displayed. The corresponding sections of the phantoms are marked by the red squares in (a,d). For these regions, phase difference images are shown in Figs. 2(c) and (f). They were retrieved by means of Eq. (19), where the first order Fourier terms have been extracted via the fringe scanning method. For that, $42$ single phase steps have been simulated. As expected from the right hand side of Eq. (19), the phase difference image of the polystyrene sphere in Fig. 2 (f) shows two phase maps of the sphere which are separated by a distance of $2d_{\mathrm {s}}$ and of different sign. The right phase map is due to the interference of the positive first and the zeroth diffraction order, the left one due to the zeroth and the negative first order. For the Siemens star in Fig. 2 (c), the right hand side of Eq. (19) does no longer hold. Here, the non-negligible absorption leads to visible deviations from the actual phase difference.

Fig. 2. Phantoms used in the numerical simulations. The respective phase shift of the phantoms is depicted in (a) and (d). For the areas marked by the red squares in (a) and (d), exemplary simulations of the measured intensities with a mean fluence of $N \approx 110$ photons per pixel are shown in (b) and (e). Phase difference images of these regions are displayed in (c) and (f). In (c) the phase difference image is flawed by the non-negligible absorption of the Siemens star. The scale bar indicates $15\,\mathrm {\mu m}$ length in the sample plane (a,d) and $200\,\mathrm {\mu m}$ in the detector plane (b,c,e,f) respectively.

Download Full Size | PDF

The simulations were carried out for two different XTI setups. In Subsection 4.1, we first consider the case of a $\pi /2$-shifting grating. This allows a direct comparison of the SIR method with the already existing deconvolution-based techniques. In Subsection 4.2, we then investigate, by means of a $\pi$-shifting grating, a case for which the possibility of reconstruction via deconvolution does not exist. The parameters of both Talbot interferometers and the respective X-ray microscope, cf. Figure 1, are listed in Table 1.

Table 1. Overview of the parameters of both simulated XTI setups.

View Table | View all tables in this article

4.1 XTI with a $\pi /2$-shifting grating

The simulated intensity images of both phantoms, as exemplarily shown in Figs. 2(b) and (e), have been used in three different reconstruction schemes for the $\pi /2$-shifting setup. The iterative deconvolution method recently proposed by Takano et al. [20] served as a comparison to our SIR technique. The deconvolution approach by Takano et al. takes phase difference images, as the ones shown in Figs. 2(c) and (f) as input with the aim of separating the multiple visible and separated phase maps of the sample. Our proposed iterative scheme has in principle the same goal, but operates directly on the intensity fringe patterns. To apply the deconvolution method, the first order Fourier terms $I^{(1)}$ of the sample and reference measurement have to be extracted from the simulated intensity measurements first, cf. Eq. (19). This was realized once with the fringe scanning approach [4] on the full data set of $42$ phase steps and through the application of the Fourier transform method [24] on a single exposure. While an extension of our SIR method to fringe scanning data sets is conceptually possible, it was only applied to single exposures in this work. The minimization of the cost function was terminated as soon as the norm of the gradient fell below a threshold of $10^{-3}$ or after a maximum of $1000$ iterations. At this point, no more significant reductions of the cost function could be achieved. Two uniformly initialized arrays are used as the starting values for the amplitude $A_{i,j}$ and the phase $\varphi _{i,j}$. The amplitude was set to the square root of the mean measured intensity while the phase was set to zero. For the Siemens star phantom, we initially performed $500$ iterations, where a non-weighted least squares function, i.e. $\sum _{i,j} (I_{i,j} - \tilde {I}_{i,j})^2$, has been employed as the data term of our cost function instead of the negative log-likelihood. This improved the convergence of the method and prevented stagnation of the minimization in local minima or saddle points, that still significantly differed from the ground truth. The parameters of the TV regularization were set to $\beta _{\mathrm {x}} = 0.5$, $\beta _{\mathrm {y}} = 0.05$, and $\epsilon = 10^{-4}$ for both $A$ and $\varphi$. The iterative reconstructions were each run on a single socket compute node equipped with an Intel Xeon E3-1240 v5 processor. In Table 2, we list the mean run time $t_{\mathrm {mean}}$ for one iteration of the L-BFGS-B algorithm in dependence of the highest diffraction order $n_{\mathrm {max}}$ considered in the forward model of Eq. (20). In each case, the mean run time was determined by running $100$ iterations of the minimization process for the spheres phantom with its reconstruction area of $1448 \times 1448\,\mathrm {pixel}^2$. The values agree with the expected quadratic growth of $t_{\mathrm {mean}}$ with $n_{\mathrm {max}}$. To keep the computational costs acceptable, only diffraction orders up to $n_{\mathrm {max}} = 9$ were considered in the full reconstructions.

Table 2. Mean run time $t_{\mathrm {mean}}$ per iteration in dependence of the highest diffraction order $n_{\mathrm {max}}$ used in the forward model. The values were determined for a reconstruction area of $1448 \times 1448\,\mathrm {pixel}^2$.

View Table | View all tables in this article

The performance of the three reconstruction approaches is surveyed in Fig. 3 using the simulated measurements with a mean fluence of $N \approx 110$ photons per pixel for both phantoms. The reconstructed phase shifts of the polystyrene spheres are depicted in (a,b,c). Due to the $42$ times higher photon statistics, the result of the fringe scanning deconvolution (a) shows a smoother background and less distortions of the retrieved phase shift $\Phi$ in comparison with the single-shot deconvolution method (b). In (b) single statistical fluctuations of $I^{(1)}_{\mathrm {sample}}$ and $I^{(1)}_{\mathrm {ref}}$ are horizontally distributed over the width of the phase image through the deconvolution which leads to periodic artefacts. The reconstruction by the SIR method (c) is visually comparable to (a), albeit being retrieved from the same single exposure as (b). To a large part, this can be attributed to the employed regularization favoring smooth results. Further differences between the three reconstructions are revealed by inspecting the respective reconstruction error $\Delta \Phi = \Phi (\mathrm {reconstruction})-\Phi (\mathrm {phantom})$ of the phase shift shown in (d,e,f). Compared to the two deconvolution-based reconstructions (d,e), the errors for the SIR method in (f) are approximately constant over larger areas of the phase image. For instance, the reconstructed phase shift of all the spheres is slightly too large in a region around the respective center. The three reconstructions of the Siemens star are depicted in (g,h,i). Unlike the reconstruction of the spheres, the fringe scanning deconvolution method now also exhibits visible stripe artefacts. In contrast to the noise-induced stripes of (b) and (e), their cause lies in the non-negligible absorption of the Siemens star. Consequently, their occurrence is limited to the vertical extent of the Siemens star. Here, absorption by the sample leads to inconsistencies with the right-hand side of Eq. (19) and thus with the foundation of the deconvolution approach. The SIR method (i) does not suffer from effects of absorption, as they are consistently incorporated in the forward model. The reconstruction errors of the SIR method (l) are again spread out by the regularization. Due to the directional dependence of the phase sensitivity in XTI, the errors are the smaller the more the spokes are aligned along the grating bars. Larger isolated errors can be observed at the edges of the Siemens star. They are conjectured to be caused by the inexact representation of the shifted copies $A(x-nd_{\mathrm {s}})$ and $\varphi (x-nd_{\mathrm {s}})$ in the discrete parametrization of the forward model, cf. Eq. (20). Compared to the errors of the fringe scanning deconvolution (j), larger errors manifest in the center of the Siemens star for both single-shot reconstructions (k,l). This suggests a worse spatial resolution for these schemes.

Fig. 3. Reconstructions for the $\pi /2$-shifting setup with $N \approx 110$ photons per pixel. The retrieved phase shift for different methods is shown in (a,b,c) for the spheres and in (g,h,i) for the Siemens star. The reconstruction errors with respect to the phantoms are depicted in (d,e,f) and (j,k,l) respectively. The first (a,d,g,j) and second column (b,e,h,g) represent the results of the deconvolution method on the full fringe scanning data set and in conjunction with a single-shot approach respectively. The third column (c,f,i,l) displays the results of our SIR method. The scale bar indicates $15\,\mathrm {\mu m}$ length in all panels.

Download Full Size | PDF

In order to test the behavior of the three schemes at different noise levels, Poissonian noise for mean fluences ranging from $N = 10$ to $N = 10^6$ photons per pixel has been added to the simulated intensity measurements. To assess the quality of the subsequent reconstructions, we employ the root-mean-square deviation

(24)$$\mathrm{RMSD}(\Phi) = \left( \frac{1}{N_{\mathrm{px}}} \sum_{\mathrm{pixel}} \Delta\Phi^2 \right)^{-\frac{1}{2}}$$

of the retrieved phase image as error metric. Here, ten realizations have been generated and reconstructed for each mean fluence in order to estimate statistical fluctuations. Figs. 4(a) and 4(b) show the results. For the spheres phantom (a), the overall reconstruction error of the SIR method is smaller at low fluences than those of the deconvolution based approaches. The fringe-scanning deconvolution becomes superior in terms of the RMSD at a mean fluence of $N \approx 800$ photons per pixel. For the single-shot version this crossover takes place at $N \approx 2\cdot 10^4$ photons per pixel. While the RMSD is decreasing over the entire simulated fluence range for the fringe scanning deconvolution, both the single-shot deconvolution and the SIR method show a clear end of the respective decline. In case of the single-shot deconvolution, a possible explanation for this stagnation is the limited spatial resolution inherent to the Fourier transform method [24] for the extraction of the first order Fourier term. Thereby, it is assumed that the spatial frequency bandwidth of the sample does not exceed the demagnified fringe period in the sample plane. Spatial frequencies outside this range cannot be properly recovered, but rather result in erroneous values of $I^{(1)}$. Regarding the SIR method, the lack of improvement for higher fluences can be attributed to the regularization and to remaining inaccuracies of the forward model, i.e. the already mentioned interpolation of the shifted wave fronts and the omission of diffraction orders beyond $n_{\mathrm {max}}$. For the Siemens star (b), the RMSD curve of the fringe scanning deconvolution now also shows a stagnation. It is caused by the absorption of the sample, which leads to inconsistencies with the deconvolution model. Consequently, at some point, the reconstruction quality does not improve anymore despite an increasing signal-to-noise ratio in the original data. Here, the SIR method features the smallest RMSD over the considered fluence range.

Fig. 4. Root-mean-square deviation as a function of the mean fluence. (a) Results for the spheres phantom. (b) Results for the Siemens star. (c) Achieved resolution as quantified by the FRC criterion for the Siemens star. Ten realizations were generated for each photon number. The error bars indicate the respective mean value plus or minus one standard deviation.

Download Full Size | PDF

Having already mentioned the limited spatial resolution of the single-shot deconvolution method, we now come to a quantitative analysis. To this end, we employ the Fourier ring correlation (FRC), where the reconstructed phase shift is correlated with the ground truth of the respective phantom [35]. A resolution criterion of the FRC is the point at which the correlation coefficient, seen as a function of spatial frequency, falls below a predefined $1$-bit threshold. For the Siemens star phantom, this point, the so-called FRC crossover, has been computed for all fluences and all noise realizations. The results are displayed in Fig. 4 (c). Here, the spatial frequencies are given in periods per pixel, which we abbreviate with $\mathrm {px}^{-1}$. The single-shot deconvolution exhibits a stagnation of the resolution with a maximum value of $0.11\,\mathrm {px}^{-1}$. This is in accordance with the available bandwidth in the Fourier transform method [24]. For the setup’s fringe period of length $d_{\mathrm {f}} = 7\,\mathrm {px}$, it amounts to a maximal resolution of $1/(2d_{\mathrm {f}}) \approx 0.071\,\mathrm {px}^{-1}$ perpendicular to the grating bars. The higher value of the FRC criterion can be explained by a better resolution in the direction parallel to the grating bars. The fringe-scanning deconvolution and the SIR method both reach the full resolution of $0.5\,\mathrm {px}^{-1}$. While the superior behavior of the SIR method might be a result of the fact, that the Siemens star is a beneficial sample in conjunction with the edge-preserving properties of the TV regularization [36], our results suggest that the SIR method enables a significant improvement of the spatial resolution for single-shot XTI. Consequently, the reconstruction errors observed in Fig. 3 (l) for the center of the Siemens star might be avoided with different regularization strategies.

4.2 XTI with a $\pi$-shifting grating

We now turn to the setup based on a $\pi$-shifting grating. In this case, it is not formally possible to establish a deconvolution ansatz similar to Eq. (19) as for the situation of a $\pi /2$-shifting grating. Our SIR method on the other hand is not affected by such a priori restrictions. Its adaptation to the new conditions is rather a mere update of the setup parameters listed in Table 1. The optimization of the cost function is again carried out as described in Section 4.1. The retrieved images of the amplitude transmission coefficient $|S|$ and the phase shift $\Phi$ of the Siemens star are shown in Fig. 5. The employed simulated raw data again had a mean fluence of $N \approx 110$ photons per pixel. Similar to the reconstructions for the $\pi /2$-shifting setup shown in Fig. 3, slight reconstruction errors can be observed around the center of the Siemens star. More strikingly, effects from the shifted wave fields due to the grating diffraction orders remain in the image of $|S|$. We suspect, that they are a result of diffraction orders with $|n|>9$, which are not taken into account in the forward model. Overall, the SIR method still provides a quantitative reconstruction reasonably close to the ground truth with RMSD values of $0.094$ for $\Phi$ and $0.020$ for $|S|$.

Fig. 5. Results obtained for the $\pi$-shifting setup and for a mean fluence of $N \approx 110$ photons per pixel. (a) The reconstructed amplitude transmission coefficient. (b) The reconstructed phase shift. The scale bar indicates $15\,\mathrm {\mu m}$ length respectively.

Download Full Size | PDF

5. Conclusion and outlook

In this work, we reviewed the image formation for X-ray Talbot interferometry in an X-ray microscope by means of scalar diffraction theory. In this regime, the separation distance of adjacent grating diffraction orders on the detector is typically larger than the spatial resolution of the imaging system. Due to that, instead of a differential phase image, a phase difference image can be obtained, from which the phase shift of the sample can be retrieved by a deconvolution process. However, the occurrence of an accurate phase difference image is bound to certain experimental conditions. The Talbot interferometer should be based on a $\pi /2$-shifting grating and the amplitude of the wave field should be constant. Away from these conditions, the available phase contrast information depends nonlinearly on phase differences across the wave field, and deconvolution-based reconstruction techniques yield erroneous results.

Based on our theoretical treatment of the image formation, we developed a statistical image reconstruction method which follows a maximum likelihood approach. The method allows the simultaneous retrieval of both the amplitude and the phase of a wave front from a single measured intensity fringe pattern. Consequently, the amplitude transmission coefficient and the phase shift of an imaged sample can be obtained by combining the respective results of a sample and a reference measurement. We validated our method on simulated data including Poissonian noise over a wide range of mean photon numbers. Even though the employed cost function is not convex and no guarantees for proper convergence can be given, our method yielded results that were reasonably close to the respective ground truth for all simulated samples and setups.

In particular, we investigated two cases which do not fulfill the premises of a deconvolution-based reconstruction. In the first case, non-negligible absorption of the sample led to variations in the amplitude of the wave field and thereby to errors of the extracted phase difference image, cf. Figure 2 (c). Here, the statistical image reconstruction method displayed a superior reconstruction quality over a wide range of photon noise. In the second case, a Talbot interferometer based on a $\pi$-shifting grating was employed. While our method can be easily applied to this setup, the corresponding derivation of a deconvolution model is not possible from the outset.

Additionally, our results indicate further advantages of the statistical image reconstruction method for single-shot imaging with Talbot interferometry. Since it works directly on the measured intensity values, no extraction of particular Fourier orders is required. Contrary to previous single-shot schemes, this means that no window functions have to be applied, which suppress or cut off the higher spatial frequencies and hence limit the achievable spatial resolution. With the advent of hard X-ray free electron lasers, this is especially relevant for future single pulse imaging experiments where the dynamics of matter under extreme conditions can be probed with a sub-picosecond temporal resolution [37].

In conclusion, our findings point out a broader and more flexible applicability of X-ray Talbot interferometry, both for imaging with X-ray microscopes, but also in the field of wave front sensing. There it could allow the retrieval of an X-ray beam’s amplitude without having to limit the sensitivity to wave front aberrations [26]. One conceivable use case in imaging, is the application of our proposed reconstruction method in conjunction with an X-ray projection microscope, where it could enable the retrieval of amplitude and phase of Fresnel diffraction images.

Funding

Erlangen Centre for Astroparticle Physics.

Acknowledgments

We thank Matthew James Johnson for the kind replies and the quick bug fixing regarding some minor problems with JAX.

Disclosures

The authors declare no conflicts of interest.

References

1. H. F. Talbot, “Lxxvi. facts relating to optical science. no. iv,” The London, Edinburgh, Dublin Philos. Mag. J. Sci. 9(56), 401–407 (1836). [CrossRef]

2. K. Patorski, “I the self-imaging phenomenon and its applications,” in Progress in Optics, vol. 27, (Elsevier, 1989), pp. 1–108.

3. A. Momose, S. Kawamoto, I. Koyama, Y. Hamaishi, K. Takai, and Y. Suzuki, “Demonstration of x-ray talbot interferometry,” Jpn. J. Appl. Phys. 42(Part 2, No. 7B), L866–L868 (2003). [CrossRef]

4. T. Weitkamp, A. Diaz, C. David, F. Pfeiffer, M. Stampanoni, P. Cloetens, and E. Ziegler, “X-ray phase imaging with a grating interferometer,” Opt. Express 13(16), 6296–6304 (2005). [CrossRef]

5. B. J. McMorran and A. D. Cronin, “An electron talbot interferometer,” New J. Phys. 11(3), 033021 (2009). [CrossRef]

6. F. S. Yasin, T. R. Harvey, J. J. Chess, J. S. Pierce, and B. J. McMorran, “Path-separated electron interferometry in a scanning transmission electron microscope,” J. Phys. D: Appl. Phys. 51(20), 205104 (2018). [CrossRef]

7. M. Strobl, C. Grünzweig, A. Hilger, I. Manke, N. Kardjilov, C. David, and F. Pfeiffer, “Neutron dark-field tomography,” Phys. Rev. Lett. 101(12), 123902 (2008). [CrossRef]

8. I. Manke, N. Kardjilov, R. Schäfer, A. Hilger, M. Strobl, M. Dawson, C. Grünzweig, G. Behr, M. Hentschel, C. David, A. Kupsch, A. Lange, and J. Banhart, “Three-dimensional imaging of magnetic domains,” Nat. Commun. 1(1), 125 (2010). [CrossRef]

9. B. Brezger, L. Hackermüller, S. Uttenthaler, J. Petschinka, M. Arndt, and A. Zeilinger, “Matter-wave interferometer for large molecules,” Phys. Rev. Lett. 88(10), 100404 (2002). [CrossRef]

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

11. Y. Takeda, W. Yashiro, T. Hattori, A. Takeuchi, Y. Suzuki, and A. Momose, “Differential phase x-ray imaging microscopy with x-ray talbot interferometer,” Appl. Phys. Express 1, 117002 (2008). [CrossRef]

12. W. Yashiro, Y. Takeda, A. Takeuchi, Y. Suzuki, and A. Momose, “Hard-x-ray phase-difference microscopy using a fresnel zone plate and a transmission grating,” Phys. Rev. Lett. 103(18), 180801 (2009). [CrossRef]

13. W. Yashiro, S. Harasse, A. Takeuchi, Y. Suzuki, and A. Momose, “Hard-x-ray phase-imaging microscopy using the self-imaging phenomenon of a transmission grating,” Phys. Rev. A 82(4), 043822 (2010). [CrossRef]

14. H. Kuwabara, W. Yashiro, S. Harasse, H. Mizutani, and A. Momose, “Hard-x-ray phase-difference microscopy with a low-brilliance laboratory x-ray source,” Appl. Phys. Express 4(6), 062502 (2011). [CrossRef]

15. S. Berujon, H. Wang, I. Pape, K. Sawhney, S. Rutishauser, and C. David, “X-ray submicrometer phase contrast imaging with a fresnel zone plate and a two dimensional grating interferometer,” Opt. Lett. 37(10), 1622–1624 (2012). [CrossRef]

16. M. Engelhardt, J. Baumann, M. Schuster, C. Kottler, F. Pfeiffer, O. Bunk, and C. David, “High-resolution differential phase contrast imaging using a magnifying projection geometry with a microfocus x-ray source,” Appl. Phys. Lett. 90(22), 224101 (2007). [CrossRef]

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

18. H. Takano, Y. Wu, and A. Momose, “Development of full-field x-ray phase-tomographic microscope based on laboratory x-ray source,” Proc. SPIE 10391, 1039110 (2017). [CrossRef]

19. H. Takano, Y. Wu, J. Irwin, S. Maderych, M. Leibowitz, A. Tkachuk, A. Kumar, B. Hornberger, and A. Momose, “Comparison of image properties in full-field phase x-ray microscopes based on grating interferometry and zernike’s phase contrast optics,” Appl. Phys. Lett. 113(6), 063105 (2018). [CrossRef]

20. H. Takano, K. Hashimoto, Y. Nagatani, J. Irwin, L. Omlor, A. Kumar, A. Tkachuk, Y. Wu, and A. Momose, “Improvement in quantitative phase mapping by a hard x-ray microscope equipped with a lau interferometer,” Optica 6(8), 1012–1015 (2019). [CrossRef]

21. J. W. Goodman, Introduction to Fourier Optics (Roberts and Company Publishers, 2005).

22. D. Paganin, Coherent X-ray Optics (Oxford University, 2006).

23. T. J. Suleski, “Generation of lohmann images from binary-phase talbot array illuminators,” Appl. Opt. 36(20), 4686–4691 (1997). [CrossRef]

24. M. Takeda, H. Ina, and S. Kobayashi, “Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry,” J. Opt. Soc. Am. 72(1), 156–160 (1982). [CrossRef]

25. A. Schropp, R. Hoppe, V. Meier, J. Patommel, F. Seiboth, H. J. Lee, B. Nagler, E. C. Galtier, B. Arnold, U. Zastrau, J. B. Hastings, D. Nilsson, F. Uhlén, U. Vogt, H. M. Hertz, and C. G. Schroer, “Full spatial characterization of a nanofocused x-ray free-electron laser beam by ptychographic imaging,” Sci. Rep. 3(1), 1633 (2013). [CrossRef]

26. M. Seaberg, R. Cojocaru, S. Berujon, E. Ziegler, A. Jaggi, J. Krempasky, F. Seiboth, A. Aquila, Y. Liu, A. Sakdinawat, H. J. Lee, U. Flechsig, L. Patthey, F. Koch, G. Seniutinas, C. David, D. Zhu, L. Mikeš, M. Makita, T. Koyama, A. P. Mancuso, H. N. Chapman, and P. Vagovič, “Wavefront sensing at x-ray free-electron lasers,” J. Synchrotron Radiat. 26(4), 1115–1126 (2019). [CrossRef]

27. J. A. Fessler, “Penalized weighted least-squares image reconstruction for positron emission tomography,” IEEE Trans. Med. Imaging 13(2), 290–300 (1994). [CrossRef]

28. S. K. Park, R. Schowengerdt, and M.-A. Kaczynski, “Modulation-transfer-function analysis for sampled image systems,” Appl. Opt. 23(15), 2572–2582 (1984). [CrossRef]

29. J. M. Bardsley and J. Goldes, “Regularization parameter selection methods for ill-posed poisson maximum likelihood estimation,” Inverse Probl. 25(9), 095005 (2009). [CrossRef]

30. R. Acar and C. R. Vogel, “Analysis of bounded variation penalty methods for ill-posed problems,” Inverse Probl. 10(6), 1217–1229 (1994). [CrossRef]

31. C. Zhu, R. H. Byrd, P. Lu, and J. Nocedal, “Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization,” ACM Trans. Math. Softw. 23(4), 550–560 (1997). [CrossRef]

32. J. Nocedal and S. Wright, Numerical Optimization (Springer Science & Business Media, 2006).

33. J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, and S. Wanderman-Milne, “JAX: composable transformations of Python+NumPy programs,” http://github.com/google/jax (2018). Version 0.1.40.

34. A. Ritter, P. Bartl, F. Bayer, K. C. Gödel, W. Haas, T. Michel, G. Pelzer, J. Rieger, T. Weber, A. Zang, and G. Anton, “Simulation framework for coherent and incoherent x-ray imaging and its application in talbot-lau dark-field imaging,” Opt. Express 22(19), 23276–23289 (2014). [CrossRef]

35. M. Van Heel and M. Schatz, “Fourier shell correlation threshold criteria,” J. Struct. Biol. 151(3), 250–262 (2005). [CrossRef]

36. D. Strong and T. Chan, “Edge-preserving and scale-dependent properties of total variation regularization,” Inverse Probl. 19(6), S165–S187 (2003). [CrossRef]

37. A. Schropp, R. Hoppe, V. Meier, J. Patommel, F. Seiboth, Y. Ping, D. G. Hicks, M. A. Beckwith, G. W. Collins, A. Higginbotham, J. S. Wark, H. J. Lee, B. Nagler, E. C. Galtier, B. Arnold, U. Zastrau, J. B. Hastings, and C. G. Schroer, “Imaging shock waves in diamond with both high temporal and spatial resolution at an xfel,” Sci. Rep. 5(1), 11089 (2015). [CrossRef]

	Talbot interferometer			X-ray microscope
Setup	$d$	$m_{t}$	$R_{1}$	$f$	$z_{1}$	$z_{2}$
$π / 2$ -shifting	$10 μ m$	$7 / 2$	$3.98 m$	$204.1 m m$	$206.4 m m$	$18.32 m$
$π$ -shifting	$10 μ m$	$11 / 8$	$1.37 m$	$140.6 m m$	$142.2 m m$	$12.62 m$

	Talbot interferometer			X-ray microscope
Setup	$d$	$m_{t}$	$R_{1}$	$f$	$z_{1}$	$z_{2}$
$π / 2$ -shifting	$10 μ m$	$7 / 2$	$3.98 m$	$204.1 m m$	$206.4 m m$	$18.32 m$
$π$ -shifting	$10 μ m$	$11 / 8$	$1.37 m$	$140.6 m m$	$142.2 m m$	$12.62 m$

Maximum likelihood reconstruction for grating-based X-ray microscopy

Abstract

1. Introduction

2. Theoretical background

3. Reconstruction method

4. Results

4.1 XTI with a $\pi /2$-shifting grating

4.2 XTI with a $\pi$-shifting grating

5. Conclusion and outlook

Funding

Acknowledgments

Disclosures

References

Cited By

Figures (5)

Tables (2)

Equations (24)

Optics Express