## Abstract

We demonstrate phase contrast enhancement of X-ray computed tomography derived from propagation based imaging. In this method, the absorption and phase components are assumed to be correlated, allowing for phase retrieval from a single image. Experimental results are shown for liquid samples. Signal-to-noise ratio is greatly enhanced relative to pure attenuation based imaging.

© 2014 Optical Society of America

## 1. Introduction

Tomographic imaging using hard X-rays is a highly active field, as there are few other methods to noninvasively examine thick objects that are opaque to other forms of radiation (visible light, sound, etc.). Conventional X-ray computed tomography recovers the attenuation of an object, which typically provides poor contrast between materials with similar atomic number. In applications such as security inspection, there are workarounds; for example, nitrogen is a major component of many explosives, and has a unique signature which can be detected using X-ray scatter imaging [1]. For light materials and non-nitrogenous explosives, there is no easy workaround—hence the inconvenience of having to throw out water containers before security screening at airports. For medical X-ray imaging, contrast agents are injected to enhance contrast in tissues of interest at the price of discomfort to the patient, dangers from the procedure, and nephrotoxicity of the contrast agent used [2, 3].

Phase information can provide significantly enhanced contrast for light materials without the use of additional preparations [4–6]. While measurement of phase is not as straightforward as attenuation-based imaging, propagation-based techniques provide a relatively simple way to reconstruct the phase of an object. By using images taken at multiple propgation distances, there is no need for optical elements or highly conditioned X-ray sources [7,8]. Even the spatial coherence of the x–ray source becomes less of an issue, since in actuality it is the optical path length (OPL) of the specimen that we seek to recover, not the phase *per se* [9].

In this paper, we adopt the transport of intensity equation (TIE), which relates the measured intensity to the phase [10]:

*I*(

*x*,

*y*,

*z*) is the intensity image at the propagation distance

*z*,

*ϕ*(

*x*,

*y*) is the phase map (equivalently, OPL) of the object, and

*k*is the wavenumber. By solving the TIE at multiple angles of exposure, the refractive index distribution of an object can be obtained through tomographic reconstruction. In this way, TIE tomographic reconstruction can be thought of as a two-step problem.

Direct inversion of the TIE tomographic reconstruction problem is ill-posed for two reasons:

- (1) due to the existence of a zero in the TIE transfer function (the projection component); and
- (2) undersampling of the object index of refraction distribution in radon space in cases where only a small number of projections are acquired (the tomographic component).

For TIE, the traditional choice of regularizer to relieve the ill–posedness is Tikhonov [11–13]. However, this type of regularization also strongly deteriorates the low-frequency signal in the phase image. Myers *et al.* propose inversion of the TIE using prior knowledge that the sample consists of a single material of known refractive index [14]. For tomography, the traditional regularized solution is the filtered backprojection (FBP), which utilizes the ramp filtered dual of the Radon transform. The drawback to this technique is that the quality of FBP reconstructions depends on having a large number of projection angles. Recently, iterative solvers have been shown to relax these requirements, specifically the solvers which exploit the sparsity of the expected tomographic reconstruction [15, 16].

We have previously demonstrated the effectiveness of a single-step TIE based method for compressive X-ray phase tomography of weakly attenuating samples, wherein the TIE and tomographic reconstruction steps are computed in a combined transfer function, and solved using iterative methods [17]. However, this method places a severe restriction on the class of objects suitable for imaging.

In this paper, we extend this method to account for samples which consist of only light materials, based on the phase-attenuation duality for light materials [18]. The formulation of the TIE solver exploiting the duality is shown in Section 2. We also demonstrate experimentally in Section 3 the use of duality to resolve the refractive index of water from that of chemically similar (in terms of Z number) but potentially hazardous materials, such as hydrogen peroxide and acetone; this particular application is potentially significant for security in airports and other controlled spaces.

## 2. Reconstruction process

A schematic diagram of the imaging geometry is presented in Fig. 1. A microfocus source with a spectrally weighted mean wavelength *λ* is located at the plane *z* = −*z*_{0}, and a detector is located at the plane *z* = *d*. The sample, centered at *z* = 0, is characterized by a complex refractive index *n*(*x*, *y*, *z*; *λ*) = 1 − *δ*(*x*, *y*, *z*; *λ*) + *iβ*(*x*, *y*, *z*; *λ*), where 1 − *δ*(*x*, *y*, *z*; *λ*) and *β*(*x*, *y*, *z*; *λ*) are the real and imaginary parts of the refractive index, respectively. We assume that the geometry of our experiment is such that the beam passing through the object is approximated by a plane wave oriented along the optical axis, and that the interaction between the sample and the field follows the projection approximation.

