1. Introduction

X-ray Phase-Contrast Tomography (PCT) describes a wide range of techniques for reconstructing the three-dimensional distribution of complex refractive index in a sample from measurements that display some form of phase contrast [1–24]. PCT techniques resemble conventional absorption-based computed tomography (CT) techniques in that measurements are taken as the X-ray source traverses a (typically circular) path about the object. A method, which may be as simple as free space propagation of the transmitted X-ray beam as it passes from the object to the detector, is employed at each projection such that the recorded image(s) display both absorption and phase contrast. Consequently, the reconstruction process is not solely reliant on absorption, but also exploits the refraction that occurs as the X-ray beam passes through the object. Such use of phase contrast allows imaging of structures which display little or no absorption contrast, such as the interface between different types of soft tissue [25,26].

Henceforth, we restrict ourselves to discussion of propagation-based PCT with paraxial illumination, in which phase contrast is achieved through free-space propagation between the object and detector [2,3,6–8,11–13,15,18,22,25,27]. The transport of intensity equation (TIE) can be used to reconstruct the phase in the object plane (Fig.1) from two images taken at different propagation distances, provided the difference between these distances is sufficiently small (see e.g. [18], together with references therein). A “monomorphous” object, such as an object composed of a single material [11,28], or an object composed entirely of light (Z<10) elements being imaged with high (60–500 keV) energy X-rays [13], is defined to be an object in which the imaginary part of the complex refractive index is proportional to the real part of the refractive index decrement, with the proportionality coefficient being independent of position inside the object. If we assume a monomorphous object, and assume that the interaction between the radiation and the object can be described by the projection approximation, the phase and the logarithm of intensity in the object plane are proportional to each other. Thus, a reduced amount of information is imprinted on the transmitted wavefield, and phase retrieval may be performed from a single propagated image at each sample orientation [28–30]. Under the further assumption that the object is weakly absorbing, it has been shown (and experimentally verified) that the combination of phase retrieval using the TIE at each projection and cone-beam CT reconstruction using the well-known Feldkamp-Davis-Kress (FDK) algorithm [31,32] yields a phase-and-amplitude contrast tomography (PACT) algorithm for the reconstruction of the full three-dimensional complex refractive index of a weakly absorbing monomorphous object illuminated by paraxial, divergent, partially-coherent polychromatic X-rays, such as those from a lab-based source [19]. Note that, if the radiation is monochromatic, the assumption of weak absorption may be dropped [11, 28]. Furthermore, for both monochromatic and polychromatic illumination, the CT reconstruction process is stable with respect to high frequency noise, unlike conventional absorption-based CT [8,19,22,33]. As is typically the case in propagation-based PCT of single-material objects, strong a priori information regarding the object composition is used to simplify the phase retrieval stage, but no use is made of this information to seek similar simplifications in the second, tomographic stage of the reconstruction.

Typically, hundreds of projections are required to achieve a PCT reconstruction at a useful resolution [1]. We previously noted that single-material (i.e. “binary”) objects, namely objects composed of a single material with known refractive index but unknown spatial distribution, are a subset of monomorphous objects. It is well established in the context of absorption-based tomography that by exploiting the single-material nature of the object, a reconstruction can be obtained from far fewer projections than is otherwise possible, with the resolution of the reconstruction being chiefly limited by the resolution of each projection rather than the number of projections [34–43]. The uniqueness and stability properties of this binary tomography (BT) problem are well studied [44–48], and algorithms have been published which use a variety of methods to solve the problem [34–43, 49].

In this paper we derive an algorithm for binary cone-beam polychromatic PACT using few projections, and verify through both numerical and experimental testing that this exploitation of a priori knowledge in the tomographic stage of the reconstruction allows significant reductions in the scanning time and absorbed dose when compared to existing PACT methods. In section 2 we use an existing PACT algorithm to motivate the union of BT and phase retrieval, and derive an iterative algorithm for Binary PACT. This algorithm is based on existing iterative filtered backprojection (IFBP) algorithms for many-material tomography. Section 3 contains numerical simulations demonstrating the applicability of our algorithm to a variety of objects. These simulations are then used to establish a recommended minimum number of projections for the subsequent experiment, and to investigate the effects of noise on the reconstruction. Finally, in section 4 we present an experimental reconstruction of nylon fibers originally performed from a dataset of 720 images using a PACT algorithm [19], this time using a Binary PACT algorithm and only 20 of the original 720 images. The methods used in this paper, and BT methods in general, generalize easily to the case of multiple-material objects for which the projected thickness of each material can be retrieved in the object plane. This approach to tomographic reconstruction is quite different from, and complementary to, conventional CT.

