## Abstract

We propose an algorithm for tomographic reconstruction of the refractive index map of an object translated across a fan-shaped X-ray beam. We adopt a forward image model valid under the non-paraxial condition, and use a unique mapping of the acquired projection images to reduce the computational cost. Even though the imaging setup affords only a limited angular coverage, our algorithm provides accurate refractive index values by employing the positivity and piecewise-smoothness constraints.

© 2013 OSA

## 1. Introduction

The highly penetrating nature of X-rays [1] has been utilized in a variety of applications, including non-destructive evaluation, medical imaging [2, 3] and baggage inspection [4]. The major source of contrast in the existing X-ray imaging techniques is derived from incoherent processes such as Compton scattering and photoelectric absorption. This type of absorption or attenuation contrast is sensitive to the atomic number of the sample under interrogation [5]. On the other hand, X-ray phase imaging (XPI) relies on the contrast due to non-homogeneous electron density distribution [6]. Therefore, it can differentiate soft tissues (e.g., tumor and healthy tissue) and detect improvised homemade explosives with high sensitivity. Several XPI techniques have been demonstrated using synchrotron [7–10] and, more recently on a lab scale, using micro-focus X-ray sources [11, 12] or traditional X-ray tubes [13–16]. For a recent review of XPI techniques, the readers are referred to [17].

X-ray phase tomography requires varying the angle of illumination that is incident on the specimen [18–21]. This may be achieved by rotating the imaging chain around the object, as in medical X-ray computed tomography, or by rotating the object in the case where the source is hard to move (for example, in synchrotron imaging). In either case, a rotational gantry or stage must be installed and the imaging chain must be carefully aligned with the object being imaged. For high-throughput applications such as luggage inspection at the airport, moving-object geometry is preferred [22, 23]. There are medical applications as well, e.g., breast tomosynthesis [24], where projections are available only from a linear track on one side of the object. Phase tomography of a moving object has been demonstrated in the ultrasound [25], terahertz [26], and optical regimes [27], but not in the X-ray regime. We recently reported a forward image model for X-ray phase imaging of thick and large objects [28]. It does not rely on the projection or paraxial approximation, and thus can be applied to the object at a far off-axis location. In this study, adopting the model, we propose an iterative algorithm for X-ray phase tomography of the object translated across an extended X-ray beam. We demonstrate that an artifact-free, depth-resolved refractive index map of a moving object can be obtained in an X-ray setup. We show that there exists an interesting analogy between our problem and plane-wave diffraction tomography (PW-DT).

## 2. Forward model for an object translated in a fan-beam X-ray system

Consider an object being translated across a fan-shaped X-ray beam in the geometry shown in Fig. 1. Suppose that the shape and material property of the object do not vary along the y-axis (i.e., perpendicular to the cross-sectional plane that is being imaged). In this 2-D model, we assume that a fan-shaped X-ray bean is emitted from an ideal line source. Let us place the source on the optical axis at *z* = -*R*_{1} and the 2-D object be translated along the line *z* = 0. The detector, located at *z* = *R*_{2}, records the intensity of transmitted X-ray beam after the object. In this setup, the detector is translated with the object so as to capture the diffracted beam from the object at varying locations. The detector shift *x _{i}* is determined by the magnification factor

*M*= (

*R*

_{2}+

*R*

_{1})⁄

*R*

_{1}and the object shift

*x*;

_{o}*x*=

_{i}*Mx*, where

_{o}*x*=

_{o}*a*(

*t*) is the location of the object at time

*t*. For a weakly absorbing object, the phase profile can be reasonably inferred from a single intensity image [29]. If the object absorption is non-negligible, one can derive the phase profile from a sequence of projection images acquired for different object-to-detector distance [30, 31]. Once the phase profile at the detector plane in Fig. 1 is derived, the iterative algorithm presented in this paper can be applied irrespective of whether the object is weakly or heavily absorbing. The same algorithm can be applied to the phase profiles acquired with other modalities, e.g. grating-based [13, 14], coded-aperture [15], or wire mesh-based [16] system. An exception has to be made for thick metal or thick cortical bone that completely stops the X-rays. In that case, no information about phase or absorption can be derived from the projection data in the region where there is no X-ray penetration.

