Contrast enhancement in X-ray phase contrast tomography ∗
Adam Pan,1, Ling Xu,1 Jon C. Petruccelli,2 Rajiv Gupta,3 Bipin Singh,4 and George Barbastathis5,6,7 1 Department
of Health Sciences and Technology, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA 2 Department of Physics, University at Albany – State University of New York, 1400 Washington Avenue, Albany, NY 12222, USA 3 Department of Imaging, Massachusetts General Hospital, 55 Fruit Street, Boston, MA 02215, USA 4 Radiation Monitoring Devices, Inc., 44 Hunt Street, Watertown, MA 02472, USA 5 Department of Mechanical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA 6 Singapore–MIT Alliance for Research and Technology (SMART) Centre, 1 Create Way, Singapore 138602, Singapore 7 Shanghai Jiao Tong University – University of Michigan Joint Institute 800 Dong Chuan Road, Minhang District, Shanghai 200240, China ∗ [email protected]
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 OCIS codes: (100.5070) Phase retrieval; (340.7440) X-ray imaging; (100.3010) Image reconstruction techniques.
References and links 1. S. Singh and M. Singh, “Explosives detection systems (EDS) for aviation security,” Signal Process. 83, 3155 (2003). 2. K. M. Hasebroock and N. J. Serkova, “Toxicity of MRI and CT contrast agents,” Expert Opin. Drug Metab. Toxicol. 5, 403–416 (2009). 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). 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). 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). 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). 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). 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). 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).
#212233 - $15.00 USD (C) 2014 OSA
Received 20 May 2014; revised 2 Jul 2014; accepted 9 Jul 2014; published 17 Jul 2014 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018020 | OPTICS EXPRESS 18020
10. M. Reed Teague, “Deterministic phase retrieval: a green?s function solution,” J. Opt. Soc. Am. 73, 1434–1441 (1983). 11. A. V. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” J. Opt. Soc. Am. A 19, 472–480 (2002). 12. A. Burvall, U. Lundstrm, P. A. C. Takman, D. H. Larsson, and H. M. Hertz, “Phase retrieval in x-ray phasecontrast imaging suitable for tomography,” Opt. Express 19, 10359–10376 (2011). 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). 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). 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). 16. T. Goldstein and S. Osher, “The split bregman method for l1-regularized problems,” SIAM J. Imaging Sci. 2, 323–343 (2009). 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). 18. X. Wu, H. Liu, and A. Yan, “X-ray phase-attenuation duality and phase retrieval,” Opt. Lett. 30, 379–381 (2005). 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). 20. D. Paganin, Coherent X-Ray Optics (Oxford University Press, 2006). 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). 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). 23. L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D: Nonlinear Phenomena 60, 259–268 (1992).
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]: h i ∂ I(x, y, z) ∂ ∂ k = ∇⊥ · I(x, y, z)∇⊥ φ (x, y) , ∇⊥ = , , (1) ∂z ∂x ∂y where 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
#212233 - $15.00 USD (C) 2014 OSA
Received 20 May 2014; revised 2 Jul 2014; accepted 9 Jul 2014; published 17 Jul 2014 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018020 | OPTICS EXPRESS 18021
Sample (n) x (rotational axis)
X-ray detector x´
X-ray microfocus source
zs θ z y´
y ys z=-z0
z=0
z=d
Fig. 1. Imaging geometry for TIE tomography.
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 = −z0 , 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 #212233 - $15.00 USD (C) 2014 OSA
Received 20 May 2014; revised 2 Jul 2014; accepted 9 Jul 2014; published 17 Jul 2014 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018020 | OPTICS EXPRESS 18022
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(x0 , y0 ), d 0 λ γ(λ ) 2 2φ (x, y) I(Mx, My) 0 0 = 1+ ∇⊥ exp , (2) g(x , y ) = I0 4π γ(λ ) where x0 = Mx and y0 = My correspond to spatial coordinates in the detector plane, M = (z0 + d)/z0 is the magnification at the detector plane, I0 is the incident intensity without the object, and d 0 = 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 ZZ
φ (x, y; θ ) =
δ (x, ys , zs ; λ )D(y − ys cos θ + zs sin θ )dys dzs ,
(3)
where we use 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 n o F2D I(x,y) I γ(λ ) 0 −1 (4) ln F2D φ (x, y) = 1 + πzλ γ(λ ) |u2 + v2 | 2 −1 {|w| F1D {φ (x, y, θ )}} , δ (x, ys , zs ; λ ) = F1D
(5)
where u, v, and w are spatial frequency variables, and F1D and F2D 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 F −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 g = PRn ≡ An
#212233 - $15.00 USD (C) 2014 OSA
(6)
Received 20 May 2014; revised 2 Jul 2014; accepted 9 Jul 2014; published 17 Jul 2014 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018020 | OPTICS EXPRESS 18023
where the operator 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: nˆ = arg min ||n||TV such that g = An
(7)
n
p where the TV norm is defined as ||n||TV = ∑ (∇x n)2 ) + (∇y n)2 ) + (∇z n)2 ), and ∇x , ∇y , and n
∇z are the finite difference operators in Cartesian coordinates [23]. We use the two-step iterative shrinkage/thresholding algorithm (TwIST) to solve the minimization [22]. 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 sourceto-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 signalto-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 1.2 1.1
0.8
1
0.6
0.4
Normalized Intensity
1
0.9 0.8 0.7
0.2 0
0.6
5 mm
0.5
(a)
(b)
Fig. 2. A projection image of the sample (a), with its profile (b) showing significant white noise at the detector level. The liquids occupy the upper portions of the tubes, and air bubbles occupy the lower portions of the tubes.
#212233 - $15.00 USD (C) 2014 OSA
Received 20 May 2014; revised 2 Jul 2014; accepted 9 Jul 2014; published 17 Jul 2014 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018020 | OPTICS EXPRESS 18024
−10
−10
x 10
x 10 6
6
2 0 −2 −4
Absorption coefficient (β)
4 4
2 0 −2 −4
−6 0
−6
5 mm −8
−8
(a)
(c) −7
−7
x 10
x 10 0.5
0
−1 −1.5 −2
0
5 mm
Real refractive index (δ)
0 −0.5
−0.5 −1 −1.5 −2 −2.5
−2.5
−3
(b)
(d)
Fig. 3. Reconstruction results. (a) Filtered backprojection, (b) compressive PAD TIE, (c) plot along line in (a), and (d) plot along line in (b). These reconstructions show cross sections of the Eppendorf tubes containing the liquids taken along the plane coinciding with the line in Fig. 2(a)
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 #212233 - $15.00 USD (C) 2014 OSA
Received 20 May 2014; revised 2 Jul 2014; accepted 9 Jul 2014; published 17 Jul 2014 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018020 | OPTICS EXPRESS 18025
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.
#212233 - $15.00 USD (C) 2014 OSA
Received 20 May 2014; revised 2 Jul 2014; accepted 9 Jul 2014; published 17 Jul 2014 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018020 | OPTICS EXPRESS 18026
Adam Pan,1, Ling Xu,1 Jon C. Petruccelli,2 Rajiv Gupta,3 Bipin Singh,4 and George Barbastathis5,6,7 1 Department
of Health Sciences and Technology, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA 2 Department of Physics, University at Albany – State University of New York, 1400 Washington Avenue, Albany, NY 12222, USA 3 Department of Imaging, Massachusetts General Hospital, 55 Fruit Street, Boston, MA 02215, USA 4 Radiation Monitoring Devices, Inc., 44 Hunt Street, Watertown, MA 02472, USA 5 Department of Mechanical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA 6 Singapore–MIT Alliance for Research and Technology (SMART) Centre, 1 Create Way, Singapore 138602, Singapore 7 Shanghai Jiao Tong University – University of Michigan Joint Institute 800 Dong Chuan Road, Minhang District, Shanghai 200240, China ∗ [email protected]
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 OCIS codes: (100.5070) Phase retrieval; (340.7440) X-ray imaging; (100.3010) Image reconstruction techniques.
References and links 1. S. Singh and M. Singh, “Explosives detection systems (EDS) for aviation security,” Signal Process. 83, 3155 (2003). 2. K. M. Hasebroock and N. J. Serkova, “Toxicity of MRI and CT contrast agents,” Expert Opin. Drug Metab. Toxicol. 5, 403–416 (2009). 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). 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). 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). 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). 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). 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). 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).
#212233 - $15.00 USD (C) 2014 OSA
Received 20 May 2014; revised 2 Jul 2014; accepted 9 Jul 2014; published 17 Jul 2014 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018020 | OPTICS EXPRESS 18020
10. M. Reed Teague, “Deterministic phase retrieval: a green?s function solution,” J. Opt. Soc. Am. 73, 1434–1441 (1983). 11. A. V. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” J. Opt. Soc. Am. A 19, 472–480 (2002). 12. A. Burvall, U. Lundstrm, P. A. C. Takman, D. H. Larsson, and H. M. Hertz, “Phase retrieval in x-ray phasecontrast imaging suitable for tomography,” Opt. Express 19, 10359–10376 (2011). 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). 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). 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). 16. T. Goldstein and S. Osher, “The split bregman method for l1-regularized problems,” SIAM J. Imaging Sci. 2, 323–343 (2009). 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). 18. X. Wu, H. Liu, and A. Yan, “X-ray phase-attenuation duality and phase retrieval,” Opt. Lett. 30, 379–381 (2005). 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). 20. D. Paganin, Coherent X-Ray Optics (Oxford University Press, 2006). 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). 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). 23. L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica D: Nonlinear Phenomena 60, 259–268 (1992).
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]: h i ∂ I(x, y, z) ∂ ∂ k = ∇⊥ · I(x, y, z)∇⊥ φ (x, y) , ∇⊥ = , , (1) ∂z ∂x ∂y where 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
#212233 - $15.00 USD (C) 2014 OSA
Received 20 May 2014; revised 2 Jul 2014; accepted 9 Jul 2014; published 17 Jul 2014 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018020 | OPTICS EXPRESS 18021
Sample (n) x (rotational axis)
X-ray detector x´
X-ray microfocus source
zs θ z y´
y ys z=-z0
z=0
z=d
Fig. 1. Imaging geometry for TIE tomography.
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 = −z0 , 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 #212233 - $15.00 USD (C) 2014 OSA
Received 20 May 2014; revised 2 Jul 2014; accepted 9 Jul 2014; published 17 Jul 2014 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018020 | OPTICS EXPRESS 18022
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(x0 , y0 ), d 0 λ γ(λ ) 2 2φ (x, y) I(Mx, My) 0 0 = 1+ ∇⊥ exp , (2) g(x , y ) = I0 4π γ(λ ) where x0 = Mx and y0 = My correspond to spatial coordinates in the detector plane, M = (z0 + d)/z0 is the magnification at the detector plane, I0 is the incident intensity without the object, and d 0 = 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 ZZ
φ (x, y; θ ) =
δ (x, ys , zs ; λ )D(y − ys cos θ + zs sin θ )dys dzs ,
(3)
where we use 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 n o F2D I(x,y) I γ(λ ) 0 −1 (4) ln F2D φ (x, y) = 1 + πzλ γ(λ ) |u2 + v2 | 2 −1 {|w| F1D {φ (x, y, θ )}} , δ (x, ys , zs ; λ ) = F1D
(5)
where u, v, and w are spatial frequency variables, and F1D and F2D 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 F −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 g = PRn ≡ An
#212233 - $15.00 USD (C) 2014 OSA
(6)
Received 20 May 2014; revised 2 Jul 2014; accepted 9 Jul 2014; published 17 Jul 2014 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018020 | OPTICS EXPRESS 18023
where the operator 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: nˆ = arg min ||n||TV such that g = An
(7)
n
p where the TV norm is defined as ||n||TV = ∑ (∇x n)2 ) + (∇y n)2 ) + (∇z n)2 ), and ∇x , ∇y , and n
∇z are the finite difference operators in Cartesian coordinates [23]. We use the two-step iterative shrinkage/thresholding algorithm (TwIST) to solve the minimization [22]. 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 sourceto-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 signalto-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 1.2 1.1
0.8
1
0.6
0.4
Normalized Intensity
1
0.9 0.8 0.7
0.2 0
0.6
5 mm
0.5
(a)
(b)
Fig. 2. A projection image of the sample (a), with its profile (b) showing significant white noise at the detector level. The liquids occupy the upper portions of the tubes, and air bubbles occupy the lower portions of the tubes.
#212233 - $15.00 USD (C) 2014 OSA
Received 20 May 2014; revised 2 Jul 2014; accepted 9 Jul 2014; published 17 Jul 2014 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018020 | OPTICS EXPRESS 18024
−10
−10
x 10
x 10 6
6
2 0 −2 −4
Absorption coefficient (β)
4 4
2 0 −2 −4
−6 0
−6
5 mm −8
−8
(a)
(c) −7
−7
x 10
x 10 0.5
0
−1 −1.5 −2
0
5 mm
Real refractive index (δ)
0 −0.5
−0.5 −1 −1.5 −2 −2.5
−2.5
−3
(b)
(d)
Fig. 3. Reconstruction results. (a) Filtered backprojection, (b) compressive PAD TIE, (c) plot along line in (a), and (d) plot along line in (b). These reconstructions show cross sections of the Eppendorf tubes containing the liquids taken along the plane coinciding with the line in Fig. 2(a)
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 #212233 - $15.00 USD (C) 2014 OSA
Received 20 May 2014; revised 2 Jul 2014; accepted 9 Jul 2014; published 17 Jul 2014 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018020 | OPTICS EXPRESS 18025
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.
#212233 - $15.00 USD (C) 2014 OSA
Received 20 May 2014; revised 2 Jul 2014; accepted 9 Jul 2014; published 17 Jul 2014 28 July 2014 | Vol. 22, No. 15 | DOI:10.1364/OE.22.018020 | OPTICS EXPRESS 18026
