## Abstract

A method for the estimation of the phase shift from beam deflection measurements in one direction is proposed in the context of X-ray grating interferometry. The common approach, which consists in simply integrating the measurements along this direction, produces typical line artifacts due to noise and missing information. Therefore, an algorithmic method is proposed, based on independence of measurement noise to differential phase and on the prior knowledge that the unmeasured phase derivatives are normally distributed over the field of view. This is shown to be equivalent to a Tikhonov regularization by a difference operator. A computationally tractable formulation of the optimal solution is derived. The method is demonstrated on experimental data and quantitatively evaluated by numerical simulations.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 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:

where $F_n$ is the derivative operator defined here explicitely as the following $(n - 1) \times n$ matrix: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$:

where $\textrm {vec} (\cdot )$ denotes the vector made from the concatenation of columns in the matrix operand. Equation (6) is then written as withThis 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 :

## 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.

## 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.

## Funding

Japan Science and Technology Agency (JPMJCR1765); Japan Society for the Promotion of Science (21H04530).

## Acknowledgements

The experiment was performed in SPring-8. This research was funded by Japan Science and Technology Agency (JST) (Grant No. JPMJCR1765) and Japan Society for the Promotion of Science (JSPS) (Grant No. 21H04530).

## Disclosures

The authors declare no conflicts of interest.

## Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## References

**1. **C. David, B. Nöhammer, and H. H. Solak, “Differential x-ray phase contrast imaging using a shearing interferometer,” Appl. Phys. Lett. **81**(17), 3287–3289 (2002). [CrossRef]

**2. **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]

**3. **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]

**4. **A. Momose, W. Yashiro, Y. Takeda, Y. Suzuki, and T. Hattori, “Phase tomography by X-ray Talbot interferometry for biological imaging,” Jpn. J. Appl. Phys. **45**(6A), 5254–5262 (2006). [CrossRef]

**5. **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]

**6. **F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Hard-X-ray dark-field imaging using a grating interferometer,” Nat. Mater. **7**(2), 134–137 (2008). [CrossRef]

**7. **W. Yashiro, Y. Terui, K. Kawabara, and A. Momose, “On the origin of visibility contrast in x-ray Talbot interferometry,” Opt. Express **18**(16), 16890–16901 (2010). [CrossRef]

**8. **S. Lian and H. Kudo, “Phase unwrapping with differential phase image,” Proc. SPIE **10132**, 101325N (2017). [CrossRef]

**9. **C. Kottler, C. David, F. Pfeiffer, and O. Bunk, “A two-directional approach for grating based differential phase contrast imaging using hard x-rays,” Opt. Express **15**(3), 1175–1181 (2007). [CrossRef]

**10. **M. Jiang, C. L. Wyatt, and G. Wang, “X-ray phase-contrast imaging with three 2d gratings,” Int. J. Biomed. Imag. **2008**, 1–8 (2008). [CrossRef]

**11. **T. Weitkamp, I. Zanette, C. David, J. Baruchel, M. Bech, P. Bernard, H. Deyhle, T. Donath, J. Kenntner, S. Lang, J. Mohr, B. Müller, F. Pfeiffer, E. Reznikova, S. Rutishauser, G. Schulz, A. Tapfer, and J.-P. Valade, “Recent developments in x-ray Talbot interferometry at ESRF-ID19,” Proc. SPIE **7804**, 780406 (2010). [CrossRef]

**12. **I. Zanette, T. Weitkamp, T. Donath, S. Rutishauser, and C. David, “Two-Dimensional X-Ray Grating Interferometer,” Phys. Rev. Lett. **105**(24), 248102 (2010). [CrossRef]

**13. **J. S. Golan, * Foundations of Linear Algebra* (Springer, 1995).

**14. **W. Yashiro, Y. Takeda, and A. Momose, “Efficiency of capturing a phase image using cone-beam x-ray talbot interferometry,” J. Opt. Soc. Am. A **25**(8), 2025–2039 (2008). [CrossRef]

