Phase retrieval based on deep learning in grating interferometer

Grating interferometry is a promising technique to obtain differential phase contrast images with illumination source of low intrinsic transverse coherence. However, retrieving the phase contrast image from the differential phase contrast image is difficult due to the accumulated noise and artifacts from the differential phase contrast image (DPCI) reconstruction. In this paper, we implemented a deep learning-based phase retrieval method to suppress these artifacts. Conventional deep learning based denoising requires noise/clean image pair, but it is not feasible to obtain sufficient number of clean images for grating interferometry. In this paper, we apply a recently developed neural network called Noise2Noise (N2N) that uses noise/noise image pairs for training. We obtained many DPCIs through combination of phase stepping images, and these were used as input/target pairs for N2N training. The application of the N2N network to simulated and measured DPCI showed that the phase contrast images were retrieved with strongly suppressed phase retrieval artifacts. These results can be used in grating interferometer applications which uses phase stepping method.


Results
The schematic of the proposed method is shown in Fig. 2. We obtained a set of phase stepping images with and without samples using a grating interferometer with a sinusoidal model of the signal shape. The numerical phantom was created for simulation, and we obtained phase stepping images, using the same signal model. Training of the N2N model requires many noisy images with the same expected values. We obtained the noisy DPCIs using a combination of the phase stepping images, and these images were used as network inputs and targets. The results of the Shepp-Logan phantom obtained from MATLAB are shown in Fig. 3. We compared the resulting image to the results obtained by direct integration, wavelet-Fourier (WF) filtering 14 , and an iterative method 32 . The peak-to-noise ratio (PSNR) and the structure similarity (SSIM) index are displayed in the results. There are several types of noises and artifacts in the DPCIs, including background noise and a moiré artifact due to interferometer alignment error. The direct integration result shows that these noise sources accumulate and produce very poor image quality. WF filtering can remove the horizontally accumulated artifacts, however, it is limited to removing only one directional noise, and hence the result is also of poor visual image quality. The iterative method minimizes a cost function to retrieve phase information, and it yields high quality of the retrieved phase  www.nature.com/scientificreports/ information except from horizontally blurry regions. This is because the DPCI we obtain is a one-directionally differentiated image. Therefore, the horizontal line was not retrieved well. The proposed method achieved better results compared to other methods, with even preserving the horizontal features. The network used in this study leveraged not only DPCIs, but also transmission and dark-field images. The transmission and dark-field images can help to preserve shapes of the phase information. The quantitative analysis showed that the result using the deep learning method is of the highest visual quality. We also used the "peppers image" from USC-SIPI image database which is known for many high frequency features and more complex shapes than the Shepp-Logan phantom. The results presented in Fig. 4 show that the phase information is well retrieved even in the case of a complex scene. The area in the blue box has features shown only in the transmission image, and that in the red box only in the PCI. The results show that the feature in the red box uniquely exiting in the PCI is well preserved in the proposed method while that in the blue box uniquely existing the transmission image has little effect the retrieved PCI. In these results, the N2N network can be employed on any grating interferometry dataset. For measured data, we employed images acquired from two separate X-ray grating interferometers and one neutron grating interferometer. The first X-ray grating interferometer is a conventional Talbot-Lau interferometer (TLI) in which G 1 is close to G 2 rather than G 0 3,15 . The second X-ray TLI employed the symmetric geometry in which G 1 is positioned to middle of the system 8,10,33 . A conventional TLI has the advantage to obtain high spatial resolution images, and symmetric TLI has the advantage to obtain highly sensitivity DPCI 34 . Training datasets were created through random combinations of half of the phase step images for each sample. Figure 5 shows the X-ray TLI results of a flower, a cicada, a leaf, with the first and second row obtained using the conventional geometry, and the third row is obtained from the symmetric geometry. Direct integration method results in very poor image quality, with no discernable detail in the PCI. The WF filtering technique also results in very low PCI quality. The iterative method yields a better PCI than the direct integration and WF filtering, where there are many more discernable features. However, the proposed method can provide significantly better estimates of the phase information compared to other retrieving methods. The resulting PCI has a similar shape as the transmission image, but it has different contrast and details. Therefore, it helps to obtain various information about the objects.
We also analyzed data from a neutron grating interferometer setup. The samples were coins and a step sample composed of brass and copper as shown in Figs. 6h and 7h. Usually, neutron images are noisier than X-ray images due to the low flux of neutron. Figures 6 and 7 compare the results of the various phase retrieval methods. Similar to the X-ray TLI results, the proposed method produces superior image quality and PCI estimates.
Lastly, we analyze X-ray TLI images of a circuit board. In this case, the circuit board is larger than the field-ofview (FOV). When analyzing an object that is larger than the TLI FOV, it is not possible to obtain the boundary condition of zero phase, and we don't set the air region, and the PCI may have low image quality. Figure 8 shows the retrieved PCI of the circuit board by the proposed method. The retrieved phase is of good quality even though mask images don't exist. The good quality is attributed to the fact that transmission and the dark-field images used as the inputs of the network have similar shapes to the phase image.