2. An algorithm for binary phase and amplitude computed tomography

We consider an object with complex refractive e index n(r,λ)=1-Δ(r,λ)+iβ (r,λ) and assume it to be monomorphous, i.e. we assume that the ratio ε(λ)=β(r,λ)/Δ(r,λ) is independent of position. Note that r=(r ₁,r ₂,r ₃) is a Cartesian coordinate system (see Fig. 1), λ is the wavelength, and both β(r,λ) and Δ(r,λ) are real. We assume that the object is located entirely within a sphere of radius d centered at the origin of r. As we are chiefly interested in samples which produce insufficient absorption contrast for conventional tomography, we shall further assume that the object is weakly absorbing, and hence that ε(λ)≪1. For a point source of monochromatic scalar X-rays, we have the following expression for the phase shift ϕ ₀(x,θ,λ) due to the object over the object plane, obtained under the projection approximation [18]:

Here, k=2π/λ is the wavenumber, x=(x ₁,x ₂)=(-R ₁ p ₁/p ₂,R ₁ p ₃/p ₂) is a coordinate system in the object plane, and p=(p ₁,p ₂,p ₃)=(r ₂cosθ-r ₁sinθ,ρ-r ₁cosθ-r ₂sinθ,r ₃) is a rotated coordinate system centered at the X-ray source position s(θ)=(ρcosθ,ρsinθ,0).

The radius of the circle that the source traverses in the r ₁-r ₂ plane is ρ, and R ₁=ρ+d is the distance from the source to the object plane. Note that the source and detector are both considered to move in synchrony, with the object being fixed; thus r=(r ₁,r ₂,r ₃) is a fixed coordinate system, with both x=(x ₁,x ₂) and p=(p ₁,p ₂,p ₃) being moving coordinate systems.

 

Fig 1. The experimental setup.

Download Full Size | PDF

Under these conditions the refractive index can be retrieved from the set of object-plane phase maps according to [31,32]:

where F ₁ is the one dimensional Fourier transform with ξ ₁ dual to x ₁, and B denotes the paraxial divergent-beam filtered backprojection operator. Equation 2 is the Feldkamp-Davis-Kress (FDK) algorithm for paraxial incident illumination. Implicit in the assumption of paraxial incident radiation is the assumption that terms of order d ²/R ² ₁ are negligible. In the absence of this assumption, the object-plane phase is multiplied by the factor R ₁(R ² ₁+x ² ₁+x ² ₂)^-1/2, accounting for the curvature of the wavefront in the object plane. This term can be approximated according to R ₁(R ² ₁+x ² ₁+x ² ₂)^-1/2≈1 for paraxial incident radiation. Letting each projection be sampled over a q×q pixel Cartesian lattice, it is well known that the FDK algorithm requires a number of projections p at least of order πq/2 in order to produce an accurate reconstruction with resolution q×q×q [1].

As we have assumed the absorption contrast in the images to be small, we introduce a small but non-zero propagation distance R ₂ between the object and detector planes. This propagation distance allows the contrast of the collected projection images to be enhanced in comparison to that of the corresponding contact images in the object plane, through the mechanism of propagation-based phase contrast [25]. The object-to-detector propagation distance should be sufficiently small such that the complex amplitude $U_0(\mathbf{x},\lambda )=\sqrt{I_0(\mathbf{x},\lambda )}\exp \left[i{\varphi }_0(\mathbf{x},\lambda )\right]$ in the object plane satisfies the smoothness requirement |(∇² _⊥)ⁿU ₀(x,λ)|≪(2k/R′)ⁿ|U ₀(x,λ)|, where n=1,2,3…, M=(R ₁+R ₂)/R ₁ is the magnification and R′=R ₂/M. Fresnel diffraction causes the transverse phase distribution in the object plane to manifest as transverse intensity fluctuations in the propagated image over the detector plane. The object-plane phase can then be retrieved from a single propagated image using the following variant of the transport of intensity equation [28,29,30]:

Here, $S_{R_2}(M\mathbf{x},\theta ,\lambda )$ is the spectral density in the detector plane, S_in(λ) is the spectral density incident on the object, ξ ²=ξ ² ₁+ξ ² ₂, F ₂ is the two dimensional Fourier transform with (ξ ₁,ξ ₂) dual to x=(x ₁,x ₂) and P_M(x)=P _source(x/(1-M))*P _detector(x), the PSF of the imaging system which accounts for both the finite size of the fully incoherent source and the spatial resolution of the detector (note that an asterisk denotes two-dimensional convolution). Making no further use of the single material nature of the object, Eqs (2) and (3) allow one to derive the following reconstruction formula, which was presented in reference [19]:

