Direct information retrieval after 3D reconstruction in grating-based X-ray phase-contrast computed tomography

The theoretical framework of the direct information retrieval method is presented in phase-contrast computed tomography. Numerical simulations are also performed, which reveal that the proposed method provides comparable results of the reverse projection method with short computational time. Furthermore, the compatibility with the existing data preprocessing methods and iterative reconstruction algorithms is discussed.


Introduction
X-ray phase-contrast imaging provides information about an object that would be otherwise inaccessible (Nugent et al., 1996;Fitzgerald, 2000) using conventional attenuation-based X-ray imaging. Along with improved X-ray sources and optical components over the last two decades, X-ray phasecontrast imaging has experienced a great evolution (Wilkins et al., 1996;Momose et al., 1996;Zhou & Brahme, 2008;Pfeiffer et al., 2008;Tapfer et al., 2011;Thuering et al., 2011;Bravin et al., 2013;Wu et al., 2015), especially with the introduction of grating-based X-ray phase-contrast imaging (GBPCI), which included Talbot interferometers (Momose et al., 2003;Weitkamp et al., 2005;Xi et al., 2012), Talbot-Lau interferometers (Pfeiffer et al., 2006) and non-interferometric grating-based imaging (Huang et al., 2009). Owing to its compatibility with conventional X-ray sources and the possibility of a large field of view, GBPCI has been considered as a potential imaging method in clinical medicine. In addition, phase-contrast computed tomography (PCCT) provides three-dimensional phase-contrast visualization from the retrieved two-dimensional refraction angle images at numerous tomographic viewing angles using a reconstruction algorithm. Generally, two-dimensional information retrieval is required prior to independent and discriminating absorption and phase threedimensional reconstruction. The existing data pre-processing methods and computed tomography (CT) algorithms in absorption reconstruction cannot be employed directly in phase reconstruction.
Recently, many highly efficient and ingenious information retrieval techniques (Zanette et al., 2012;Wang et al., 2014;Zhu et al., 2010;Yang et al., 2016) have been put forward to improve the performance of grating-based PCCT in terms of acquisition time and radiation dose. Following previous research in diffraction enhanced imaging Wang et al., 2007), Zhu et al. (2010) presented a fast and lowdose information retrieval method known as the reverse projection (RP) method. The phase information collapse phenomenon (Raupach & Flohr, 2011) can be avoided in this method; however, information retrieval still needs to be performed before three-dimensional reconstruction. Afterwards, Diemoz et al. (2011) proposed a simplified method, which obtained three-dimensional information without information retrieval. However, the displayed three-dimensional mixed image, containing both absorption and phase information, is not quantitative and is unable to take full advantage of the high contrast of phase-contrast images. Brendel et al. (2016) achieved three types of information tomography using iterative reconstruction without information retrieval.
Inspired by these methods, we present a direct threedimensional information retrieval method, where threedimensional reconstruction can be achieved using an absorption CT algorithm before information retrieval. Therefore, the existing data pre-processing methods and iterative reconstruction algorithms in absorption reconstruction may be introduced to phase reconstruction immediately. In addition, the computational time of the information retrieval procedure is shortened, owing to the reduced data size and the absence of the logarithm operation. Theoretical derivations and numerical simulations have been carried out to confirm its feasibility and validity.
2. Grating-based phase-contrast imaging setup and three-dimensional information retrieval method 2.1. Grating-based phase-contrast imaging setup For simplicity, we use a parallel-beam Talbot interferometer, for example, to interpret and certify the direct three-dimensional information retrieval method. As described in Fig. 1(a), a Talbot interferometer consists of a phase-shift grating (G1), an absorption grating (G2) and a detector. Here, the system rotation axis is vertical with respect to the grating lines so that the physical quantities are rotation invariants (Zhu et al., 2005). The Cartesian coordinate system (x, y, z) and (X, Y, z) are, respectively, the reference frame of the sample and the imaging system. The detected photon number at pixel (X, z), where X and z are integers, can be expressed as follows with negligible scattering of the object (Zhu et al., 2010): where z g is the displacement along the direction perpendicular to the beam axis and the grating lines between G1 and G2, ' is the tomographic viewing angle and I 0 is the detected photon number without the gratings and sample in the light path. The function S(z g ) describes the normalized background shifting curve (SC). M represents the attenuation term, which can be formulated by the integral of the absorption coefficient of the object along the light path: research papers the fractional Talbot distance D is the distance between G1 and G2, and is the refraction angle induced by the object, associated with the real-part decrement of the refractive index by The normalized background SC S(z g ) shown in Fig. 1(b) can be plotted by recording the photon number with the increase of z g when the object is absent in the light path.