While TIE is typically formulated using two intensity measurments at different positions along the optical axis [19], the phase-attenuation duality (PAD) approximation allows us to reformulate the TIE such that only one intensity measurement is necessary. We assume that electron interactions are the primary source of x-ray attenuation and phase delay, thus the real and complex components of the complex index are correlated such that *δ*(*x*, *y*, *z*; *λ*)/*β*(*x*, *y*, *z*; *λ*) = *γ*(*λ*) [18]. Under paraxial and small-wavelength approximations, the TIE gives the following relationship between the projected phase map of the object, *ϕ*(*x*, *y*), and the image of the object at the detector, *I*(*x′*, *y′*),

*x′*=

*Mx*and

*y′*=

*My*correspond to spatial coordinates in the detector plane,

*M*= (

*z*

_{0}+

*d*)/

*z*

_{0}is the magnification at the detector plane,

*I*

_{0}is the incident intensity without the object, and

*d′*=

*d/M*is the effective propagation distance given by the Fresnel scaling theorem [20] (page 397). Equation (2) is the discretized form of PAD TIE which uses a single image to retrieve phase. In a tomographic measurement, the sample is rotated about the

*x*axis such that the projected phase

*ϕ*at a certain angle

*θ*is

*D*(·) to denote the Dirac delta function to avoid confusion with the real component of the refractive index. Equation (2) along with the forward Radon transform, Eq. (3), provide us with a forward model for describing the image generated on our X-ray detector by our object, and implementation of these operators in the Fourier domain is computationally efficient [21]. To perform direct reconstruction, Eq. (2) can be inverted using the Fourier domain solution to Poisson’s equation, and Eq. (3) can be inverted using the Fourier slice theorem implementation of FBP as

*u*,

*v*, and

*w*are spatial frequency variables, and

*ℱ*

_{1D}and

*ℱ*

_{2D}both denote Fourier transforms but in two different domains: the former is for applying the Fourier slice theorem in tomography and, hence, it operates along the projection coordinate variable; whereas the latter applies to the lateral plane (

*x*,

*y*). The superscript

*ℱ*

^{−1}denotes the inverse Fourier transform correspondingly in the two cases.

To perform iterative reconstruction, we utilize a modified form of our approach in [17], which will be described briefly for convenience. The detector is assumed to consist of an *N* × *N* grid of square pixels of side length *M*Δ. Let Θ denote the number of angular projections. Then the measured projections *g* at all angles may be arranged into a real-valued vector **g** of length *N*^{2}Θ. The refractive index of the object will also be discretized into a *N* × *N* × *N* cube width side length Δ, packed into a complex valued vector **n**. Then let **P** and **R** denote operators corresponding to the discretized forms of Eqs. (2) and (3), respectively, such that

**A**is the cascade of the linear operators

**P**and

**R**. If

**n**is sparse in some basis, as it often is the case in tomography, then Eq. (6) can be solved using compressive reconstruction [15, 16, 22]. Specifically, total variation (TV) minimization has been shown to be an effective sparsity basis for tomographic reconstruction:

*, ∇*

_{x}*, and ∇*

_{y}*are the finite difference operators in Cartesian coordinates [23]. We use the two-step iterative shrinkage/thresholding algorithm (TwIST) to solve the minimization [22].*

_{z}## 3. Experimental application to discrimination of liquids

To investigate tomography and identification of liquids with similar density profiles but differing chemical compositions, we imaged four Eppendorf tubes containing water, hydrogen peroxide, acetone, and air. Water and hydrogen peroxide are chemically similar to one another, and therefore their absorption signatures are also essentially indistinguishable. Acetone and hydrogen peroxide are both common components used in explosive synthesis. The air filled tube was added as a control.

A total of 72 projections were taken, with source-to-detector distance 2.159 m and source-to-object distance 0.635 m. The X-ray microfocus source (Hammamatsu L812103) was set to 7 μm spot size, 100 kVp, and 100 μA. The spectrally weighted mean energy of 46 keV was used to determine *k* in phase reconstruction.

These experiments showed that X-ray phase imaging allowed for far higher sensitivity than standard absorption contrast (Figs. 2 and 3). In Fig. 3, we see that absorption based tomography is heavily noise-corrupted, and does not provide sufficient visual information to distinguish between water and hydrogen peroxide, while compressive phase reconstruction can clearly distinguish between the two samples. To quantify this difference, we measured a peak signal-to-noise ratio (PSNR) of 3 dB in a conventional CT, while obtaining a PSNR of 38 dB in our phase CT, with all other experimental parameters constant, a SNR gain of 35 dB.

In a Student’s T-test comparison between phase and attenuation based imaging, both methods were able to achieve discrimination with the full CT data, though phase was able to provide much higher certainty in the signal region (*P* < 0.010 for absorption, *P* < 6.79 × 10^{−7} for phase). Taking only 100 randomly sampled voxels from each CT (after segmentation), intensity based CT could only distinguish water and peroxide in 26% of tests, meaning rejection of the null hypothesis at the 5% significance level, while phase based CT could distinguish the two liquids in 100% of tests.

## 4. Discussion

We have demonstrated a technique for quantitative phase contrast tomographic reconstruction, and our results show a significant gain in contrast enhancement compared to its conventional counterpart. While typical tomography can be sensitive to detector noise, our method is robust to noise in two ways: the edge-enhancemened nature of the propagated X-ray image emphasizes material boundaries, and our application of compressive reconstruction with TV minimization favors sparse edges in the reconstructed image. We have applied our technique to detect small differences in phase content between two radiographically similar objects, showing a clear distinction between water and hydrogen peroxide. Furthermore, our approach to compressive tomographic reconstruction allows such images to be taken using only a fraction of the number of images needed for a conventional CT scan. While absorption imaging was still able to distinguish hydrogen peroxide from water given a full CT reconstruction, phase imaging is able to produce the same confidence in measurement in a much smaller sample, such as in the case of an incomplete CT reconstruction, and has significant potential applications in the radiographic detection of hazardous materials.

## Acknowledgments

We are grateful to Justin Lee, Lei Tian, and David J. Brady for helpful discussions. This work was supported by the Department of Homeland Security, Science and Technology Directorate through contract HSHQDC-11-C-00083.

## References and links

**1. **S. Singh and M. Singh, “Explosives detection systems (EDS) for aviation security,” Signal Process. **83**, 3155 (2003). [CrossRef]

**2. **K. M. Hasebroock and N. J. Serkova, “Toxicity of MRI and CT contrast agents,” Expert Opin. Drug Metab. Toxicol. **5**, 403–416 (2009). [CrossRef]

**3. **E. M. Lautin, N. J. Freeman, A. H. Schoenfeld, C. W. Bakal, N. Haramati, A. C. Friedman, J. L. Lautin, S. Braha, E. G. Kadish, and S. Sprayregen, “Radiocontrast-associated renal dysfunction: incidence and risk factors,” American Journal of Roentgenology **157**, 49–58 (1991). [CrossRef]

**4. **T. 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**, 595–598 (1995). [CrossRef]

**5. **A. Momose, T. Takeda, Y. Itai, and K. Hirano, “Phasecontrast xray computed tomography for observing biological soft tissues,” Nat Med **2**, 473–475 (1996). [CrossRef]

**6. **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**, 258–261 (2006). [CrossRef]

**7. **J. P. Guigay, M. Langer, R. Boistel, and P. Cloetens, “Mixed transfer function and transport of intensity approach for phase retrieval in the fresnel region,” Opt. Lett. **32**, 1617–1619 (2007). [CrossRef]

**8. **R. C. Chen, L. Rigon, and R. Longo, “Comparison of single distance phase retrieval algorithms by considering different object composition and the effect of statistical and structural noise,” Opt. Express **21**, 7384–7399 (2013). [CrossRef]

**9. **J. C. Petruccelli, L. Tian, and G. Barbastathis, “The transport of intensity equation for optical path length recovery using partially coherent illumination,” Opt. Express **21**, 14430 (2013). [CrossRef]

**10. **M. Reed Teague, “Deterministic phase retrieval: a green?s function solution,” J. Opt. Soc. Am. **73**, 1434–1441 (1983). [CrossRef]

**11. **A. V. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” J. Opt. Soc. Am. A **19**, 472–480 (2002). [CrossRef]

**12. **A. Burvall, U. Lundstrm, 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**, 10359–10376 (2011). [CrossRef]

**13. **T. E. Gureyev, T. J. Davis, A. Pogany, S. C. Mayo, and S. W. Wilkins, “Optical phase retrieval by use of first born- and rytov-type approximations,” Appl. Opt. **43**, 2418–2430 (2004). [CrossRef]

**14. **G. R. Myers, D. M. Paganin, T. E. Gureyev, and S. C. Mayo, “Phase-contrast tomography of single-material objects from few projections,” Opt. Express **16**, 908–919 (2008). [CrossRef]

**15. **E. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information Theory **52**, 489–509 (2006). [CrossRef]

**16. **T. Goldstein and S. Osher, “The split bregman method for l1-regularized problems,” SIAM J. Imaging Sci. **2**, 323–343 (2009). [CrossRef]

**17. **L. Tian, J. C. Petruccelli, Q. Miao, H. Kudrolli, V. Nagarkar, and G. Barbastathis, “Compressive x-ray phase tomography based on the transport of intensity equation,” Opt. Lett. **38**, 3418–3421 (2013). [CrossRef]

**18. **X. Wu, H. Liu, and A. Yan, “X-ray phase-attenuation duality and phase retrieval,” Opt. Lett. **30**, 379–381 (2005). [CrossRef]

**19. **T. E. Gureyev and S. W. Wilkins, “On x-ray phase imaging with a point source,” J. Opt. Soc. Am. A **15**, 579–585 (1998). [CrossRef]

**20. **D. Paganin, *Coherent X-Ray Optics* (Oxford University Press, 2006). [CrossRef]

**21. **S. Matej, J. Fessler, and I. Kazantsev, “Iterative tomographic image reconstruction using fourier-based forward and back-projectors,” IEEE Transactions on Medical Imaging **23**, 401–412 (2004). [CrossRef]

**22. **J. Bioucas-Dias and M. A. T. Figueiredo, “A new TwIST: Two-step iterative shrinkage/thresholding algorithms for image restoration,” IEEE Transactions on Image Processing **16**, 2992–3004 (2007). [CrossRef]

**23. **L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D: Nonlinear Phenomena **60**, 259–268 (1992). [CrossRef]