In the X-ray regime, the complex refractive index of an object can be expressed in the following form [32]:

where the variables*δ*and

*β*are related to the phase delay and absorption of the X-ray beam, respectively, due to the object under interrogation. We have recently reported a forward image model for X-ray phase imaging based on the first Rytov approximation [28]. Using the exact Mie solution as comparison, we have shown that the first Rytov approximation is valid and accurate in the X-ray regime [33], in part because the refractive index of most materials is close to one. For the object at (

*a*,0), the light field recorded on the detector pixel at (

*x*,

*R*

_{2}) can be expressed in the following form [Eq. (10) in [28]]:

*A*represents the strength of the source,

_{s}*λ*is the wavelength of the X-ray beam,

*k*= 2

*π*/

*λ*is the wavenumber. The variable

*r*is the distance between the source (

*0*,-

*R*

_{1}) and the point (

*x*,

*R*

_{2}) on the detector; $r=\sqrt{{x}^{2}+{({R}_{1}+{R}_{2})}^{2}}$. The function ${\overline{u}}^{(s)}(x;a)$ is related to the scattering potential of the object $f(X,Z)={k}^{2}(1-{(n(X,Z)/{n}_{0})}^{2})$ as follows [28]:

*w*and

*W*are defined as $w=\sqrt{{\left(1/\lambda \right)}^{2}-{\left(U+x/\lambda r\right)}^{2}}$ and $W=w-{m}_{1}/\lambda $, respectively. The variable

*m*is defined as ${m}_{1}=\sqrt{1-{(x/r)}^{2}}$.

_{1}In a real experiment, the function ${\overline{u}}^{(s)}(x;a)$ can be obtained from two measurements, one ${u}_{samp}\left(x;a\right)$ with the sample in the field-of-view and the other ${u}_{bg}\left(x;a\right)$ without the sample: ${\overline{u}}^{(s)}\left(x;a\right)=\mathrm{log}\left[{u}_{samp}\left(x;a\right)/{u}_{bg}\left(x;a\right)\right]$, where log is the natural logarithm. The background images are acquired only once and used for all the subsequent measurements. In this study, we numerically generate the phase profiles using Eqs. (2) and (3). Consider a line source emitting a fan-shaped, mono-energetic beam at 30 keV and a detector with a pixel size of 5 *μ*m. *R*_{1} and *R*_{2} are fixed at 0.5 m. The object is a cylinder with a diameter of 1 mm and complex refractive index 1 - 0.56 × 10^{−7} + i × 0.00 (i.e., *δ* = 2.56 × 10^{−7} and *β* = 0.00). Figure 2(a) shows an example phase profile, the argument of the complex function ${\overline{u}}^{(s)}(x;a)$, for *a* = 0.289 m. Figure 2(b) shows a 2-D mapping of the phase profiles for varying locations of the cylinder. The width and height of each box in Fig. 2(b) represent object translation step and detector field-of-view, respectively. The dotted box in the figure corresponds to the phase profile in Fig. 2(a). For tomographic reconstruction of an object, multiple projection images are typically recorded for varying viewing angles on the object. In the geometry considered in this study [Fig. 1], the viewing angle is varied by translating the object across the fan-shaped X-ray beam. Given the phase profiles [Fig. 2(b)] and a forward image model [Eqs. (2) and (3)], one can retrieve the refractive index map in an iterative manner. We note, however, that direct evaluation of the integral in Eq. (3) is computationally expensive.

## 3. Iterative reconstruction algorithm for X-ray phase tomography of a moving object

We introduce the following 2-D mapping for the data acquired in the previous section:

Figure 3(a) demonstrates how the mapping represented by Eq. (4) works on ${\overline{u}}^{(s)}\left(x;a\right)$. The mapping operation: (i) selects a row in Fig. 2(b) using a box whose height is determined by the pixel size of the detector, (ii) translates each box horizontally so that the centroid of each box is aligned on the line*a*= 0, (iii) flips the data with respect to the line

*a*= 0, and changes the variables from

*x*to $\alpha $ and from

*a*to

*s*. Figure 3(b) shows ${U}_{s}\left(s;\alpha \right)$ after the mapping.

Applying the mapping to Eq. (3), we obtain the following expression for ${U}_{s}\left(s;\alpha \right)$:

*w*and

*W*are defined as $w=\sqrt{{\left(1/\lambda \right)}^{2}-{\left(U+\alpha /\lambda r\right)}^{2}}$, $W=w-{m}_{1}/\lambda $, respectively. The variables ${m}_{1}$, $r$, and ${\alpha}^{\prime}$ are defined in terms of $\alpha $ as ${m}_{1}=\sqrt{1-{(\alpha /r)}^{2}}$, $r=\sqrt{{\alpha}^{2}+{\left({R}_{1}+{R}_{2}\right)}^{2}}$and ${\alpha}^{\prime}=\alpha /M$, respectively. Noteworthy, the variables

*w*and

*W*in Eq. (5) do not depend on

*s*; thus, the integral in the equation can be efficiently evaluated for each

*α*using a fast Fourier transform algorithm. We note that Eq. (5) has a similar form to the complex scattered phase in plane-wave diffraction tomography (PW-DT) [34, 35]. From the analogy with PW-DT, Eq. (6) can be derived, which simply relates ${U}_{s}\left(s;\alpha \right)$ to $f(X,Z)$ in the 2-D spatial frequency plane:

*f*. We call this method Fourier mapping and direct inversion as in PW-DT [35].

The accuracy of tomographic reconstruction heavily depends on the angular coverage of the incident beam onto the object [36]. For the geometry considered in this study, the angular coverage is determined by the divergence angle of X-rays and the angle subtended by the detector trajectory at the source. The divergence angle of X-rays can be simply increased by rotating the source, but increasing the detector trajectory requires a larger system footprint. Therefore, the angular coverage cannot be increased too much in the geometry as shown in Fig. 1. In PW-DT, we have shown that an iterative reconstruction algorithm using the positivity and piecewise-smoothness constraints can effectively suppress the artifacts due to the limited angular coverage [37]. Here, we take a similar approach for X-ray phase tomography of a moving sample. The problem of retrieving the complex refractive index of an object from *N* scattered field measurements can be cast into the following form:

*μ*is the regularization parameter representing the trade-off between the data fidelity and penalty terms. $J\left(f\right)$ is the penalty functional incorporating

*a priori*information about the specimen. The objective function that minimizes Eq. (8) can be iteratively found using Eq. (9):

*τ*is the relaxation parameter determining the speed of convergence. The superscripts

*k*and

*k*+ 1 indicate the iteration number. In this study, the following constraints are used: (i) the refractive index profile is piecewise smooth; and (ii) the electron density and the attenuation coefficient are positive. The piecewise-smoothness constraint (i) can be imposed by using the penalty functional of the following form: $J\left(f\right)=\left(1/2\right){\displaystyle \int \sqrt{{\left|\nabla f\right|}^{2}+{\beta}^{2}}dV}$, where the parameter

*β*is an arbitrary small number preventing blow-up of $\nabla J$ in the case $\nabla f$ = 0. The positivity constraint (ii) is separately imposed to the real and imaginary parts of the refractive index map at every step of iteration. The iteration is terminated when the following convergence criterion is fulfilled:where

*ϵ*is a number that decides the iteration error.

## 4. Numerical tests of the proposed algorithm