In obtaining this equation, we assumed that: (i) $d<\sqrt{\frac{\pi \lambda R\prime }{\epsilon \left(\lambda \right)}}$ ; (ii) the wavelength spread is sufficiently narrow that ε(λ)/λ≅ε(λ ₀)/λ ₀, where λ ₀ is some suitable central or characteristic wavelength of the incident beam. Also, in this equation we have introduced the spectrally averaged refractive index decrement Δ_s(x)=∫^∞ 0 S_in(λ)D(λ)Δ(x,λ)dλ/I_in, the polychromatic (time-averaged) incident detected intensity I_in=∫^∞ ₀ S_in(λ)D(λ)dλ, and the registered intensity $I_{R_2}(\mathbf{x},\theta )={\int }_0^{\infty }{S}_{R_2}(\mathbf{x},\theta ,\lambda )D\left(\lambda \right)d\lambda$ in the detector plane, where D(λ) is the spectral sensitivity of the detector. Equation (4) is a PACT algorithm which allows us to reconstruct the complex refractive index of a weakly absorbing, monomorphous object illuminated by divergent, polychromatic, partially-coherent X-rays, from a single propagated image at each projection [19]. For the purposes of the present paper, this result forms a convenient benchmark against which to compare our Binary PACT algorithm.

In deriving Eq. (4) we made use of the assumption of a weakly absorbing monomorphous object to linearize the phase retrieval process (as in Eq. (3)), and allow phase retrieval from a single image (rather than a pair of propagated images) at each projection. However, the paraxial FDK algorithm (Eq. (2)) used as a basis for the polychromatic, divergent-illumination, PACT algorithm (Eq. (4)) makes no use of this a priori information about the object. It is therefore natural to seek a Binary PACT algorithm that makes full use of this information, both for phase retrieval and tomography.

As a starting point from which to derive such an algorithm, we return to Eq. (3), note that in the single material case ϕ ₀(x,θ,λ)=-kΔ(λ)T(x,θ) [18], where T(x,θ) is the projected thickness, integrate with respect to λ and re-arrange to give:

This formula allows the recovery of the projected thickness of a single material object with refractive index n(λ)=1-Δ(λ)+iβ(λ) from a single propagated image [50] (cf. [28]). Note that μ̄ and $\overline{\Delta }$ are assumed to be known quantities.

Having retrieved the projected thickness at each projection, we can exploit the a priori knowledge that the material is either present or absent at each voxel to produce a reconstruction from far fewer projections than would otherwise be required, breaking the p≈πq/2 “rule of thumb” that applies to conventional CT [34–43,49]. Having first retrieved the projected thickness at each projection using Eq. (5), we define the function $\tilde{\Delta }$ ₀(r)=(B T)(r)$\overline{\Delta }$ . An initial estimate for the spectrally averaged real part of the refractive index decrement Δ_S,0(r) is then determined via a thresholding process:

where M is chosen such that the correct total amount of material is present in the object, i.e. such that the correct number of pixels assume the value Δ_S,0(r)=Δ̄, and T is the thresholding operator. The correct amount of material can either be known a priori, or determined by integrating the projected thicknesses. We then iterate according to:

where P is the projection operator that generates a set of projected thicknesses from the estimated spectrally averaged real part of the refractive index decrement Δ_S,i(r). The parameter 0<γ_i≤1 is varied to ensure a convergent reconstruction as follows: (i) Each iteration is performed first with γ_i=1, and then tested to check that:

thus ensuring that at each step the projected thicknesses (PΔ_S,i+1)(x_l,y_j,θ_k) are closer to the retrieved thicknesses T (x_l,y_j,θ_k) than at the previous iteration. (ii) Should the test fail, the alterations made this iteration are rejected, γ_i is halved, and the iteration is attempted again. (iii) The process terminates when either (a) no successful step is possible without reducing γ_i below some pre-determined limit; or (b) the simulated phase-contrast image data, synthesized using the current iterate of the reconstructed object, is deemed to be sufficiently close to the measured phase-contrast data (i.e. within one standard deviation of the known noise level). The iterative process (Eqs (6)–(8)) can be applied to any case in which the projected thickness T (x, θ) of a material can be found at each projection. Note that unlike in many published BT algorithms, Eqs (6)–(8) make no prior assumptions regarding the smoothness or connectivity of the object, except those required for the validity of Eq. (3).

