research papers
Direct information retrieval after 3D reconstruction in gratingbased Xray phasecontrast computed tomography
^{a}National Synchrotron Radiation Laboratory, University of Science and Technology of China, Hefei, Anhui 230029, People's Republic of China, ^{b}School of Electronic Science and Applied Physics, Hefei University of Technology, Hefei, Anhui 230009, People's Republic of China, and ^{c}Institute of HighEnergy Physics, Chinese Academy of Sciences, Beijing 100049, People's Republic of China
^{*}Correspondence email: wuzhao@ustc.edu.cn, ychtian@ustc.edu.cn
Gratingbased Xray differential phasecontrast imaging has attracted a great amount of attention and has been considered as a potential imaging method in clinical medicine because of its compatibility with the traditional Xray tube source and the possibility of a large field of view. Moreover, phasecontrast computed tomography provides threedimensional phasecontrast visualization. Generally, twodimensional information retrieval performed on every projection is required prior to threedimensional reconstruction in phasecontrast computed tomography. In this paper, a threedimensional information retrieval method to separate absorption and phase information directly from two reconstructed images is derived. Theoretical derivations together with numerical simulations have been performed to confirm the feasibility and veracity of the proposed method. The advantages and limitations compared with the reverse projection method are also discussed. Owing to the reduced data size and the absence of a logarithm operation, the computational time for information retrieval is shortened by the proposed method. In addition, the hybrid threedimensional images of absorption and phase information were reconstructed using an absorption reconstruction algorithm, hence the existing data preprocessing methods and iterative reconstruction algorithms in absorption reconstruction may be utilized in phase reconstruction immediately.
Keywords: direct information retrieval method; gratingbased phasecontrast imaging; differential phasecontrast imaging; computed tomography.
1. Introduction
Xray phasecontrast imaging provides information about an object that would be otherwise inaccessible (Nugent et al., 1996; Fitzgerald, 2000) using conventional attenuationbased Xray imaging. Along with improved Xray sources and optical components over the last two decades, Xray 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 gratingbased Xray phasecontrast 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 noninterferometric gratingbased imaging (Huang et al., 2009). Owing to its compatibility with conventional Xray sources and the possibility of a large field of view, GBPCI has been considered as a potential imaging method in clinical medicine. In addition, phasecontrast computed tomography (PCCT) provides threedimensional phasecontrast visualization from the retrieved twodimensional refraction angle images at numerous tomographic viewing angles using a reconstruction algorithm. Generally, twodimensional information retrieval is required prior to independent and discriminating absorption and phase threedimensional reconstruction. The existing data preprocessing 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 gratingbased PCCT in terms of acquisition time and radiation dose. Following previous research in diffraction enhanced imaging (Zhu et al., 2006; 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 threedimensional reconstruction. Afterwards, Diemoz et al. (2011) proposed a simplified method, which obtained threedimensional information without information retrieval. However, the displayed threedimensional mixed image, containing both absorption and phase information, is not quantitative and is unable to take full advantage of the high contrast of phasecontrast 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 preprocessing 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. Gratingbased phasecontrast imaging setup and threedimensional information retrieval method
2.1. Gratingbased phasecontrast imaging setup
For simplicity, we use a parallelbeam Talbot interferometer, for example, to interpret and certify the direct threedimensional information retrieval method. As described in Fig. 1(a), a Talbot interferometer consists of a phaseshift 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 of the object along the light path:
the fractional Talbot distance D is the distance between G1 and G2, and is the refraction angle induced by the object, associated with the realpart decrement δ of the 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.
2.2. Threedimensional 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 firstorder Taylor expansion around upslope and downslope positions indicated in Fig. 1(b),
where
the axial symmetry of the SC, has been considered in the derivation. Substituting equations (2), (3) and (5) into equation (4) gives
The hybrid threedimensional data of the absorption and phase information can be reconstructed by absorption CT algorithms. A conventional filtered backprojection (FBP) method (Huang et al., 2006) is employed for reconstruction in equation (7):
where F and stand for the Fourier transform and the absolute value of the frequency, respectively. Then, the threedimensional and phasecontrast data are obtained by:
It can be found from equations (8) and (9) that the threedimensional absorption and phasecontrast images are retrieved directly from two threedimensional images reconstructed by an absorption CT algorithm. Thus the twodimensional information retrieval procedures at every tomographic viewing angle are totally abandoned and a direct threedimensional 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 preprocessing techniques and other advanced CT algorithms in absorption reconstruction to phase reconstruction.
3. 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 µm. 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 µm 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 realpart decrements of the refractive indexes of the materials are listed in Table 1 (refer to the CSIRO website, https://www.tsimaging.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 upslope and downslope positions. We obtained threedimensional absorption and phasecontrast images by the RP method and the proposed method. In the RP method, the twodimensional phase and absorption information were separated first and then reconstructed with different CT algorithms independently. In the proposed method, two hybrid threedimensional 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 blackgrid plane depicted in Fig. 1(c) of the reconstructed mixed images at the upslope and downslope 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 and the opposite derivative of the realpart decrement of the as revealed in equation (6). Using the information retrieval formulae (8) and (9), the absorption and phasecontrast images in Figs. 2(c) and 2(d) can be retrieved directly.
As illustrated in Figs. 3(a)–3(b), the reconstructed absorption and phasecontrast 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 theoretical and retrievedphase 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 fieldprogrammable gate arrays (FPGAs) will be sharply reduced, similar to the speedups of 26.9 for FPGAs in Choi & Cong (2016) and 22 for GPUs and 35.8 for FPGAs in Després & Jia (2017). Assuming the computational time of CT reconstruction is speeded up 20 times, we can predict that more than onethird of the total computing time of the RP method will be saved by employing the DIR method in medical 16slice CT data with a size of 2880 × 816 × 16 pixels.

4. Discussions
In this paper, we derive the direct threedimensional 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 ≠ 0.
Compared with the previous method, the proposed method has two advantages: computational efficiency, which has been discussed in §3; 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 absorptionbased computed tomography. In the lineratio ring artefacts correction method, the sensitivity difference between adjacent detector elements can be extracted from the ratios of the sum at the upslope and downslope of adjacent detectors along all 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 lineratio ring artefacts correction method is invalid in that case. Furthermore, taking the physical quantity of as a whole, the absorptionbased iterative reconstruction algorithms are directly generalized to phasecontrast tomography due to the rotation invariance. Hence, interior tomography and shortscan tomography in phasecontrast 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 noninterferometric 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 threedimensional 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 realpart decrement of the et al. (2011).
The realpart decrement of the free of stripe artefacts can be recovered with an iterative method proposed by Thüring(4) A twostep 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 (Wei et al., 2017). The proposed method can also realise singleshot imaging with this grating.
5. Conclusions
In summary, we have presented a direct threedimensional information retrieval method in gratingbased Xray differential phasecontrast computed tomography. This method obtains the
and the derivative of the realpart decrement of the immediately from two threedimensional reconstructed images. Compared with the previous method, the proposed method improves the computational efficiency with comparative image quality and introduces the existing data preprocessing methods and the iterative reconstruction algorithms in absorption reconstruction to phase reconstruction directly. The presented method also lends itself to other Xray differential phasecontrast imaging, such as analyzerbased imaging and edge illumination.Funding information
The following funding is acknowledged: the China Postdoctoral Science Foundation (grant No. 2017M612097 to Zhao Wu); the Fundamental Research Funds for the Central Universities (grant No. WK2310000065 to Zhao Wu); the National Research and Development Projects for Key Scientific Instruments (grant No. CZBZDYZ20140002).
References
Arcadu, F., Marone, F. & Stampanoni, M. (2017). J. Synchrotron Rad. 24, 205–219. CrossRef IUCr Journals Google Scholar
Bravin, A., Coan, P. & Suortti, P. (2013). Phys. Med. Biol. 58, R1–R35. Web of Science CrossRef PubMed Google Scholar
Brendel, B., von Teuffenbach, M., Noël, P. B., Pfeiffer, F. & Koehler, T. (2016). Med. Phys. 43, 188–194. CrossRef Google Scholar
Choi, Y. K. & Cong, J. (2016). IEEE Trans. Biomed. Circuits Syst. 10, 754–767. CrossRef Google Scholar
Després, P. & Jia, X. (2017). Phys. Med. 42, 76–92. Google Scholar
Diemoz, P. C., Coan, P., Zanette, I., Bravin, A., Lang, S., Glaser, C. & Weitkamp, T. (2011). Opt. Express, 19, 1691–1698. CrossRef Google Scholar
Fitzgerald, R. (2000). Phys. Today, 53, 23–26. Web of Science CrossRef Google Scholar
Huang, Z. F., Kang, K. J., Li, Z., Zhu, P. P., Yuan, Q. X., Huang, W. X., Wang, J. Y., Zhang, D. & Yu, A. M. (2006). Appl. Phys. Lett. 89, 041124. CrossRef Google Scholar
Huang, Z. F., Kang, K. J., Zhang, L., Chen, Z. Q., Ding, F., Wang, Z. T. & Fang, Q. G. (2009). Phys. Rev. A, 79, 013815. CrossRef Google Scholar
Kim, Y., Baek, J. & Hwang, D. (2014). Opt. Express, 22, 13380–13392. CrossRef Google Scholar
Momose, A., Kawamoto, S., Koyama, I., Hamaishi, Y., Takai, K. & Suzuki, Y. (2003). Jpn. J. Appl. Phys. 42, L866–L868. Web of Science CrossRef CAS Google Scholar
Momose, A., Takeda, T., Itai, Y. & Hirano, K. (1996). Nat. Med. 2, 473–475. CrossRef CAS PubMed Web of Science Google Scholar
Nugent, K. A., Gureyev, T. E., Cookson, D. F., Paganin, D. & Barnea, Z. (1996). Phys. Rev. Lett. 77, 2961–2964. CrossRef PubMed CAS Web of Science Google Scholar
Pfeiffer, F., Bech, M., Bunk, O., Kraft, P., Eikenberry, E. F., Brönnimann, C., Grünzweig, C. & David, C. (2008). Nat. Mater. 7, 134–137. Web of Science CrossRef PubMed CAS Google Scholar
Pfeiffer, F., Weitkamp, T., Bunk, O. & David, C. (2006). Nat. Phys. 2, 258–261. Web of Science CrossRef CAS Google Scholar
Raupach, R. & Flohr, T. G. (2011). Phys. Med. Biol. 56, 2219–2244. Web of Science CrossRef PubMed Google Scholar
Tapfer, A., Bech, M., Pauwels, B., Liu, X., Bruyndonckx, P., Sasov, A., Kenntner, J., Mohr, J., Walter, M., Schulz, J. & Pfeiffer, F. (2011). Med. Phys. 38, 5910–5915. Web of Science CrossRef PubMed Google Scholar
Thuering, T., Modregger, P., Grund, T., Kenntner, J., David, C. & Stampanoni, M. (2011). Appl. Phys. Lett. 99, 041111. CrossRef Google Scholar
Thüring, T., Modregger, P., Pinzer, B. R., Wang, Z. T. & Stampanoni, M. (2011). Opt. Express, 19, 25545–25558. Google Scholar
Wang, L., Guan, Y., Liang, Z.., Guo, L., Wei, C., Luo, R., Liu, G. & Tian, Y. (2017). J. Synchrotron Rad. 24, 490–497. CrossRef IUCr Journals Google Scholar
Wang, M., Zhu, P. P., Zhang, K., Hu, X. F., Huang, W. X., Cen, Y. W., Yuan, Q. X., Yu, X. L. & Wang, J. Y. (2007). J. Phys. D Appl. Phys. 40, 6917–6921. CrossRef Google Scholar
Wang, Z. L., Gao, K., Wang, D. J., Wu, Z., Chen, H., Wang, S. H. & Wu, Z. Y. (2014). Opt. Lett. 39, 877–879. CrossRef Google Scholar
Wei, C. X., Wu, Z., Wali, F., Luo, R. H. & Tian, Y. C. (2017). J. Phys. Conf. Ser. 849, 012036. CrossRef Google Scholar
Weitkamp, T., Diaz, A., David, C., Pfeiffer, F., Stampanoni, M., Cloetens, P. & Ziegler, E. (2005). Opt. Express, 13, 6296–6304. Web of Science CrossRef PubMed Google Scholar
Wilkins, S. W., Gureyev, T. E., Gao, D., Pogany, A. & Stevenson, A. W. (1996). Nature (London), 384, 335–338. CrossRef CAS Web of Science Google Scholar
Wu, Z. (2017). Med. Phys. 44, 2038. CrossRef Google Scholar
Wu, Z., Gao, K., Chen, J., Wang, D. J., Wang, S. H., Chen, H., Bao, Y., Shao, Q. G., Wang, Z. L., Zhang, K., Zhu, P. P. & Wu, Z. Y. (2015). Med. Phys. 42, 742–749. Google Scholar
Wu, Z., Gao, K., Wang, Z. L., Ge, X., Chen, J., Wang, D. J., Pan, Z. Y., Zhang, K., Zhu, P. P. & Wu, Z. Y. (2013). Med. Phys. 40, 031911. CrossRef Google Scholar
Xi, Y., Kou, B., Sun, H., Qi, J., Sun, J., Mohr, J., Börner, M., Zhao, J., Xu, L. X., Xiao, T. & Wang, Y. (2012). J. Synchrotron Rad. 19, 821–826. Web of Science CrossRef CAS IUCr Journals Google Scholar
Yang, Y., Xie, H. Q., Cai, W. X., Mao, H. & Tang, X. Y. (2016). Med. Phys. 43, 2855–2869. CrossRef Google Scholar
Zanette, I., Bech, M., Rack, A., Duc, G. L., Tafforeau, P., David, C., Mohr, J., Pfeiffer, F. & Weitkamp, T. (2012). Proc. Natl Acad. Sci. USA, 26, 10119–10204. Google Scholar
Zhou, S. & Brahme, A. (2008). Phys. Med. 24, 129–148. Web of Science CrossRef PubMed Google Scholar
Zhu, P. P., Wang, J. Y., Yuan, Q. X., Huang, W. X., Shu, H., Gao, B., Hu, T. D. & Wu, Z. Y. (2005). Appl. Phys. Lett. 87, 264101. CrossRef Google Scholar
Zhu, P. P., Yuan, Q. X., Huang, W. X., Wang, J. Y., Shu, H., Chen, B., Liu, Y. J., Li, E. R. & Wu, Z. Y. (2006). J. Phys. D Appl. Phys. 39, 4142–4147. CrossRef Google Scholar
Zhu, P. P., Zhang, K., Wang, Z. L., Liu, Y. J., Liu, X. S., Wu, Z. Y., McDonald, S. A., Marone, F. & Stampanoni, M. (2010). Proc. Natl Acad. Sci. USA, 107, 13576–13581. CrossRef Google Scholar
This is an openaccess article distributed under the terms of the Creative Commons Attribution (CCBY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.