In this section, we test the reconstruction algorithm proposed in the previous section using the data numerically generated with Eqs. (2) and (3). Figure 4(a) shows the object spectrum $\tilde{f}\left(U,W\right)$ obtained with applying Eq. (6) to the cylinder data [Fig. 3]. From the figure, it can be clearly seen that a significant portion of the object spectrum is missing, and the refractive index map is elongated along the *z* direction in Fig. 4(b). The refractive index histogram [Fig. 4(c)] indicates that the refractive index is seriously underestimated; the true *δ* value is 2.56 × 10^{−7}, while the peak in Fig. 4(c) is at around 1.05 × 10^{−7}. Figure 4(d) shows $\tilde{f}\left(U,W\right)$ after 500 iterations of the algorithm [Eq. (9)]. In the figure, the region that was previously empty is filled with new data, and the refractive index map in Fig. 4(e) clearly shows the round shape of the cylinder. The refractive index histogram in Fig. 4(f) indicates a peak at around 2.56 × 10^{−7} (true value).

Next, we tested our algorithm using a more realistic phantom. The phantom is a water cylinder containing four cylinders of different materials: polypropylene, Mylar, Teflon, and PMMA (polymethyl methacrylate). The diameter of the water cylinder is 200 mm, and that of the four inner cylinders is 40 mm. The refractive index value of each material is listed in Table 1 [32]. Figures 5(a) and 5(b) show the refractive index map of original phantom and that after 1000 iterations, respectively. Compared to the result without regularization [Fig. 5(c)], Fig. 5(b) clearly shows the boundaries of outer and inner cylinders. Figure 5(d) shows the refractive index value (*δ*) averaged over each inner cylinder region at different iteration steps. Here, the refractive value was normalized by the true *δ* value for each material. The four materials have different speeds of convergence, the origin of which is not clear yet. Nonetheless, after 1000 iterations, the *δ* values of the four materials match with the true *δ* values within 10% accuracy. Considering that the data were collected only for the small detector movement that corresponds to the angular coverage of (−30° to 30°), this accuracy is remarkable. The accuracy may be further improved by utilizing more information about the specimen [39].

## 5. Conclusion

In this study, we derived an iterative algorithm for X-ray phase tomography of an object translated across a fan-shaped X-ray beam. The developed algorithm will be useful in a high-throughput X-ray phase tomography system continually interrogating objects on a conveyor belt or using multiple off-axis sources.

## Acknowledgments

This work was supported by the Department of Homeland Security’s Science and Technology Directorate through contract HSHQDC-11-C-0083, the National Research Foundation of Singapore through the Singapore-MIT Alliance for Research and Technology Centre, and the DARPA AXiS program (Grant No. N66001-11-4204, P.R. No. 1300217190). The authors thank Prof. George Barbastathis and Dr. Niyom Lue for useful discussion.

## References and links

**1. **A. Stanton, “Wilhelm Conrad Röntgen on a new kind of rays: translation of a paper read before the Würzburg Physical and Medical Society, 1895,” Nature **53**, 274–276 (1896).

**2. **C. A. Helms, *Fundamentals of Skeletal Radiology* (Saunders, 2005).

**3. **E. D. Pisano, M. J. Yaffe, and C. M. Kuzmiak, *Digital Mammography* (Lippincott Williams & Wilkins, 2004).

**4. **H. Vogel and D. Haller, “Luggage and shipped goods,” Eur. J. Radiol. **63**(2), 242–253 (2007). [CrossRef] [PubMed]

**5. **J. H. Hubbell, Wm. J. Veigele, E. A. Briggs, R. T. Brown, D. T. Cromer, and R. J. Howerton, “Atomic form factors, incoherent scattering functions, and photon scattering cross sections,” J. Phys. Chem. Ref. Data **4**(3), 471–538 (1975). [CrossRef]

**6. **D. M. Paganin, *Coherent X-ray Optics* (Oxford University, 2006).