3. Reconstructions from simulated data

We now present reconstructions performed on simulated data generated for three 2mm×2mm phantoms (see Fig. 2), representing nylon fibers, a vascular structure, and a sparse splatter pattern. The simulations assume a coherent plane wave with wavelength 0.148 nm incident on a binary object with complex refractive index n=1-Δ+iβ, where Δ=6.5×10^-6 and β=1.3×10^-8. This refractive index has been chosen to approximate that of Carbon, a typical “light”, weakly-absorbing material. The “detector” is assumed to lie 1cm away from the object plane for all phantoms.

 

Fig. 2. Phantoms used for the numerical simulations: a vascular structure (left), a splatter pattern (centre), and a nylon phantom (right). Each phantom is 2 mm×2 mm=600×600 pixels in size.

Download Full Size | PDF

 

Fig. 3. Reconstructions with 0 (left), 5 (centre) and 10 (right) percent noise for each of the three phantoms. From left to right, the errors in the reconstructions for the nylon phantom were 0.0079, 0.0095, 0.012. The errors in the vascular reconstructions are 0.014, 0.027 and 0.037, and the errors in the sparse splatter reconstructions are 0.0015, 0.0047 and 0.0073, again reading left to right.

Download Full Size | PDF

Noise is added using a Gaussian distribution with standard deviation equal to a set percentage of the average signal at the detector. Figure 3 shows reconstructions at several levels of noise for each of the simulated phantoms at 600×600 pixel resolution, assuming a detector and phantom pixel size of 3.3 µm, and 20 equally spaced projections. The propagation of the wavefield from the object plane to the detector plane was modelled using a fast-Fourier-transform-based implementation of the angular spectrum formulation of wavefield diffraction, which furnishes a rigorous solution to the Helmholtz equation for a forward-propagating complex monochromatic scalar wavefield, as it travels from plane to parallel plane _{in vacuo} (see e.g. [18]). At the spatial resolution adopted for our simulations, the vessels in the vascular pattern are ~10 or ~15 pixels wide, and the smallest spots in the splatter pattern consist of 9 pixels in total. We express the error as the ratio of the number of incorrectly determined pixels to the total number of pixels in the image. This error metric has been selected as it makes no distinction between pixels based on whether the correct binary state is “on” or “off”. When these simulations were repeated in the absence of noise, at resolutions varying from 100×100 to 600×600 pixels (with correspondingly coarser detector resolutions and the same number of projections), the edges of the retrieved object deviated from those in the original phantoms by two or fewer pixels, resulting in lower calculated errors at higher resolutions. Figure 4 illustrates the behavior of the error with respect to the number of projections for the three phantoms. As expected, the error decreases monotonically as the number of projections is increased. Also note that for low (<12) numbers of projections the behavior of the error depends strongly on object geometry, and high error leads to unreliable reconstructions in the case of the vascular phantom. In the absence of noise, 20 projections were found to be sufficient to lower the error below 1.5% for all phantoms. Consequently, 20 projections were used in the experimental data set in section 4.

 

Fig 4. Error vs. number of projections. The nylon, vascular and splatter phantoms are represented by the solid, dashed, and dotted lines respectively.

Download Full Size | PDF

4. Reconstructions from experimental data

An experimental tomographic data set was collected using a laboratory-based X-ray source, the X-ray ultra microscope (XuM) [11]. Focusing a 30 keV electron beam onto a 500 nm thick tantalum foil target produced a divergent polychromatic X-ray beam with a source size of 0.2 µm, and a mean photon energy of 8.4 keV. The sample consisted of four nylon wires 100 µm, 240 µm, 330 µm and 420 µm in diameter, enclosed by a Perspex cylindrical shell with 2 mm outer diameter and 100 µm thick walls. The nylon wires were strung such that they were parallel to the axis of the cylinder. Using a source-to-sample distance R ₁=25mm and source-to-detector distance ce R ₁+R ₂=259mm, 720 images of the sample were obtained, with a total exposure time of more than 15 hours. The magnification (10.4×) present in this imaging system as a consequence of the free-space source-to-detector propagation increases the effective spatial resolution. The images were equally spaced over the range θ∈[0,2π), and were collected using a direct detection, deep depletion CCD. The X-rays had a mean wavelength of λ ₀=1.48Å (corresponding to 8.4 keV), so that, for both nylon and Perspex, ε(λ ₀)≅0.0016. Note that the quasi-monochromatic approximation was inapplicable due to the presence of significant bremsstrahlung radiation in the X-ray beam [11]. Figure 5 is a sample image from the data set. Fringing due to the propagation induced phase-contrast is clearly visible around the edges in the image, indicating that a reconstruction performed using conventional absorption-based CT methods will contain so-called “edge enhancement” artifacts.

 