Discussion
The phase retrieval in grating interferometer is one of the important parts for material analysis. However, its use is limited for many reasons, such as motor control error, non-uniform phase coefficient, etc. which often generate severe image artifacts. In this study, we implemented the deep learning method to retrieve the phase contrast image produced by grating interferometers. We showed that this method has several advantages compared other  www.nature.com/scientificreports/ commonly employed phase retrieval methods. First, one can obtain a high-quality PCI with fast acquisition time if the network is previously trained. Second, the N2N deep learning network does not require pseudo-images for training since the N2N does not require clean images and can operate on real data directly. The proposed method can use noisy DPCIs using combination of the phase stepping images, and it is possible to train networks which are fitted for real data. However, there are some limitations to this method. First, the noise and artifacts in the images are not perfectly removed, and some percentage of them still remain in the images. In the published paper 31 , the condition of using N2N is that the dataset is sufficiently large and the conditional expectation between paired noisy images is zero. It means that the noises between each image pair should be zero-mean and independent of each other. In practice, the images used in this paper do not satisfy this condition. However, this condition can be approximately satisfied when the phase stepping number is very large. But due to the radiation dose problem, the value of "very large" and the amount of allowed radiation form a trade-off relationship. Second, we should set mask image before training. The mask images in this paper were obtained manually using threshold to the images. Nonetheless, the manual threshold method can be replaced by advanced and automatic methods for creating mask images.

Conclusion
The grating interferometer can obtain the information about the objects through various types of the images, and it is important to analyze the information in the images using proper methods. In this study, we implemented a deep learning-based phase retrieval method for grating interferometry. This method can be used both with simulation and measured data; it produces better image quality than conventional methods and has the advantage of being applicable to any grating interferometer employing the phase stepping method regardless of the radiation sources such as X-ray and neutron. We believe that this method will be an effective tool to retrieve the PCI from the DPCI, and it will improve the analysis of complex objects.

