1. Introduction

Among the techniques for X-ray phase imaging, grating interferometers have become popular in the recent years, offering advantages like the possibility to function with a compact laboratory source and broad energy bandwidth [1–7]. Based on the diffraction of X-rays by a periodic grating, the deflection of the rays by refraction at the imaged object is measured. In this technique, an image $\Phi _u$ of the derivative of the object phase shift along one direction $\vec {u}$ is obtained.

This study deals with the estimation of the phase shift image $\Phi$ from measurements of its derivative. Obtaining a clear phase image is an important addition to X-ray grating interferometry, especially when comparisons with the absorption or visibility contrasts are needed. The common approach of directly integrating $\Phi _u$ along $\vec {u}$ by a cumulative sum or by linear filtering in frequency space gives a poor estimation of the phase shift, with typical line artifacts in the $\vec {u}$ direction. This is explained by the unknown integration constants that are missing from the differential phase image and also by the noisy measurements. Clearly, the problem of estimating $\Phi$ from $\Phi _u$ is ill-posed. Some solutions to this problem have been proposed in the literature [8–12], often involving a modification of the experimental setup to measure the second coordinate of the differential phase, thus allowing an easier estimation. The consecutive acquisition of differential phase in both directions, in order to obtain a complete information about the gradient, has been performed [9]. An experimental setup based on two-dimensional gratings has then been proposed [10], allowing the simultaneous measurement of the complete phase gradient. Two-dimensional gratings have also been constructed and used experimentally for phase imaging [11,12]. We propose here a regularization which avoids line artifacts. In contrast with the other methods, the solution proposed here is algorithmic, requiring the measurement of the differential phase only in one direction, and making use of statistical assumptions to compensate for the missing information.

Although addressed to physics-oriented readers, the presented solution uses notions of linear algebra. Textbooks like [13] may provide the reader more information on concepts used here.

2. Formulation

We consider $\Phi$ and $\Phi _u$ as matrices with sizes $m \times n$ and $m \times (n - 1)$ respectively. The differential direction $\vec {u}$ is assumed horizontal, along the second coordinate of the matrix. The observation model for the phase shift can be written as follows:

The proposed approach to improve the estimation of $\Phi$ is by including prior information about its statistics. Because measurements of the derivative are only in one direction, we assume a statistical distribution for the unmeasured vertical derivative $F_m \Phi$, along the first coordinate of the $\Phi$ matrix. In reality this distribution really depends on the sample that is measured and shall be adapted to the application. For the imaging of natural samples, we propose a simple model where the unmeasured vertical differential phase follows a centered Gaussian distribution of variance ${\sigma }_p^2$. This decreases the likelihood of strong line artifacts present in phase images resulting from naive integration of noisy differential phase. More realistic models of prior distribution can be considered for specific applications, although the phase estimation will then require numerical optimization methods. With our simple Gaussian model for the vertical differential phase distribution, the prior density of the phase is:

Using Bayes law, the posterior distribution for the phase is proportional to the product of Gaussian densities $P \left (\Phi _u | \Phi \right )$ and $P \left ( \Phi \right )$, so that the maximum a posteriori estimate is obtained as

This minimization problem has the form of a Tikhonov regularization [16]. Derivating the minimized term in Eq. (5) and equating to zero gives the following condition on solution $\hat {\Phi }$:

This last equation in $\hat {\Phi }$ is known as the Sylvester equation [17]. An analytic solution can be found by writing it with Kronecker products and using the following identity for any matrix product $PQR$:

This can be clearly identified as its eigendecomposition $V_{M} D_{M} V_{M}^\textrm {T}$, from which the pseudoinverse is $M^{+} = V_{M} D_{M}^{+} V_{M}^\textrm {T}$. The estimate $\hat {\Phi }$ becomes

Using Eq. (7) and the fact that $D_M$ is diagonal, we obtain the final expression of $\hat {\Phi }$, that is used for its actual computation:

As stated before, the obtained estimate is the minimum norm solution of Eq. (5). It can be easily shown that the obtained $\hat {\Phi }$ has therefore a zero mean value. Although this mean value is unrealistic, an estimate of the absolute phase shift may not be retrieved without additional prior information and the use of a more complex iterative algorithm.

3. Simulation

To quantitatively assess the performance of the proposed phase estimation, we have performed simulations using the Barbara and Pepper images classically used in the image processing community, as well as the Shepp Logan phantom (seen in Fig. 1). Measurements are simulated by taking the horizontal derivative of each image and adding independent Gaussian noise with a known variance. The definition of signal to noise ratio (SNR) used in simulated noisy images is the ratio between the noiseless image dynamic range and the additional random noise standard deviation. An estimation of the original image is then performed in the two following ways. First, a basic integration is done by filtering in reciprocal space by $1/(2 \pi i \omega _u)$. Second, the proposed method is applied, using the prior knowledge of signal and noise variances from the simulated data (${\sigma }_p^2$ and ${\sigma }_o^2$). Results are compared by the peak signal to noise ratio (PSNR) between estimated and original images :

 

