Compressive X–ray phase tomography based on the transport of intensity equation
Abstract
We develop and implement a compressive reconstruction method for tomographic recovery of refractive index distribution for weakly attenuating objects in a microfocus X–ray system. This is achieved through the development of a discretized operator modeling both the transport of intensity equation and X–ray transform that is suitable for iterative reconstruction techniques.
Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
RMD, Inc., Watertown, MA, USA
Singapore-MIT Alliance for Research and Technology (SMART) Centre, Singapore 138602, Singapore
Corresponding author:
Traditional tomography with hard X–rays recovers the attenuation of an object. Attenuation does not always provide good contrast when imaging objects made of materials with low electron density, e.g. soft tissues. In these cases, richer information is often contained in the phase, i.e. the optical thickness of the sample [1, 2, 3]. Propagation based techniques are particularly suitable for X–ray phase imaging because they allow phase to be recovered from intensity images taken at multiple propagation distances without the need for optical elements [4, 5]. Here we adopt the transport of intensity equation (TIE) which relates the measured intensity to the Laplacian of the phase under a weakly–attenuating sample approximation. Implementing TIE at many angles while rotating the object allows tomographic reconstruction of the refractive index distribution.
TIE tomographic reconstruction first requires a suitable forward model, consisting of 1) a projection of refractive index through the sample and 2) modeling diffraction after the sample with the TIE. Recovering the phase then amounts to inverting these operations on the measured data. A straightforward method of reconstruction is to invert the forward model in two steps [7]. For the first step, the TIE can be solved by a Poisson equation solver. The TIE is ill–posed because the transfer function relating the intensity measurement to the phase tends to zero as the spatial frequency decreases. As a result, reconstructions are often corrupted by significant low–frequency noise, requiring regularization. Tikhonov regularization is most commonly employed to reduce these artifacts [8, 7]. For the second step, a standard tomographic reconstruction is carried out, e.g. using the filtered back–projection (FBP) method. In order for FBP to yield a result free from high–frequency “streaking” artifacts, projections must be taken at many angles. It is often desirable to use fewer projections in order to reduce dose or acquisition time, in which case the tomographic inversion problem is underdetermined. Rather than solve these two inverse problems independently, the forward and inverse models may be adapted into a single–step operation combining TIE and tomography [6, 9, 10, 11]. However, a single step inversion still requires many projection angles in order to avoid artifacts. Iterative solvers have recently been proposed to reduce these artifacts when attempting a reconstruction from a small number of projections. Myers et. al. propose inversion of TIE tomography measurements to obtain a sample distribution using prior knowledge that the sample consists of a single material of known refractive index [12]. Sidky et. al. propose inverting only the tomographic measurement to obtain a boundary–enhanced image [13]. A similar approach has also been proposed to combine a different phase retrieval technique [contrast transfer function (CTF)] and algebraic tomographic reconstruction [14]. In this Letter, we design a forward model that combines TIE and tomography operations in a single, discretized linear operator in the Fourier domain and develop an iterative reconstruction method for recovery of refractive index of a weakly–attenuating object.
A schematic diagram of the imaging geometry is shown in Fig. 1. A point source with mean wavelength is located at the plane of . This is a good approximation for a table–top microfocus X–ray source. The (weakly attenuating) sample is characterized by a real–valued refractive index , where the origin of the cartesian coordinates is located within the object. We further assume that the the object is small enough that the beam passing through it can be approximated as a plane wave oriented along the axis and that the interaction between the sample and the field can be treated using the projection approximation, i.e. phase delay imparted upon the field passing through the sample is . The intensity is recorded by an area detector located at . Although the TIE for microfocus sources is usually formulated using two measurements of intensity at different positions along the optical axis [8], because the sample is assumed to be weakly attenutating, we consider a fixed detector position with two images taken with and without the sample in place, and , respectively. Under the paraxial and small–wavelength approximations, the relationship between and is given by \be g(x,y;θ) ≡ 2πλd’[1-I(M\rsx, M\rsy)I\mi(M\rsx, M\rsy)]= \delxy^2ϕ(x,y;θ), \eewhere is the gradient operator in the plane, is a magnification factor, and is the effective propagation distance. Equation Compressive X–ray phase tomography based on the transport of intensity equation is a finite difference form of the TIE that uses two images taken with and without the sample in place in order to recover the projected phase of a weakly attenuating sample. Taken together, Eqs. (Compressive X–ray phase tomography based on the transport of intensity equation) and (1) specify the continuous model describing how the refractive index of a weakly attenuating sample creates visible signature in the intensity under the assumption that the object subtends a small solid angle as viewed from the source. In a tomographic measurement, the sample is rotated about the –axis such that the projected phase for a given rotation angle is \beϕ(x,y;θ)=∬n(x,y_\rs,z_\rs)δ(y-y_\rscosθ+z_\rssinθ)\mdy_\rs\mdz_\rs,. \ee
In order to construct an iterative process for the retrieval of phase from measured intensity, one must first represent these two operations in a suitable discretized form. Let us assume our detector consists of an grid of square pixels of side length and let denote the number of angular projections. Then the measured at all projection angles may be arranged into a real–valued vector of length . The refractive index of the object may be discretized into a 3D volume consisting of cubic voxels of side length and arranged into a vector . Let denote an matrix representing the discretized form of Eq. (Compressive X–ray phase tomography based on the transport of intensity equation) and denote the matrix representing Eq. (1) such that \be\bg= \bP\bR\bn≡A \bn, \eewhere represents a matrix combining both and into a single operation.
An effective approach to invert this equation and solve for if our data is undersampled is to adapt compressed sensing theory by assuming can be expressed as sparse, i.e. it contains only a small number of nonzero coefficients in some specified basis [15]. Since the sample of interest often consists of regions with constant refractive indices, we choose the following compressive reconstruction model \be^\bn=arg min_\bn∥\bn∥_\TV such that \bg= A\bn, \eewhere the TV function is our sparsity basis, defined as , and , , and are the finite difference operators in the three Cartesian spatial dimensions. We adapt the two–step iterative shrinkage/thresholding algorithm (TwIST) [16] to solve the minimization.
Although in its final form takes the vector to the measurement vector , it is more natural to think about construction of in terms of and arranged as matrices of dimensions and , respectively, as illustrated in Fig. 2. For notational simplicity, rather than referring to the entries of these matrices by their indices, we use the discretized values of position or spatial frequency associated with that entry. In what follows we denote the 2D discrete Fourier transform (DFT) matrix which applies a DFT from the Cartesian coordinates and to spatial frequencies and as and its inverse DFT (IDFT) as .
The discrete projection operator can be implemented in either the spatial domain [12, 13, 14] or the Fourier domain [17]. We adopt the Fourier domain method here since it is robust to discretization error and noise, and computationally much more efficient [18]. The Fourier transform of each projection can be computed using the projection–slice theorem: 1) is applied to , 2) radial slices through the –axis, each corresponding to a given projection angle, are taken in the Fourier domain by an operator , 3) the spatial distributions of projected phase over all projection angles are obtained by a 1D DFT over these slices from the (radial) spatial frequency coordinate to the spatial coordinate , an operation denoted by . The steps in this construction are illustrated graphically in Fig. 2(a).
The TIE operator is simply the discrete Laplacian operator acting on the projected phases, which can be implemented efficiently in the Fourier domain as \be\bP= \delxy^2 = \FT_x,y^-1 H∘\FT_x,y, \eewhere the transfer function is an matrix with each slice along the dimension identical to , and denotes entry–wise multiplication.
There is a subtlety regarding the implementation of these discrete operators to ensure that the output has the correct dimensionality and sample spacing, and that the operator produces results that are stable under iteration. First, recall that is constructed from measurements over two spatial dimensions at a set of projection angles such that the spatial coordinates are sampled over a Cartesian grid with sample spacing . Recall also that we define over voxels of side length . All DFTs and IDFTs are implemented over Cartesian grids in space or spatial frequencies of sample spacing and , respectively, which allows implementation via fast Fourier transforms (FFT). However, attention to resampling is required when implementing the operator . In order to construct a grid over from planar slices through data on the grid , must incorporate interpolation and then resampling (gridding). We implement this by adapting the non-uniform FFT (NUFFT) algorithm [17] such that the NUFFT operator . The TIE and projection operators may then be combined and simplified (eliminating a 1D DFT–IDFT pair) as \beA = \FT_x,y^-1 H∘\FT_x\FT_NU. \ee produces an output with the proper grid spacing and ensures the stability required for the iterative algorithms utilized in compressive reconstruction.
A microfocus source (Hamamatsu L8121–03) located at m, was operated at 40kVp and 100mA to produce a circular focal spot of 5m in diameter. The resulting X–ray beam has central wavelength nm. Intensity images were taken with a custom designed EMCCD based X–ray camera, where a 150m thick RMD CsI:Tl scintillator was fiber optically coupled to an EMCCD [Andor, ixon, 512512 pixels (), 16m pixel size] with a 6:1 taper [19]; the effective pixel size at the scintillator is 96m. The scintillator was at m (, m). For a beetle sample, a single image was taken at every 5 degrees () with 6.7 seconds exposure time during the tomographic measurement,. The discretized sample therefore consists of voxels, while the data vector has entries. The intensity of the incident beam was calibrated by taking a single background image without the sample in place. A single projection measurement of is shown in Fig. 3, which shows the range of attenuation through the sample. Notice that the maximum attenuation of the incident light is approximately of its initial value, which is not strictly negligible. However, it has been previously observed that phase retrieval using the weak–attenuation formula is quite robust to values of attenuation even beyond those seen in this sample [9, 12].
During the reconstruction, we first compute for each angle. Reconstruction results using two different methods are compared in Fig. 4 for the refractive index deviation from the surrounding air . In (a), the Fourier domain TIE solver with Tikhonov regularization chosen to provide optimal results is used to compute phase projections at each angle, and then the FBP method with a Ram–Lak filter is applied for the tomographic inversion. The results suffer from severe streaking artifacts due to missing samples between slices in the Fourier domain. Low–frequency artifacts (blurring) around edges are also observable but are less severe as compared to a single projected phase reconstruction. This is likely due to denser sampling around the origin resulting from the intersection of the Fourier slices. Both artifacts can be greatly suppressed using compressive reconstruction with TV minimization, Eq. (1), whose results are shown in (b), since TV minimization favors large structures with sharp edges. Plots of the along the dashed lines in (a) and (b) are illustrated in Fig. 4(c) for the FBP method and in (d) for the compressive reconstruction. Notive the significant reduction in high frequency artifacts in from (c) to (d). 3D renderings of the compressive reconstruction of the refractive index are shown in Fig. 4(c) and (d).
Compressed sensing works best when the measurement is “incoherent,” i.e. the sparse information in the unknown is evenly spread out in the measurement [15]. This is achieved in our model by the projection operator . The incoherence of the measurement could be improved through the use of source coding [20] or coded aperture [21] which is beyond the scope of our current work.
This work was supported by the Department of Homeland Security, Science and Technology Directorate through contract HSHQDC-11-C-00083. We are grateful to David J. Brady, Justin Lee and Rajiv Gupta for helpful discussions.
References
- [1] A. Momose, T. Takeda, Y. Itai, and K. Hirano, Nat. Med. 2, 473 (1996).
- [2] F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, Nat. Phys. 2, 258 (2006).
- [3] T. Davis, D. Gao, T. Gureyev, A. Stevenson, and S. Wilkins, Nature 373, 595 (1995).
- [4] J. P. Guigay, M. Langer, R. Boistel, and P. Cloetens, Opt. Lett. 32, 1617 (2007).
- [5] R. C. Chen, L. Rigon, and R. Longo, Opt. Express 21, 7384 (2013).
- [6] A. V. Bronnikov, J. Opt. Soc. Am. A 19, 472 (2002).
- [7] A. Burvall, U. Lundström, P. A. C. Takman, D. H. Larsson, and H. M. Hertz, Opt. Express 19, 10359 (2011).
- [8] T. E. Gureyev and S. W. Wilkins, J. Opt. Soc. Am. A 15, 579 (1998).
- [9] A. Groso, R. Abela, and M. Stampanoni, Opt. Express 14, 8103 (2006).
- [10] T. E. Gureyev, G. R. Myers, Y. I. Nesterets, D. M. Paganin, K. M. Pavlov, and S. W. Wilkins, in Proc. SPIE, Vol. 6318, 63180V (2006).
- [11] T. W. Baillie, T. E. Gureyev, J. A. Schmalz, and K. M. Pavlov, Opt. Express 20, 16913 (2012).
- [12] G. R. Myers, D. M. Paganin, T. E. Gureyev, and S. C. Mayo, Opt. Express 16, 908 (2008).
- [13] E. Y. Sidky, M. A. Anastasio, and X. Pan, Opt. Express 18, 10404 (2010).
- [14] A. Kostenko, K. J. Batenburg, A. King, S. E. Offerman, and L. J. van Vliet, Opt. Express 21, 12185 (2013).
- [15] E. Candès, J. Romberg, and T. Tao, IEEE T. Inform. Theory 52, 489 (2006).
- [16] J. M. Bioucas-Dias and M. A. T. Figueiredo, IEEE T. Image Process 16, 2992 (2007).
- [17] J. Fessler and B. Sutton, IEEE T. Signal Proces 51, 560 (2003).
- [18] S. Matej, J. A. Fessler, and I. G. Kazantsev, IEEE T. Med. Imaging 23, 401 (2004).
- [19] V. V. Nagarkar, I. Shestakova, V. Gaysinskiy, B. Singh, B. W. Miller, and H. Bradford Barber, Nucl. Instrum. Methods A 563, 45 (2006).
- [20] J. Petruccelli, L. Tian, and G. Barbastathis, in Computational Optical Sensing and Imaging (Optical Society of America, 2013), paper CTu2C.3.
- [21] K. MacCabe, K. Krishnamurthy, A. Chawla, D. Marks, E. Samei, and D. Brady, Opt. Express 20, 16310 (2012).