**15. **V. Revol, C. Kottler, R. Kaufmann, U. Straumann, and C. Urban, “Noise analysis of grating-based x-ray differential phase contrast imaging,” Rev. Sci. Instrum. **81**(7), 073709 (2010). [CrossRef]

**16. **A. N. Tikhonov and V. Y. Arsenin, * Solutions of ill-posed problems* (Winston, 1977).

**17. **R. H. Bartels and G. W. Stewart, “Solution of the matrix equation AX + XB = C [F4],” Commun. ACM **15**(9), 820–826 (1972). [CrossRef]

**18. **S. Yokozeki and T. Suzuki, “Shearing Interferometer Using the Grating as the Beam Splitter,” Appl. Opt. **10**(7), 1575–1580 (1971). [CrossRef]

**19. **W. Yashiro, D. Noda, and K. Kajiwara, “Sub-10-ms X-ray tomography using a grating interferometer,” Appl. Phys. Express **10**(5), 052501 (2017). [CrossRef]

**20. **W. Yashiro, R. Ueda, K. Kajiwara, D. Noda, and H. Kudo, “Millisecond-order X-ray phase tomography with compressed sensing,” Jpn. J. Appl. Phys. **56**(11), 112503 (2017). [CrossRef]

**21. **W. Yashiro, C. Kamezawa, D. Noda, and K. Kajiwara, “Millisecond-order X-ray phase tomography with a fringe-scanning method,” Appl. Phys. Express **11**(12), 122501 (2018). [CrossRef]

**22. **J. H. Bruning, D. R. Herriott, J. E. Gallagher, D. P. Rosenfeld, A. D. White, and D. J. Brangaccio, “Digital wavefront measuring interferometer for testing optical surfaces and lenses,” Appl. Opt. **13**(11), 2693–2703 (1974). [CrossRef]

**23. **G. W. Faris and R. L. Byer, “Three-dimensional beam-deflection optical tomography of a supersonic jet,” Appl. Opt. **27**(24), 5202–5212 (1988). [CrossRef]

**24. **W. Voegeli, K. Kajiwara, H. Kudo, T. Shirasawa, X. Liang, and W. Yashiro, “A Multi-beam X-ray optical system for high-speed tomography,” Optica **7**(5), 514–517 (2020). [CrossRef]

**25. **T. Shirasawa, X. Liang, W. Voegeli, E. Arakawa, K. Kajiwara, and W. Yashiro, “High-speed multi-beam X-ray imaging using a lens coupling detector system,” Appl. Phys. Express **13**(7), 077002 (2020). [CrossRef]

**26. **W. Yashiro, S. Harasse, K. Kawabata, H. Kuwabara, T. Yamazaki, and A. Momose, “Distribution of Unresolvable Anisotropic Microstructures Revealed in Visibility-Contrast Images using X-ray Talbot Interferometry,” Phys. Rev. B **84**(9), 094106 (2011). [CrossRef]

**27. **J. Davis, D. Gao, T. E. Gureyev, A. W. Stevenson, and S. W. Wilkins, “Phase-contrast imaging of weakly absorbing materials using hard X-rays,” Nature **373**(6515), 595–598 (1995). [CrossRef]

**28. **D. Chapman, W. Thomlinson, R. E. Johnston, D. Washburn, E. Pisano, N. Gmür, Z. Zhong, R. Menk, F. Arfelli, and D. Sayers, “Diffraction enhanced x-ray imaging,” Phys. Med. Biol. **42**(11), 2015–2025 (1997). [CrossRef]

**29. **A. Ando, H. Sugiyama, A. Maksimenko, W. Pattanasiriwisawa, K. Hyodo, and Z. Xiaowei, “A new optics for dark-field imaging in X-ray region ‘Owl.’,” Jpn. J. Appl. Phys **40**(Part 2, No. 8A), L844–L846 (2001). [CrossRef]

**30. **A. Yoneyama, T. T. Lwin, and M. Kawamoto, “Fast diffraction-enhanced imaging using continuous sample rotation and analyzer crystal scanning,” J Synchrotron Rad. **27**(2), 468–471 (2020). [CrossRef]