**7. **A. Momose, T. Takeda, and Y. Itai, “Blood vessels: depiction at phase-contrast X-ray imaging without contrast agents in the mouse and rat-feasibility study,” Radiology **217**(2), 593–596 (2000). [PubMed]

**8. **S. W. Wilkins, T. E. Gureyev, D. Gao, A. Pogany, and A. Stevenson, “Phase-contrast imaging using polychromatic hard X-rays,” Nature **384**(6607), 335–338 (1996). [CrossRef]

**9. **U. Bonse and F. Beckmann, “Multiple-beam X-ray interferometry for phase-contrast microtomography,” J. Synchrotron Radiat. **8**(1), 1–5 (2001). [CrossRef] [PubMed]

**10. **M. Ando, A. Maksimenko, H. Sugiyama, W. Pattanasiriwisawa, K. Hyodo, and C. Uyama, “Simple X-ray dark-and bright-field imaging using achromatic Laue optics,” Jpn. J. Appl. Phys. **41**(Part 2, No. 9A/B), L1016–L1018 (2002). [CrossRef]

**11. **Y. S. Kashyap, P. S. Yadav, T. Roy, P. S. Sarkar, M. Shukla, and A. Sinha, “Laboratory-based X-ray phase-contrast imaging technique for material and medical science applications,” Appl. Radiat. Isot. **66**(8), 1083–1090 (2008). [CrossRef] [PubMed]

**12. **Z. Zaprazny, D. Korytar, V. Ac, P. Konopka, and J. Bielecki, “Phase contrast imaging of lightweight objects using microfocus X-ray source and high resolution CCD camera,” J. Instrum. **7**(03), C03005 (2012). [CrossRef]

**13. **T. Weitkamp, C. David, O. Bunk, J. Bruder, P. Cloetens, and F. Pfeiffer, “X-ray phase radiography and tomography of soft tissue using grating interferometry,” Eur. J. Radiol. **68**(3Suppl), S13–S17 (2008). [CrossRef] [PubMed]

**14. **F. Pfeiffer, M. Bech, O. Bunk, P. Kraft, E. F. Eikenberry, Ch. Brönnimann, C. Grünzweig, and C. David, “Hard-X-ray dark-field imaging using a grating interferometer,” Nat. Mater. **7**(2), 134–137 (2008). [CrossRef] [PubMed]

**15. **A. Olivo and R. Speller, “A coded-aperture technique allowing x-ray phase contrast imaging with laboratory sources,” Appl. Phys. Lett. **91**(7), 074106 (2007). [CrossRef]

**16. **H. Wen, E. E. Bennett, M. M. Hegedus, and S. C. Carroll, “Spatial harmonic imaging of X-ray scattering--initial results,” IEEE Trans. Med. Imaging **27**(8), 997–1002 (2008). [CrossRef] [PubMed]

**17. **A. Bravin, P. Coan, and P. Suortti, “X-ray phase-contrast imaging: from pre-clinical applications towards clinics,” Phys. Med. Biol. **58**(1), R1–R35 (2013). [CrossRef] [PubMed]

**18. **A. C. Kak and M. Slaney, *Principles of Computerized Tomographic Imaging* (Society for Industrial and Applied Mathematics, 1988).

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

**20. **F. Pfeiffer, C. Kottler, O. Bunk, and C. David, “Hard X-ray phase tomography with low-brilliance sources,” Phys. Rev. Lett. **98**(10), 108105 (2007). [CrossRef] [PubMed]

**21. **S. Mayo, T. Davis, T. Gureyev, P. Miller, D. Paganin, A. Pogany, A. Stevenson, and S. Wilkins, “X-ray phase-contrast microscopy and microtomography,” Opt. Express **11**(19), 2289–2302 (2003). [CrossRef] [PubMed]

**22. **G. Donges and R. Dietrich, “Baggage inspection system,” U. S. patent 4,759,047 (1988).