Fig. 1. Example of simulated measurements and estimation by the two methods. The initial phase image is 512 $\times$ 512 pixels, with a dynamic range of 256. (a) Horizontal derivative of the modified Shepp-Logan phantom, with additive Gaussian noise ($\textrm {SNR } = 100$). (b) Integrated image by filtering in reciprocal space, showing strong line artifacts. (c) Image estimated with the proposed regularization.

Download Full Size | PDF

 

Fig. 2. Evolution of the PSNR between original zero mean images and estimates from differential measurements in function of the measurement SNR, for the basic integration and the proposed regularized method.

Download Full Size | PDF

4. Application to the experimentally obtained image

The proposed regularized phase estimation has also been demonstrated on experimental data acquired by an X-ray grating interferometer with two gratings (Talbot interferometer [18]), using a white synchrotron X-ray source from a bending magnet source at BL28B2 in Spring-8, Japan. The experimental setup is basically the same as those reported in [19–21]. A gold $\pi /2$ phase grating for the X-ray energy of 25 keV and a gold absorption grating with a thickness of 40 $\mathrm{\mu} m$ was used in the interferometry. The pitch of the two gratings was 5.3 $\mathrm{\mu} m$. Their lines are alined in the horizontal direction to take advantage of high-spatial coherence of the beam in the vertical direction. The distance between the two gratings was set to 283 mm, which corresponds to a Talbot order of 0.5 for the X-ray energy of 25 keV. A sCMOS camera (Hamamatsu ORCA-Flash 4.0, $2048 \times 2048$ pixels) and a 10 $\mathrm{\mu} m$ thick Gadox scintillator screen were used. The effective pixel size was 6.7 $\mathrm{\mu} m$. A 5-step fringe scanning [22] were performed to obtain a differential phase image with an exposure time of 1 ms for each step.

The imaged specimen was a cherry. A region of the obtained differential phase image is presented in Fig. 3(a), where the differential direction is taken to be horizontal. Note that, while a white synchrotron X-ray beam was used here, the finite band width does not affect the quantitativeness of the obtained differential phase from the criterion reported in [ [20]]. Figures 3(b) and (c) are phase images obtained from Fig. 3 after the subtraction of a smooth background. The simple integral along the horizontal axis, Fig. 3(b), was actually computed by dividing the Fourier transform of $\Phi _u$ by $2 \pi i \omega _u$, with $\omega _u$ the horizontal spatial frequency coordinate. Typical strong line artifacts are seen. The proposed regularized phase estimate $\Phi$. in Fig. 3(c) successfully avoids such artifacts. It was obtained using the formula in Eq. (13). ${\sigma }_\textrm {o}^{2}$ was measured from the differential phase measurements in a region with no sample. ${\sigma }_\textrm {p}^{2}$ was heuristically assumed equal to the variance of $\Phi _u$ minus ${\sigma }_\textrm {o}^{2}$. Figure 3(e) shows the phase image obtained by the proposed regularization method with vascular bundles clearly visualized by a further background correction. As references, the background-corrected absorption and normalized-visibility images obtained from the same experiment are displayed in Fig. 3(d) and (f). The three images show similar structures but different details in soft tissue regions.

 

Fig. 3. Images of a $1,100 \times 900$-pixel region of a cherry obtained by X-ray Talbot interferometry. (a) Differential phase. (b) Phase estimation obtained by integration performed in the reciprocal space. (c) Result of the proposed regularized phase estimation. (d-f) Comparison of Absorption (d), phase obtained by the proposed regularization method (e), and normalized visibility (f) images with vascular bundles more visualized.

Download Full Size | PDF

5. Summary

In summary, we have proposed a regularized estimation of the phase shift from the one-dimensional beam deflection measurements in X-ray grating interferometry. This makes use of an assumption on the distribution of phase gradient values, as well as the measurement noise. Using eigendecomposition, the maximum a posteriori solution is given in a tractable form. Results on simulated and experimental data have shown a strong improvement of the phase estimation compared to the simple integral method.

The proposed method is suitable for the grating-based X-ray tomography with the tomographic axis parallel to the direction of differentiation. When the tomographic axis is perpendicular to the grating lines, tomograms can be directly obtained from differential phase images [23], but when it is parallel to the grating lines, integration of the differential phase images are essential. The direction of the grating lines is determined by the spatial coherence and the horizontal direction is better for an X-ray Talbot interferometer using a white X-ray beam from a third generation synchrotron X-ray source [19–21]. For this reason, the integration is often desired for the cases where, for example, the sample has to be rotated around the vertical direction because it otherwise deforms by gravity. Another example is the multi-beam X-ray imaging system recently proposed where the tomographic axis is also in the vertical direction [24,25]. Furthermore, visibility tomography for a sample with anisotropic small-angle X-ray scattering requires the integration for X-ray phase tomography [7,26]. Note that the proposed method can also be applied to other imaging techniques based on beam deflection measurements including the X-ray diffraction enhanced imaging (DEI) [27–30] and those using visible light. Thus, the proposed method can be widely used in optical beam deflection imaging techniques.