Three-dimensional information retrieval method
Carrying out the logarithmic operation on both sides of equation (1) yields: It is worth noting that the logarithm of the SC can be well approximated to the first-order Taylor expansion around upslope and down-slope positions indicated in Fig. 1 ; the axial symmetry of the SC, has been considered in the derivation. Substituting equations (2), (3) and (5) into equation (4) gives ln S p 2 4 À ln I X; z; AEðp 2 =4Þ; ' Â Ã The hybrid three-dimensional data of the absorption and phase information can be reconstructed by absorption CT algorithms. A conventional filtered backprojection (FBP) method  is employed for reconstruction in equation (7): where F and j j stand for the Fourier transform and the absolute value of the frequency, respectively. Then, the three-dimensional absorption coefficient and phase-contrast data are obtained by: It can be found from equations (8) and (9) that the threedimensional absorption and phase-contrast images are retrieved directly from two three-dimensional images reconstructed by an absorption CT algorithm. Thus the twodimensional information retrieval procedures at every tomographic viewing angle are totally abandoned and a direct three-dimensional information retrieval method can be accessible in theory. Furthermore, we predict that this method has the potential to make a direct introduction of the established data pre-processing techniques and other advanced CT algorithms in absorption reconstruction to phase reconstruction.

Numerical simulations
Numerical simulations have been performed to validate the feasibility and veracity of the proposed method. Perfect gratings and alignment are assumed in the simulations. The energy was designed at 25 keV and the pitches of the phase shift (/2) grating and analyzer grating were 6 mm. We placed G2 at the first fractional Talbot distance where D = 0.3629 m downstream from G1 along the beam axis. The detector with pixel size of 100 mm was located closely after the grating G2. The visibility of the background SC was fixed at 0.3, referring to the reported value in a laboratory experiment (Zhu et al., 2010). The incoming photon number I 0 was 10000 per pixel and Poisson noise was introduced. The simulated phantom shown in Fig. 1(c) is a cube of polyethylene with 255 Â 255 Â 255 volume elements. The inner structure consists of a spherical shell of polycarbonate with an inner diameter of 0.64 cm and an external diameter of 1.28 cm. The attenuation coefficients and the real-part decrements of the refractive indexes of the materials are listed in Table 1 (refer to the CSIRO website, https://www.ts-imaging.net/Services/Simple/ ICUtilXdata.aspx).
In the simulations, two sets of 360 projection images have been collected evenly around a circle, when the grating G2 was positioned at the up-slope and down-slope positions. We obtained three-dimensional absorption and phase-contrast images by the RP method and the proposed method. In research papers Table 1 Attenuation coefficients and real part decrements of polyethylene and polycarbonate. Fig. 1 the RP method, the two-dimensional phase and absorption information were separated first and then reconstructed with different CT algorithms independently. In the proposed method, two hybrid three-dimensional images were reconstructed by an absorption CT algorithm and then the threedimensional information was separated using equations (8) and (9). The coronal slices in the black-grid plane depicted in Fig. 1(c) of the reconstructed mixed images at the up-slope and down-slope positions are illustrated in Figs. 2(a) and 2(b). We can see that they are almost identical with the exception of the boundary of the annulus where the phase signal is nonzero. This is because both images contain the same attenuation coefficient and the opposite derivative of the real-part decrement of the refractive index as revealed in equation (6). Using the information retrieval formulae (8) and (9), the absorption and phase-contrast images in Figs. 2(c) and 2(d) can be retrieved directly.

Colour in
As illustrated in Figs. 3(a)-3(b), the reconstructed absorption and phase-contrast images by the RP method are almost identical to their counterparts obtained by the direct retrieval method. In Fig. 3(c), we present the profiles of the theoreticaland retrieved-phase information, obtained by the RP method and the direct retrieval method, respectively, at the white dashed line plotted in Fig. 3(b). The good agreement between the two methods confirms the feasibility and veracity of the proposed method. The similar noise variances of 50 Â 50 pixels in the white square in Fig. 3(b) (1.81 Â 10 À9 and 1.83 Â 10 À9 ) confirm the equivalence of the direct retrieval method and the RP method.
The simulations were repeated five times to compare the computational complexity between the RP method and the proposed method. The computational time is recorded in Table 2, where DIR and RP stand for the presented direct information method and the reverse projection method, respectively. It reveals that the time for the information retrieval is significantly shortened by the proposed method. This may be explained by the reduced data size, from 360 Â 363 Â 255 pixels in the RP method to 255 Â 255 Â 255 pixels in the DIR method, and the absent logarithm operation. It is worth noting that the total computing efficiency becomes significant when the ratio between the tomographic viewing angle and the number of pixels in the width increases. Moreover, the computational time of CT reconstruction executed on graphic processing units (GPUs) or field-programmable gate arrays (FPGAs) will be sharply reduced, similar to the speedups of 26.9 for FPGAs in Choi & Cong (2016)  Results of the direct three-dimensional information retrieval method; coronal slices at the black grid plane are depicted in Fig. 1(c). Reconstructed images at the up-slope (a) and down-slope (b); retrieved absorption image (c) and differential phase-contrast image (d).
the computational time of CT reconstruction is speeded up 20 times, we can predict that more than one-third of the total computing time of the RP method will be saved by employing the DIR method in medical 16-slice CT data with a size of 2880 Â 816 Â 16 pixels.

Discussions
In this paper, we derive the direct three-dimensional information retrieval method in parallel beam geometry, which can also be employed in the case of fan beam geometry. Taking equiangular fan beam geometry as an example, we introduce a polar coordinate (r, , z) to describe the imaging system. When the sample is placed downstream of the grating G1, and the distance between sample and detector is R, the detected photon number is (Wu et al., 2013): The absorption term and refraction angle can be expressed as Further derivations can be completed by referring to equations (4)-(9), replacing equation (7) with a fan beam FBP reconstruction algorithm (Wu, 2017). One can also complete the derivations for the case where the sample is in front of the grating G1, referring to our previous paper (Wu et al., 2013). It should be noted that this method is not suitable for cone beam geometry because of the change in the direction of partial derivatives of at the plane z 6 ¼ 0.
Compared with the previous method, the proposed method has two advantages: computational efficiency, which has been discussed in x3; and compatibility with the existing data preprocessing methods, such as the ring correction method (Kim et al., 2014) and iterative reconstruction algorithms (Arcadu et al., 2017;Wang et al., 2017) in conventional absorption-based computed tomography. In the line-ratio ring artefacts correction method, the sensitivity difference between adjacent detector elements can be extracted from the ratios of the sum at the up-slope and down-slope of adjacent detectors along all   Table 2 Comparison of the computational time between the three-dimensional direct information retrieval method and the reverse projection method. projection angles. In contrast, the information is derived from values at two different detector positions in the scanning mode mentioned by Zhu et al. (2010), so the line-ratio ring artefacts correction method is invalid in that case. Furthermore, taking the physical quantity of x; y; z ð ÞÇC@ x; y; z ð Þ=@z as a whole, the absorption-based iterative reconstruction algorithms are directly generalized to phase-contrast tomography due to the rotation invariance. Hence, interior tomography and shortscan tomography in phase-contrast imaging may be handled immediately.
We would like to point out several limitations of the proposed method: (1) Perfect gratings and alignment are required, which limits the field of view or reduces imaging performance. With the development of grating fabrication and large periodic gratingbased non-interferometric imaging (Huang et al., 2009), these problems may be mitigated.
(2) The proposed direct information retrieval method is only discussed for the case where the object is free of the scattering information. As a matter of fact, if the object possesses high scattering information but weak phase information, a direct three-dimensional information retrieval method can also be carried out at peak and valley positions for attenuation and scattering information.
(3) The retrieved phase information is the derivative of the real-part decrement of the refractive index. The real-part decrement of the refractive index free of stripe artefacts can be recovered with an iterative method proposed by Thü ring et al. (2011).
(4) A two-step scanning mode is required in the presented method; however, there is only one step in the RP method. Recently, a newly designed grating was introduced in an applicable imaging method without mechanical phase stepping . The proposed method can also realise single-shot imaging with this grating.

Conclusions
In summary, we have presented a direct three-dimensional information retrieval method in grating-based X-ray differential phase-contrast computed tomography. This method obtains the absorption coefficient and the derivative of the real-part decrement of the refractive index immediately from two three-dimensional reconstructed images. Compared with the previous method, the proposed method improves the computational efficiency with comparative image quality and introduces the existing data pre-processing methods and the iterative reconstruction algorithms in absorption reconstruction to phase reconstruction directly. The presented method also lends itself to other X-ray differential phase-contrast imaging, such as analyzer-based imaging and edge illumination.

Funding information
The following funding is acknowledged: the China Postdoctoral Science Foundation (grant No. 2017M612097 to Zhao