**23. **R. Dietrich, “Baggage inspection system,” U. S. patent 4,783,794 (1988).

**24. **J. M. Park, E. A. Franken Jr, M. Garg, L. L. Fajardo, and L. T. Niklason, “Breast tomosynthesis: present considerations and future applications,” Radiographics **27**(Suppl 1), S231–S240 (2007). [CrossRef] [PubMed]

**25. **D. Nahamoo, S. Pan, and A. C. Kak, “Synthetic aperture diffraction tomography and its interpolation-free computer implementation,” IEEE Trans. Sonics Ultrason. **31**(4), 218–229 (1984). [CrossRef]

**26. **T. Yasuda, T. Yasui, T. Araki, and E. Abraham, “Real-time two-dimensional terahertz tomography of moving objects,” Opt. Commun. **267**(1), 128–136 (2006). [CrossRef]

**27. **N. Lue, W. Choi, G. Popescu, K. Badizadegan, R. R. Dasari, and M. S. Feld, “Synthetic aperture tomographic phase microscopy for 3D imaging of live cells in translational motion,” Opt. Express **16**(20), 16240–16246 (2008). [CrossRef] [PubMed]

**28. **Y. Sung, C. J. R. Sheppard, G. Barbastathis, M. Ando, and R. Gupta, “Full-wave approach for X-ray phase imaging,” Opt. Express **21**(15), 17547–17557 (2013). [CrossRef] [PubMed]

**29. **A. Burvall, U. Lundström, P. A. C. Takman, D. H. Larsson, and H. M. Hertz, “Phase retrieval in X-ray phase-contrast imaging suitable for tomography,” Opt. Express **19**(11), 10359–10376 (2011). [CrossRef] [PubMed]

**30. **M. R. Teague, “Deterministic phase retrieval: a Green’s function solution,” J. Opt. Soc. Am. **73**(11), 1434–1441 (1983). [CrossRef]

**31. **N. Streibl, “Phase imaging by the transport equation of intensity,” Opt. Commun. **49**(1), 6–10 (1984). [CrossRef]

**32. **B. Henke, E. Gullikson, and J. C. Davis, “X-ray Interactions: Photoabsorption, Scattering, Transmission, and Reflection at E = 50-30,000 eV, Z = 1-92,” At. Data Nucl. Data Tables **54**(2), 181–342 (1993). [CrossRef]

**33. **Y. Sung and G. Barbastathis, “Rytov approximation for x-ray phase imaging,” Opt. Express **21**(3), 2674–2682 (2013). [CrossRef] [PubMed]

**34. **E. Wolf, “Three-dimensional structure determination of semi-transparent objects from holographic data,” Opt. Commun. **1**(4), 153–156 (1969). [CrossRef]

**35. **Y. Sung, W. Choi, C. Fang-Yen, K. Badizadegan, R. R. Dasari, and M. S. Feld, “Optical diffraction tomography for high resolution live cell imaging,” Opt. Express **17**(1), 266–277 (2009). [CrossRef] [PubMed]

**36. **P. Charbonnier, L. Blanc-Feraud, G. Aubert, and M. Barlaud, “Deterministic edge-preserving regularization in computed imaging,” IEEE Trans. Image Process. **6**(2), 298–311 (1997). [CrossRef] [PubMed]

**37. **Y. Sung and R. R. Dasari, “Deterministic regularization of three-dimensional optical diffraction tomography,” J. Opt. Soc. Am. A **28**(8), 1554–1561 (2011). [CrossRef] [PubMed]

**38. **M. Bertero and P. Boccacci, *Introduction to Inverse Problems in Imaging* (Taylor & Francis, 1998).

**39. **L. Ritschl, F. Bergner, and M. Kachelrieß, “A new approach to limited angle tomography using the compressed sensing framework,” Proc. SPIE **7622**, 76222H, 76222H-9 (2010). [CrossRef]