Methods
Talbot-Lau interferometer. The regular TLI uses three gratings: a source grating G 0 , a phase grating G 1 , an analyzer grating G 2 . The TLI uses diffraction that occurs when X-ray beam passes through the G 1 , and this diffraction pattern is repeated with certain distance by Talbot effect. However, the intensity pattern is usually smaller than the detector pixel size, the G 2 grating is used to make moiré fringe. This pattern is periodically   www.nature.com/scientificreports/ where I ref , I sam denote the intensity without and with samples, and a 0 , a 1, θ denote DC offset, amplitude, and the phase coefficient, respectively. M is a total phase stepping number, and k = 1, 2, …, M. η k is a mechanical error term that occurs the grating is moved to perform phase stepping. And we used least-square method 35 for sine curve fitting shown as where I and I denote the fitted images and the measured images, respectively. The transmission, DPCI, dark-field images can be calculated from the fitted parameters as follows: Phase retrieval method. While the transmission image and the dark-field image can be obtained directly from the phase stepping images, however, the PCI we can obtain is one-dimensionally differentiated. The mathematical expression of phase retrieval is shown as where φ and ϕ denote a PCI and a DPCI, respectively, and c is a constant given specification of grating interferometer system, n denotes the noise. If the DPCI is noise-free, the PCI is retrieved using simple one-dimensional integration and there are no noises in the PCI. However, in most cases, there are moiré artifacts due to the motion of the grating and noises in DPCIs, and these are accumulated during integration. Therefore, we do not use onedimensional integration directly, and need to optimize the images using other methods.
where D denotes a one-dimensional first-order derivative operator, and D x is a derivative operator in horizontal direction. And w is the regularization term. This cost function can be reformulated as this unconstrained form 17 .
where TV is total variation operator of the images, and M is a mask function that denotes prior information for air region. The first term of regularization term is total variation (TV) term, which is for denoising in PCIs, and the second term is binary term, which is for DC offset calibration. The selection of these parameters is critical for image quality. In this study, we manually set the mask images using threshold to the transmission, DPCI, dark-field images. And we used alternating direction method of multipliers (ADMM) for phase retrieval with 1 = 0.001 and 2 = 1, respectively. In this paper, we changed the cost-function to the deep learning form: where f θ x i denotes the deep learning network outputs. We used the following cost-function for phase retrieval. We used the same parameters to Eq. (10).
Hardware setup. X-ray grating interferometer. We installed two types X-ray TLI setups in a laboratory.
In case of the first setup, we used an open-type rotating anode X-ray tube with mean energy 28.0 keV, and an energy-integrating flat-panel type detector with pixel size of 49.6 μm. And all gratings have period of 6.0 μm. The all distance between three gratings were set to 610.0 mm. The second X-ray grating interferometer setup includes a closed-type X-ray tube with mean energy 23.0 keV, and a flat-panel detector which has same specification as the first setup. The G 0 , G 1 , G 2 gratings have period of 10.0 μm, 3.2 μm, 4.8 μm, respectively. The distance between G 0 and G 1 is 300 mm and the distance between G 1 and G 2 is 144.0 mm.
Scientific Reports | (2022) 12:6739 | https://doi.org/10.1038/s41598-022-10551-y www.nature.com/scientificreports/ Neutron grating interferometer. Neutron Talbot-Lau grating interferometer was installed at the cold neutron imaging beam line NG6 at NIST Center for Neutron Research (NCNR) 36 . The wavelength of the source is 0.44 nm with polychromatic beam. The detector is an Andor, a scientific complementary metal-oxide semiconductor (sCMOS) camera which has 50 mm lens, 2160 × 2560 pixels, cooling temperature of -30 °C, and 200 um LiF:ZnS scintillator with the effective pixel size of 51.35 μm (Certain trade names and company products are mentioned in the text or identified in an illustration in order to adequately specify the experimental procedure and equipment used. In no case does such identification imply recommendation or endorsement by the National Institute of Standards and Technology nor does it imply that the products are necessarily the best available for the purpose). The G 0 and G 1 , G 2 gratings have same period of 50 μm. The distances between all gratings were set to 4260 mm. Each phase stepping image was obtained for 15 s.

Dataset preparation.
We obtain a set of the phase stepping images to obtain information: phase contrast, transmission, dark-field images. The number of phase stepping can be varied and should be optimized for parameters such as image quality, data acquisition time, noise and so on. Given that we have a certain number of phase stepping, in conventional cases, all phase stepping images are utilized for generating the DPCI, but a DPCI can also be calculated by combination of the phase stepping images. In this way, we can generate a number of DPCIs that have statistical variations, and these images can be used for N2N. In this paper, the transmission and the dark-field images were generated in a conventional way, and the noisy DPCI image pairs were generated by using a combination of the phase stepping images per each iteration and used for training process. where Image is the normalized images. And we added Poisson noise in the entire of the images. The phase stepping numbers are 8, 12, respectively. We used least-square method for data fitting. The obtained images were randomly flipped for data argumentation.
Network training. We designed a deep learning network in the form of a dual stream network shown in Fig. 9 37 . The dual steam network consists of base/detail layers for denoising. The role of the base layer is to prevent image quality degradation by making base of the result image. In this study, we used the transmission and the dark-field images as the base layer input because they can be obtained using grating interferometer, are not differentiated images, and have a similar shape to the PCI. And we used the DPCI as an input of the detail layer. The role of the detail layer is to obtain detailed features in the result image. Each stream has U-net type network with depth 2, kernel size 3, and rectified linear unit (ReLU) 23 . The features have initially 64 features, and it is doubled as the network is deeper. The features from each stream are added and connected to three convolutional blocks. And the result image can be obtained from the last convolutional layer. The cost function in this study calculates combined losses which are L 2 loss between differentiated images shown and TV loss, mask image loss shown in Eq. (11). During the training process, the ADAM algorithm was used as an optimizer 38 . The initial learning rate was set to 0.0001 and it decayed √ t times with t-th epochs. The network was trained for 100 epochs, and 100 times iterated in one epoch. For every iteration, the reference and the sample images were randomly selected from the half of the total phase stepping number to the total phase stepping number and extracted 64 batches with pixel size of 64 × 64. Network training was performed on a Tensorflow deep learning framework 39 .
Ethics declaration. The authors declare no humans were directly used in this study.

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.