Fig. 5. A phase-contrast image typical of those in the experimental data set.

Download Full Size | PDF

Figure 6 shows a two dimensional cross-section through a numerical reconstruction of β(r,λ ₀) performed with the PACT algorithm (Eq. (4)), utilizing the full data set of 720 projections. Figure 7 shows a Binary PACT reconstruction of β(r,λ ₀) from 20 projections equally spaced across the interval θ∈[0,π), using a numerical implementation of Eq. (7), with ε(λ ₀)=0.0016 (λ ₀=1.48Å).

 

Fig. 6. Reconstruction from experimental data using the PACT algorithm (720 projections).

Download Full Size | PDF

Note that the Perspex tube was made of a different material to the nylon phantom it contained. In order to ensure that the data satisfied the assumption of a binary object, the signal resulting from the Perspex tube was simulated and subtracted from the measured data at each projection. Subsequently, we reconstructed the binary nylon phantom contained within the tube. The spatial resolution of both reconstructions was 3.86µm per cubic voxel. The reduction in the number of projections used lowered the scanning time required to collect the data set from ~15 hours to ~25 minutes, and required 1/36^th of the absorbed dose. Figure 8 shows one dimensional cross-sections through both the PACT and Binary PACT reconstructions. Note that the sharp edges produced by the Binary PACT algorithm are well within the several-pixel-wide regions corresponding to edges in the PACT reconstruction.

 

Fig. 7. Reconstruction from experimental data using Binary PACT (20 projections).

Download Full Size | PDF

 

Fig. 8. Cross-sections through both the PACT reconstruction from 720 projections (solid line), and the Binary PACT reconstruction from 20 images (dotted line).

Download Full Size | PDF

Note that the observed X-ray attenuation reached almost 50 % in some areas of the sample. It was difficult to avoid significant absorption, as our X-ray microscope becomes less efficient at higher X-ray energies. This means that the weak absorption condition, ε(λ)≪1, was violated in the experiment. Nevertheless, one can see from Figs.7 and 8 that a good quantitative accuracy has been achieved in the reconstruction using the proposed algorithm, which suggests that the phase-retrieval formula Eq. (5) has a broader region of validity than the one indicated by the presented derivation. This conjecture is supported by some recent independent results, e.g. [24, 51].

5. Conclusion

We have derived and tested a method for the high-resolution three-dimensional cone-beam polychromatic imaging of weakly-absorbing binary objects from few projections. This method preserves the beneficial properties of cone-beam PACT, such as increased sensitivity to weakly-absorbing features and increased spatial resolution due to the geometric magnification of the imaging system, and incorporates the significant reductions in scanning time and absorbed dose typically provided by BT methods. In section 2, we presented a method for Binary PACT based on existing methods for iterative filtered backprojection. In section 3 we used computer simulations to demonstrate the stability of this reconstruction algorithm with respect to reasonable levels of noise for several types of model object, and established a recommended minimum number of projections for the subsequent experiment. Finally, section 4 gave an experimental implementation in which we reconstructed the three dimensional complex distribution of refractive index in a binary sample from 20 cone-beam polychromatic phase-contrast images acquired on an X-ray ultra Microscope.

The techniques discussed in this paper suggest two immediate extensions to non-binary objects. The first extension involves an adjustment of the thresholding process to accommodate a two- or three-material object, and the use of the more general, two- or three-image variant of the transport of intensity equation [18]. In the limiting case of an object composed of an infinite number of possible materials, with an unknown concentration of each, this reconstruction process becomes IFBP [52]. Consequently, we expect the required number of projections to increase as we increase the number of materials composing the object. The second option is to further explore the potential for cancellations between the phase retrieval and tomographic portions of the algorithm, for example by extending the phase retrieval technique such that the projected thickness of each material in the object can be retrieved at each projection. This “separation” of materials can already be performed in absorption CT by taking images on either side of absorption edges. This avenue of extension represents an alternative in which more information is retrieved at each angular position, but fewer projections are taken.