## Abstract

In X-ray computed tomography (CT) increased information requirements (e.g. increased resolution) typically lead to a concurrent increase in the required number of viewing angles, scanning time and delivered dose. We demonstrate that using phase-contrast imaging it is possible to “dissect” two- and three-material objects into their component materials, which in combination with binary tomographic techniques allows us to satisfy increased information requirements *without* taking the usual images at additional viewing angles. This imaging scheme reduces the scanning time and dose delivered to samples by at least an order of magnitude when compared to conventional X-ray CT. The effects of noise on our reconstruction scheme are investigated for simulated data. Finally, a slice through a glass tube filled with silica and water is reconstructed from 18 projection images taken on an X-ray ultra Microscope (XuM).

©2008 Optical Society of America

## 1. Introduction

X-ray computed tomography (CT) is a widely used method for the non-destructive, three-dimensional imaging of samples. A vast range of both macro and micro-scale applications exist in fields as diverse as medicine, materials science, biology, etc. Images showing the attenuation of the incident illumination as it passes through the sample are collected as the X-ray source and detector traverse some path about the object [1]. In conventional tomography, these absorption-contrast images are then used to generate a three dimensional map of the linear attenuation coefficient of the sample [1]. Consequently, materials with very similar linear attenuation coefficients are difficult to resolve under realistic imaging conditions [2,3]. Although features composed of such materials do not produce measurable absorption contrast, implementation of a suitable phase-contrast imaging mechanism at each projection angle may allow reconstruction of the complex refractive index of the sample [1,4-27]. In this context “phase-contrast imaging” may be defined as any form of imaging that is sensitive to the real part of the refractive index of a sample. Using such an imaging modality, materials may then be distinguished on the basis of their refractive, as well as attenuative properties. This combination of phase retrieval and tomographic imaging is termed phase-contrast tomography (PCT).

In this paper we consider a subset of PCT methods which employ propagation-based phase retrieval techniques [2,4,5,8-10,14-16,18,21,25,28]. The detector is placed a sufficient distance away from the object such that Fresnel diffraction converts the phase variations in the contact plane (which are a measure of the projected refractive properties of the sample), into recorded intensity variations at the detector (propagation-based phase contrast). Given two such measurements taken at different suitable propagation distances the Transport of Intensity equation (TIE) can then be used to reconstruct the contact-plane intensity and phase [21,29]. It is known that in the case of a single-material object the proportionality between the real and imaginary parts of the refractive index allows phase retrieval to be performed from only a single image taken at each projection angle [30].

The *a priori* knowledge that the object is composed of a single material can also be used to simplify conventional attenuation-based CT reconstruction. Methods for binary tomography (BT) exploit the “on-off” nature of the refractive index of a single-material object, and allow reconstruction from far fewer projection angles than are otherwise required [31-40]. Methods for BT are well explored in the context of absorption tomography [31-41], as are the stability and solubility conditions of the BT reconstruction problem [42-46].

A previous publication explored the natural synergy between methods for single-material phase retrieval and methods for the binary tomography of single-material objects [47]. In this paper we further explore the relationship between propagation-based phase retrieval and binary tomography by showing that phase retrieval, together with appropriate use of *a priori* information, allows methods for binary tomography to be applied to objects composed of two and three materials. This is accomplished by retrieving the projected thickness of each of the component materials as part of the phase retrieval process, effectively “dissecting” a more complex object into two or three interlocking binary objects.

In section 2 we review some relevant work concerning single-material phase retrieval and binary tomography before presenting our algorithms for the tomography of two- and three-material objects. In section 3 we test these methods on simulated data with varying levels of noise. Finally, in section 4 we apply these algorithms to a PCT data set of a glass tube containing silica grains surrounded by water dosed with 10% potassium iodide (KI) by weight. The distribution of silica and water in a slice through the sample is then reconstructed from 18 of the original 180 projections, with attendant dramatic reductions in required scanning time and absorbed dose.

## 2. Algorithms for phase-and-amplitude computed tomography of binary, ternary, and quaternary objects

#### 2.1 A review of algorithms for phase-and-amplitude computed tomography of binary objects

We begin by reviewing the tomography of binary objects. This problem is well explored and serves as a foundation to the related problems of ternary and quaternary tomography.

A binary object is composed of only a single material and void [see Fig. 1(a)] such that the complex refractive index [21] *n*(**r**)=1-Δ(**r**)+*iβ*(**r**) of the object can assume only the values 1-Δ_{1}+*iβ*
_{1} and 1, where **r**=(*r*
_{1},*r*
_{2},*r*
_{3}) is a 3D Cartesian coordinate system used to describe the object and both Δ(**r**) and *β*(**r**) are real. Tomographic reconstruction of a binary object from projection measurements is typically a two-stage process. The first stage in our analysis involves the recovery of the projected thickness

at various projection angles *θ*∈[0,*π*), where **x**=(*x*
_{1}, *x*
_{2}) is a Cartesian coordinate system in the contact plane, *P* denotes the projection operator along the direction of X-ray propagation, and the projection angle *θ* is measured anticlockwise from the *r*
_{1} axis in the *r*
_{1}-*r*
_{2} plane (see Fig. 2).

In the second stage, the 3D distribution of complex refractive index *n*(**r**) is reconstructed from the projected thicknesses *T*
_{1}(**x**,*θ*).

Assuming the projection approximation is valid, when the object is illuminated by monochromatic plane-wave X-rays with incident intensity *I _{in}* and wave-number

*k*, the intensity

*I*

_{0}(

**x**,

*θ*) and phase

*ϕ*

_{0}(

**x**,

*θ*) in the contact plane are related to the projected thickness according to [21]:

Note that under the projection approximation absorption is governed entirely by the imaginary part *β*(**r**) of the complex refractive index, and refraction is governed entirely by the real decrement Δ(**r**). If the object produces strong absorption contrast, then an image taken in the contact plane at a given projection angle allows the trivial recovery of the projected thickness at that projection angle via a simple re-arrangement of Eq. (2):

In cases where absorption contrast in the contact plane is weak, the experimental setup can be altered to include a sufficient propagation distance *R* between the object and detector (see Fig. 2). As the X-rays propagate through this region of free space they undergo Fresnel diffraction. Perturbations in the contact-plane phase manifest as changes in the intensity recorded at the detector. The “downstream” recorded image *I _{R}*(

**x**,

*θ*) then exhibits a degree of propagation-induced phase contrast, allowing imaging to proceed based on the refractive

*and*attenuative properties of the object [21]. Provided the propagation distance

*R*is sufficiently short that the complex amplitude ${U}_{0}(\mathbf{x},\theta )=\sqrt{{I}_{0}(\mathbf{x},\theta )}\mathrm{exp}\left[i{\varphi}_{0}(\mathbf{x},\theta )\right]$ in the contact plane satisfies the smoothness condition [21,48]

where ${{\nabla}_{\perp}}^{2}={{\partial}_{{x}_{1}}}^{2}+{{\partial}_{{x}_{2}}}^{2}$ is the gradient operator in the contact plane, the free-space propagation may be modeled by the following finite difference form of the Transport of Intensity equation (TIE) [29]:

The TIE is the continuity equation associated with the parabolic equation of paraxial wave optics. It is known that for binary objects with some absorption (*i.e.*
*β*
_{1}≠0) the TIE reduces to

upon substitution of Eqs. (2) and (3) into Eq. (6) [30]. Given the propagated intensity distribution at a given projection angle, Eq. (7) may then be solved for the projected thickness using existing numerical techniques such as the full multigrid algorithm [49], or fast Fourier transform (FFT) based methods [21,30]:

where **F**
_{2} is the two-dimensional Fourier transform with (*ξ*
_{1},*ξ*
_{2}) dual to (*x*
_{1},*x*
_{2}) and *ξ*
^{2}=*ξ*
_{1}
^{2}+*ξ*
_{2}
^{2}. The behavior of the retrieved thickness with respect to noise in the recorded intensity is well understood. In particular, the retrieved thickness is stable with respect to high frequency noise in *I _{R}*(

**x**,

*θ*), and the presence of absorption helps mitigate the effects of noise at lower Fourier frequencies [29,30,50]. The stability of the reconstruction with respect to low frequency noise can be further improved by choosing the optimal propagation distance R that satisfies the validity conditions for the TIE (Eq. 5) [29,50,51]. Methods for determining this optimal propagation distance and other, more sophisticated methods for dealing with low frequency noise when using the TIE exist but are beyond the scope of this paper [29,50]. The synergy between single-material phase retrieval algorithms and methods for tomographic reconstruction of binary objects has also been explored for polychromatic, partially-coherent, illumination from a small lab-based source [47].

Having solved for the projected thickness at a variety of projection angles we must then use this information to reconstruct the 3D distribution of complex refractive index within the object. Conventional methods for tomographic reconstruction such as filtered backprojection typically require information at hundreds of projection angles in order to produce a reconstruction at a useful resolution [1]. In contrast, methods developed specifically for “binary tomography” can significantly reduce the required number of projection angles by making use of the *a priori* knowledge that the complex refractive index is either 1-Δ_{1}+*iβ*
_{1} or 1 in any given voxel [31-40]. Note that these methods implicitly assume that the object is binary, and the voxel size is sufficiently small that the percentage (and hence overall effect) of partially-filled voxels (due to edges which do not exactly follow the voxellation scheme) is negligible. In this paper we will use the method for binary tomography presented previously in reference [47]. Though other more sophisticated methods exist, this method has the advantage of being easy to implement numerically, as it incorporates conventional filtered backprojection. This method, unlike most alternative methods, also makes no assumptions about the smoothness or connectivity of the object.

We make an initial estimate $\stackrel{\u02d8}{\Delta}$
_{0}(**r**) of the real decrement to the refractive index using filtered backprojection:

where **F**
_{1} is the one-dimensional Fourier transform and ** B** represents the filtered backprojection operator [1]. Note that the subscript on $\stackrel{\u02d8}{\Delta}$

_{0}(

**r**) denotes that this is the estimated real decrement to the refractive index at the zeroeth iteration of the reconstruction process. We then apply the known constraints to the real refractive index via a thresholding procedure:

where the numeric subscript on $\stackrel{\u02d8}{\Delta}$
_{T, 0}(**r**) again denotes the zeroeth iteration of the reconstruction procedure, and *M* is a constant chosen such that the correct amount of material is present in the object, *i.e.* such that:

for all *θ*. This zeroeth order estimate of the refractive index is then iteratively refined:

At each iteration the parameter 0<*γ _{i}*≤1 is optionally varied to ensure a convergent reconstruction. If necessary, the step is performed for several values of

*γ*, beginning at 1 and halving until a value of

_{i}*γ*is found for which:

_{i}ensuring that at each step the projected thicknesses (** P**
$\stackrel{\u02d8}{\Delta}$

_{T, i+1})(

**x**,

*θ*) are closer to the retrieved thicknesses

*T*

_{1}(

**x**,

*θ*) than at the previous iteration. The process terminates if no successful step is possible without reducing

*γ*below some pre-determined lower limit, or if (

_{i}**$\stackrel{\u02d8}{\Delta}$**

*P*_{T, i+1})(

**x**,

*θ*) is deemed to be sufficiently close to

*T*

_{1}(

**x**,

*θ*) (i.e. within one standard deviation of the known noise level). It has been demonstrated (using both simulated and experimental data) that this method typically requires measurements at 20 or fewer projection angles in order to produce an accurate reconstruction [47].

#### 2.2 Algorithms for phase-and-amplitude computed tomography of ternary objects

In this section we extend the method presented in section 1.1 for 3D imaging of binary objects to objects composed of two materials and void, which we term “ternary” objects (see Fig. 1(b)). The refractive index of a ternary object can assume the values 1-Δ_{1}+*iβ*
_{1}, 1-Δ_{2}+*iβ*
_{2} and 1. Note that we have assigned the subscripts “1” and “2” in order to differentiate between the two materials, and all equations are equally valid if all such subscripts are exchanged. One possible course of action would be to take measurements at additional projection angles and “relax” the second (tomographic) stage of the reconstruction process to accommodate three possible values of the refractive index when thresholding (Eq. 10), rather than two. Instead, we consider our ternary object to be composed of two non-overlapping binary objects (that may be superimposed in projection) such that:

The *first* step of our reconstruction process is then altered to recover the projected thicknesses of materials one and two (*T*
_{1}(**x**,*θ*) and *T*
_{2}(**x**,*θ*) respectively) at each projection angle. Due to the additional degree of freedom of the ternary object problem in comparison to the binary object considered in section 2.1, in the absence of additional *a priori* information phase retrieval forms an integral part of this reconstruction method and is not used solely to compensate for poor absorption contrast as in Eq. (7).

Assuming a ternary object that satisfies the projection approximation, illuminated by monochromatic plane-wave X-rays, the contact plane intensity and phase obey the respective equations:

We record two images at each projection angle, one in the contact plane and one at some small distance *R* downstream. If *R* is again chosen to satisfy the validity condition (Eq. (5)) for the TIE (Eq. (6)), then given knowledge of *I _{R}*(

**x**,

*θ*) and

*I*

_{0}(

**x**,

*θ*) one can solve the TIE for ϕ

_{0}(

**x**,

*θ*) [21, pp. 297-301]. This solution is known to be stable with respect to high frequency noise in the recorded intensities, and relatively unstable with respect to low frequency noise [29,50]. As mentioned previously, methods for dealing with this instability exist, but are beyond the scope of this paper. If experimental constraints make it impractical to measure the intensity in the contact plane one can instead measure the intensity at two propagation distances

*R*

_{1}and

*R*

_{2}=

*R*

_{1}+

*R*, where

*R*again satisfies the validity condition for the TIE. The TIE can then be solved for the propagated phase ${\varphi}_{{R}_{1}}(\mathbf{x},\theta )$ . Given ${I}_{{R}_{1}}(\mathbf{x},\theta )$ and ${\varphi}_{{R}_{1}}(\mathbf{x},\theta )$ , we have complete knowledge of the forward-propagating coherent X-ray wavefield over the plane

*R*

_{1}and can then use the inverse wavefield propagator [21,52-55] to solve for

*I*

_{0}(

**x**,

*θ*) and

*ϕ*

_{0}(

**x**,

*θ*). As we now know both ln[

*I*

_{0}(

**x**,

*θ*)/

*I*] and

_{in}*ϕ*

_{0}(

**x**,

*θ*), knowledge of Δ

_{1}, Δ

_{2},

*β*

_{1}and

*β*

_{2}allows us to solve Eqs. (15) and (16) as a system of linear equations for

*T*

_{1}(

**x**,

*θ*) and

*T*

_{2}(

**x**,

*θ*). This solution will be stable provided the ratios

*β*

_{1}/Δ

_{1}and

*β*

_{2}/Δ

_{2}are sufficiently different (

*i.e.*provided the condition number of the system of linear equations is sufficiently close to unity).

In some cases the total projected thickness of the object:

may be known *a priori*. If the sample is known to contain no internal voids then *A*(**x**,*θ*) can be found by measuring the surface of the object using, for example, laser profilometry [56]. Note that even if there are no voids within the sample, it still constitutes a ternary (rather than binary) object due to the (in general) non-trivial shape of the interface between the object and the surrounding void. Reconstruction may then proceed using only a single image at each projection angle. If images taken in the contact plane display sufficient absorption contrast, then *T*
_{1}(**x**,*θ*) and *T*
_{2}(**x**,*θ*) can be recovered by solving the system of linear Eqs. (15) and (17), and phase-contrast imaging is no longer necessary. If the total projected thickness is known *a priori* but images taken in the contact plane display insufficient absorption contrast, then we again introduce some propagation distance *R* between the contact plane and detector and take a single propagated image at each projection angle. Equations (15)-(17) can be re-arranged into:

which can be substituted into the TIE (Eq. (6)) to give

$$B(\mathbf{x},\theta )=\mathrm{exp}\left[-2k\left({\beta}_{2}-{\mathrm{\Delta}}_{2}\kappa \right)A(\mathbf{x},\theta )\right],$$

$$C(\mathbf{x},\theta )=\mathrm{exp}\left[-2k\left({\beta}_{1}-{\beta}_{2}\right){T}_{1}(\mathbf{x},\theta )-2k{\mathrm{\Delta}}_{2}\kappa A(\mathbf{x},\theta )\right],$$

where *κ*=(*β*
_{1}-*β*
_{2})/(Δ_{1}-Δ_{2}), and *B*(**x**,*θ*) is known. Equation (20) is an elliptical differential equation that can be solved for *C*(**x**,*θ*) using standard numerical methods such as the full multigrid algorithm [49]. Recovery of *T*
_{1}(**x**,*θ*) then proceeds as

Recall that Eqs. (15), (16), (17) and (6) are identical upon the exchange of subscripts 1 and 2, as our numbering scheme for materials one and two is entirely arbitrary. Consequently, Eq. (20) remains correct if all “material” subscripts are exchanged, allowing recovery of *T*
_{2}(**x**,*θ*).

Given knowledge of both *T*
_{1}(**x**,*θ*) and *T*
_{2}(**x**,*θ*), the iterative reconstruction process defined in Eqs. (9)-(13) can be employed to separately reconstruct *n*
_{1}(**r**) and *n*
_{2}(**r**). We have satisfied the increased information requirements for tomography of ternary objects by collecting additional information at each projection angle. Imaging of ternary and quaternary (see section 2.3) objects may thus be performed using measurements at no more projection angles than are required to reconstruct a binary object. This approach is complementary to conventional CT methods, in which the required additional information is typically collected by increasing the number of projection angles and the resolution of the measurements taken at each projection angle [1].

#### 2.3 An algorithm for phase-and-amplitude computed tomography of quaternary objects

The extension of our reconstruction methods from ternary to quaternary objects [see Fig. 1(c)] proceeds as follows. In a quaternary object (*i.e.* an object composed of three materials and void), the refractive index may assume the values 1-Δ_{1}+*iβ*
_{1}, 1-Δ_{2}+*iβ*
_{2}, 1-Δ_{3}+*iβ*
_{3} and 1. Following the path we took in section 2.2, we describe the quaternary object as three non-overlapping binary objects with complex refractive indices *n*
_{1}(**r**), *n*
_{2}(**r**) and *n*
_{3}(**r**) such that:

Under the projection approximation, assuming illumination by monochromatic plane-wave X-rays:

In order to reconstruct the projected thicknesses *T*
_{1}(**x**,*θ*), *T*
_{2}(**x**,*θ*) and *T*
_{3}(**x**,*θ*), we require knowledge of the total projected thickness in addition to images *I*
_{0}(**x**,*θ*) and *I _{R}*(

**x**,

*θ*) taken in the contact plane and at some small distance

*R*downstream at every projection angle. We define the total projected thickness:

We employ the TIE (Eq. (6)) to solve for *ϕ*
_{0}(**x**,*θ*) given *I*
_{0}(**x**,*θ*) and *I _{R}*(

**x**,

*θ*), and then solve the linear system of simultaneous Eqs. (23)-(25) for

*T*

_{1}(

**x**,

*θ*),

*T*

_{2}(

**x**,

*θ*) and

*T*

_{3}(

**x**,

*θ*). Again, the system of linear equations will be stable provided the condition number is sufficiently close to unity. This requirement is violated if any two materials have almost identical complex refractive indices, or if the ratio

*β*/Δ

_{n}*is almost identical for*

_{n}*all three*materials. As in section 2.2, ${I}_{{R}_{1}}(\mathbf{x},\theta )$ and ${I}_{{R}_{2}}(\mathbf{x},\theta )$ can be used where measurement of

*I*

_{0}(

**x**,

*θ*) and

*I*(

_{R}**x**,

*θ*) is not practical. An appropriate method for tomographic reconstruction of binary objects (see, for example, Eqs. (9)-(13)) from projected thicknesses is then employed to reconstruct the 3D distributions of complex refractive index

*n*

_{1}(

**r**),

*n*

_{2}(

**r**) and

*n*

_{3}(

**r**). The retrieved thicknesses will be stable with respect to noise in

*I*

_{0}(

**x**,

*θ*) and

*A*(

**x**,

*θ*) due to the stability of the aforementioned system of linear simultaneous equations, and unstable with respect to low frequency noise in the contrast function

*I*(

_{R}**x**,

*θ*)-

*I*

_{0}(

**x**,

*θ*), in the usual manner for the TIE [29,50]. Phase contrast is again vital to the reconstruction process, as knowledge of the contact plane phase is required to solve the system of linear Eqs. (23)-(25).

As knowledge of two intensity distributions recorded at different propagation distances allows us to calculate the contact plane intensity and phase, and hence complete knowledge of the X-ray wavefield, further intensity measurements at each projection angle contribute no further information and do not allow us to generalize this imaging method to quinary (four material and void) objects unless some other parameter, such as the wavelength of the illuminating X-rays, is altered [57].

#### 3. Reconstruction of simulated ternary and quaternary objects from noisy data

In this section we test the methods presented in sections 2.2 and 2.3 for tomography of ternary and quaternary objects, using simulated data with varying levels of noise. In all cases images are simulated at 20 projection angles about the object, equally spaced over the range *θ*∈[0°,179°].

At each projection angle the contact plane intensity and phase are generated using the projection approximation [Eqs. (18)-(19), (23)-(24)]. When simulated images are required some distance downstream from the contact plane, free-space propagation is modeled using a numerical implementation of the Kirchhoff integral. The noise introduced into the simulated images obeys Poisson statistics. The mean number of photons incident on the detector per pixel is chosen such that the mean standard deviation of the error is equal to a given percentage of the mean intensity. The projected thicknesses of the various materials are reconstructed using methods appropriate to the object: in section 3.1 we use Eq. (20), and in section 3.2 we use the method outlined in section 2.3. A slice (we have chosen *x*
_{2}=0) through the reconstructed thickness distributions is taken, and sequential applications of our method for binary tomography [Eqs. (9)-(13)] are used to reconstruct a 2-D slice through the center of the object. When reconstructing the distribution of each material from the projected thicknesses we constrain our reconstruction algorithm [Eqs. (9)-(13)] so that any voxel already “filled” is off-limits to subsequent materials. Reconstructions are thus performed in order of decreasing *β _{n}* as strongly absorbing materials tend be retrieved more accurately in the presence of noise, and provide a better “guide” to subsequent materials. Simulations are performed for 0, 5, and 10 percent noise, and errors are expressed as the ratio of the total number of incorrect pixels to the total number of correct pixels in the final reconstructed image.

#### 3.1 Simulated ternary object

We assume incident monochromatic plane-wave X-rays of wavelength 0.149 nm. The simulated ternary phantom consists of an irregularly shaped mass with various inclusions, the smallest of which are 3 pixels wide [see Fig. 3(a)]. We examine this phantom for three different compositions: A weakly absorbing mass (carbon) with weakly absorbing inclusions (water), a strongly absorbing mass (aluminum) with weakly absorbing inclusions (Mylar), and a weakly absorbing mass (Mylar) with strongly absorbing inclusions (aluminum). Mylar (C_{10}H_{8}O_{4}) is a hard plastic chosen to (very) roughly approximate human soft tissue as a mix of carbon, oxygen and hydrogen.

Aluminum is used as its absorptive properties are very similar to those of human bone (see, e.g. [58]). Projected thicknesses are retrieved from a single propagated image (*R*=4 mm) and known total projected thickness distributions using Eq. (20). Figure 4 shows the reconstructed distributions of each material for varying compositions and noise levels, and Table 1 the errors in each reconstruction. The gray levels in Fig. 4 are “false color” and are not indicative in any way of the relative refractive index of the two materials. Note that some or all of the three pixel wide inclusions are visible for both the Al/Mylar and Mylar/Al phantoms at all levels of noise. Unsurprisingly, the presence of strong absorption in either material improves the behavior of the reconstructions with respect to noise. We also note that the 10% Poisson noise introduced for the Carbon/H_{2}O phantom is of similar magnitude to the contrast function ln[*I*
_{0}(**x**,*θ*)/*I _{in}*].

#### 3.2 Simulated quaternary object

This phantom [see Fig. 3(b)] is a three-material phantom based on the Shepp-Logan head phantom [59]. Weakly absorbing brain matter of two different densities is surrounded by bone. Incident illumination is assumed to be parallel, with a wavelength of 0.044 nm. At each projection angle two simulated images are produced, one in the contact plane and one 7 m downstream. The phase in the contact plane is then retrieved from these simulated intensity distributions using the TIE [Eq. (6)], and the system of simultaneous Eqs. (22)-(24) then solved for the projected thickness of each of the three materials. Figure 5 shows the reconstructed distributions of each material for varying levels of noise. The 10% Poisson noise is again of similar magnitude to the contrast function. As in Figs. 1, 3 and 4, the gray levels are not indicative of the relative refractive indices of the various materials.

## 4. Reconstruction of a ternary object from experimental data

Passing from simulated to experimental data, X-ray ultra Microscope (XuM) [14] images were collected. The sample was a ~500 µm diameter glass capillary containing sand (silica) grains and water with 10% by weight potassium iodide as a contrast agent. 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 at full-width half-maximum of 0.2 µm, and a mean photon energy of 8.4 keV. The XuM utilized a point projection geometry, with a source-sample distance of 8.09 mm and source-detector distance of 259 mm. Images were collected using a direct detection, deep depletion CCD at 1 degree intervals over the range *θ*∈[0°,179°]. Figure 6 is a sample image from the data set. Fringing due to the propagation induced phase-contrast is visible around the edges in the image (see inset in Fig. 6).

Filtered backprojection of ln[*I _{R}*(

**x**,

*θ*)/

*I*] at all 180 projection angles allows us to approximately reconstruct the three dimensional distribution of imaginary refractive index of the sample. Note that as we possess only a single propagated image per projection angle it is impossible to exactly retrieve the contact plane intensity and phase without some

_{in}*a priori*information. Consequently, the propagation-based phase contrast in the collected images leads to some edge enhancement in the reconstructed slice [see Fig. 7(a)]. As the diameter and composition of the glass tube containing the sample was known, it could be artificially removed from each image, leaving a ternary object composed only of silica, water, and an air bubble.

Single-material phase retrieval (Eq. (8) and binary tomography [Eqs. (9)-(13)] were used to approximately reconstruct the position of the air bubble from images taken at 18 equally spaced projection angles. The refractive index of the “single” material was chosen to be the averaged refractive indices of 10% potassium iodide and silica. Knowledge of the position of the air bubble then leads to knowledge of the total projected thickness, and hence enough information to perform two-material phase-retrieval using Eq. (20). Finally, Eqs. (9)-(13) were used to reconstruct the distributions of silica and water within the tube [see Fig. 7(b)] from data collected at only 18 of the 180 projection angles, equally spaced over the range *θ*∈[0°,179°]. The divergence angle of the incident beam is ~4°, which leads to a maximum lateral shift of 8 pixels at the edge of the sample. However, the slice chosen for reconstruction lies in approximately the same plane as the source path, and the region of interest is towards the center of the sample, where a smaller lateral shift will take place. Note that the divergent nature of the incident illumination and the presence of weak (though non-negligible) bremsstrahlung radiation in the incident beam violates the assumption of monochromatic plane X-rays made in Eqs. (20) and (9)-(13). We thus expect some minor inaccuracies in our reconstruction [Fig. 7(b)], which nevertheless agree well with the “full” reconstruction [Fig. 7(a)]. The original experimental data was quite noisy; it is also worth noting that as Fig. 7(b) only uses 1/10^{th} of the available image set, the effective signal to noise ratio is
$\sqrt{10}\approx 3.2$
times lower for the purposes of tomographic reconstruction.

## 5. Conclusion

We have demonstrated that application of phase-retrieval techniques together with appropriate use of *a priori* knowledge allows us to apply our method for binary tomography to ternary and quaternary objects without the need to collect images at additional projection angles. We have tested the robustness of the reconstruction algorithms derived in section 2 using simulated data with varying levels of noise, and have demonstrated using experimental data that our algorithms allow fast, low-dose, three-dimensional imaging of a ternary sample from few projections. Extension of this method to quinary (four material and void) and higher order objects requires the introduction of some additional degree of freedom into the imaging system. Motivated by several recent papers, we have suggested that this could be achieved by varying the wavelength of the illuminating X-rays. The results in this paper further support the notion that recording additional information at each projection angle can reduce the required number of viewing angles. This in turn suggests that there exists some optimal balance between these two “quantities” which depends on the nature of the *a priori* information one possesses about the object.

## Acknowledgements

DMP acknowledges funding from the Australian Research Council, via the Discovery-Projects programme. GRM, DMP and TEG acknowledge useful discussions with Dr. Karen Siu and Dr. Marcus Kitchen. TEG and SCM wish to thank XRT Ltd. for encouragement of this work.

## References and links

**1. **F. Natterer, *The Mathematics of Computerized Tomography* (Wiley, New York, 1986).

**2. **R. Fitzgerald, “Phase-sensitive x-ray imaging,” Phys. Today **53**, 23–26 (2000). [CrossRef]

**3. **J. Řeháček, Z. Hradil, J. Peřina, S. Pascazio, P. Facchi, and M. Zawisky, “Advanced neutron imaging and sensing,” Adv. Imaging Electron. Phys. **142**, 53–157 (2006). [CrossRef]

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

**5. **A. J. Devaney, “Reconstructive tomography with diffracting wave-fields,” Inv. Problems **2**, 161–183 (1986). [CrossRef]

**6. **A. Momose, T. Takeda, and Y. Itai, “Phase-contrast X-ray computed tomography for observing biological specimens and organic materials,” Rev. Sci. Instrum. **66**, 1434–1436 (1995). [CrossRef]

**7. **U. Bonse and F. Busch, “X-ray computed microtomography (µCT) using synchrotron radiation (SR),” Prog. Biophys. Mol. Biol. **65**, 133–169 (1996). [CrossRef] [PubMed]

**8. **C. Raven, A. Snigirev, I. Snigireva, P. Spanne, A. Souvorov, and V. Kohn, “Phase-contrast microtomography with coherent high-energy synchrotron x rays,” Appl. Phys. Lett. **69**, 1826–1828 (1996). [CrossRef]

**9. **P. Cloetens, W. Ludwig, J. Baruchel, D. Van Dyck, J. Van Landuyt, J. P. Guigay, and M. Schlenker, “Holotomography: Quantitative phase tomography with micrometer resolution using hard synchrotron radiation x rays,” Appl. Phys. Lett. **75**, 2912–2914 (1999). [CrossRef]

**10. **A. V. Bronnikov, “Reconstruction formulas in phase-contrast tomography,” Opt. Commun. **171**, 239–244 (1999). [CrossRef]

**11. **F. A. Dilmanian, Z. Zhong, B. Ren, X. Y. Wu, L. D. Chapman, I. Orion, and W. C. Thomlinson, “Computed tomography of x-ray index of refraction using the diffraction enhanced imaging method,” Phys. Med. Biol. **45**, 933–946 (2000). [CrossRef] [PubMed]

**12. **K. M. Pavlov, C. M. Kewish, J. R. Davis, and M. J. Morgan, “A variant on the geometrical optics approximation in diffraction enhanced tomography,” J. Phys. D: Appl. Phys. **34**, A168–A172 (2001). [CrossRef]

**13. **Y. I. Nesterets, T. E. Gureyev, and S. W. Wilkins, “General reconstruction formulas for analyzer-based computed tomography,” Appl. Phys. Lett. **89**, 264103 (2006). [CrossRef]

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

**15. **M. A. Anastasio, D. Shi, Y. Huang, and G. Gbur, “Image reconstruction in spherical-wave intensity diffraction tomography,” J. Opt. Soc. Am. A **22**, 2651–2661 (2005). [CrossRef]

**16. **X. Wu and H. Liu, “X-ray cone-beam phase tomography formulas based on phase-attenuation duality,” Opt. Express **13**, 6000–6014 (2005). [CrossRef] [PubMed]

**17. **A. Momose, W. Yashiro, Y. Takeda, Y. Suzuki, and T. Hattori, “Phase tomography by X-ray Talbot interferometry for biological imaging,” Jpn. J. Appl. Phys. **45**, 5254–5262 (2006). [CrossRef]

**18. **T. E. Gureyev, D. M. Paganin, G. R. Myers, Ya. I. Nesterets, and S. W. Wilkins, “Phase-and-amplitude computer tomography,” Appl. Phys. Lett. **89**, 034102 (2006). [CrossRef]

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

**20. **M. Engelhardt, J. Baumann, M. Schuster, C. Kottler, F. Pfeiffer, O. Bunk, and C. David, “High-resolution differential phase contrast imaging using a magnifying projection geometry with a microfocus x-ray source,” Appl. Phys. Lett. **90**, 224101 (2007). [CrossRef]

**21. **D. M. Paganin, *Coherent X-Ray Optics* (Oxford University Press, New York, 2006). [CrossRef]

**22. **G. R. Myers, S. C. Mayo, T. E. Gureyev, D. M. Paganin, and S. W. Wilkins, “Polychromatic cone-beam phase-contrast tomography,” Phys. Rev. A **76**, 045804 (2007). [CrossRef]

**23. **T. E. Gureyev, Y. I. Nesterets, K. M. Pavlov, and S. W. Wilkins, “Computed tomography with linear shift-invariant optical systems,” J. Opt. Soc. Am. A **24**, 2230–2241 (2007). [CrossRef]

**24. **J. Miao, C.-C. Chen, C. Song, Y. Nishino, Y. Kohmura, T. Ishikawa, D. Ramunno-Johnson, T.-K. Lee, and S. H. Risbud, “Three-dimensional GaN-Ga_{2}O_{3} core shell structure revealed by X-ray diffraction microscopy,” Phys. Rev. Lett. **97**, 215503 (2006). [CrossRef] [PubMed]

**25. **A. V. Bronnikov, “Phase contrast CT: Fundamental theorem and fast image reconstruction algorithms,” Proc. SPIE **6318**, 63180Q (2006). [CrossRef]

**26. **F. Pfeiffer, O. Bunk, C. Kottler, and C. David, “Tomographic reconstruction of three-dimensional objects from hard X-ray differential phase contrast projection images,” Nucl. Instrum. Methods Phys. Res. A **580**, 925–928 (2007). [CrossRef]

**27. **A. Groso, R. Abela, and M. Stampanoni, “Implementation of a fast method for high resolution phase contrast tomography,” Opt. Express **14**, 8103–8110 (2006). [CrossRef] [PubMed]

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

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

**30. **D. Paganin, S. C. Mayo, T. E. Gureyev, P. R. Miller, and S. W. Wilkins, “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” J. Microsc. **206**, 33–40 (2002). [CrossRef] [PubMed]

**31. **G. T. Herman and A. Kuba, eds., *Discrete Tomography: Foundations, Algorithms and Applications* (Birkhäuser, Boston, 1999).

**32. **F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction (SIAM, Philadelphia, 2001). [CrossRef]

**33. **R. J. Gardner and P. Gritzmann, “Discrete tomography: Determination of finite sets by x-rays,” Trans. Amer. Math. Soc. **349**, 2271–2295 (1997). [CrossRef]

**34. **P. Gritzmann, D. Prangenberg, S. de Vries, and M. Wiegelmann, “Success and failure of certain reconstruction and uniqueness algorithms in discrete tomography,” Int. J. Imaging Syst. Technol. **9**, 101–109 (1998). [CrossRef]

**35. **P. Fishburn, P. Schwander, L. Shepp, and R. J. Vanderbei, “The discrete Radon transform and its approximate inversion via linear programming,” Discrete Appl. Math. **75**, 39–61 (1997). [CrossRef]

**36. **P. Gritzmann, S. de Vries, and M. Wiegelmann, “Approximating binary images from discrete x-rays,” SIAM J. Optimization **11**, 522–546 (2000). [CrossRef]

**37. **S. Weber, T. Schüle, J. Hornegger, and C. Schnörr, “Binary tomography by iterating linear programs from noisy projections,” Lecture Notes in Comput. Sci. **3322**, 38–51 (2004). [CrossRef]

**38. **T. Schüle, C. Schnörr, S. Weber, and J. Hornegger, “Discrete tomography by convex-concave regularization and D.C. programming,” Discrete Appl. Math. **151**, 229–243 (2005). [CrossRef]

**39. **K. J. Batenburg, “A network flow algorithm for binary image reconstruction from few projections,” Lecture Notes in Comput. Sci. **4245**, 86–97 (2006). [CrossRef]

**40. **L. Hajdu and R. Tijdeman, “An algorithm for discrete tomography,” Linear Algebra Appl. **339**, 119–128 (2001). [CrossRef]

**41. **A. Alpers, H. F. Poulsen, E. Knudsen, and G. T. Herman, “A discrete tomography algorithm for improving the quality of 3DXRD grain maps,” J. Appl. Cryst. **39**, 582–588 (2006). [CrossRef]

**42. **A. Alpers and P. Gritzmann, “On the degree of ill-posedness in discrete tomography,” preprint, 2004.

**43. **A. Alpers, P. Gritzmann, and L. Thorens, “Stability and instability in discrete tomography,” Lecture Notes in Comput. Sci. **2243**, 175–186 (2001). [CrossRef]

**44. **A. Alpers and P. Gritzmann, “On stability, error correction, and noise compensation in discrete tomography,” SIAM J. Discrete Math. **20**, 227–239 (2006). [CrossRef]

**45. **R. J. Gardner, P. Gritzmann, and D. Prangenberg, “On the computational complexity of reconstructing lattice sets from their X-rays,” Discrete Math. **202**, 45–71 (1999). [CrossRef]

**46. **P. Gritzmann and S. de Vries, “Reconstructing crystalline structures from few images under high resolution transmission electron microscopy,” in Mathematics: Key Technology for the Future, W. Jäger and H.-J. Krebs, eds. (Springer-Verlag, Berlin, 2003) pp. 441–459. [CrossRef]

**47. **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] [PubMed]

**48. **T. E. Gureyev, “Transport of intensity equation for beams in an arbitrary state of temporal and spatial coherence,” Optik **110**, 263–266 (1999).

**49. **W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical recipes in C, The art of scientific computing*, 2nd ed. (Cambridge University Press, Cambridge, 1995).

**50. **D. M. Paganin, A. Barty, P. J. McMahon, and K. A. Nugent, “Quantitative phase—amplitude microscopy. III. The effects of noise,” J. Microsc. **214**, 51–61 (2004). [CrossRef] [PubMed]

**51. **D. M. Paganin, “Studies in Phase Retrieval,” PhD Thesis, University of Melbourne (1999).

**52. **E. Wolf and J. R. Shewell, “The inverse wave propagator,” Phys. Lett. A25, 417–418 (1967); see also E. Wolf and J. R. Shewell, “Errata,” Phys. Lett. A 26, 104 (1967). [CrossRef]

**53. **J. R. Shewell and E. Wolf, “Inverse diffraction and a new reciprocity theorem,” J. Opt. Soc. Am. A **58**, 1596–1603 (1968). [CrossRef]

**54. **E. Lalor, “Inverse wave propagator,” J. Math. Phys. **9**, 2001–2006 (1968). [CrossRef]

**55. **G. C. Sherman, “Diffracted wave fields expressible by plane-wave expansions containing only homogenous components,” Phys. Rev. Lett. **21**, 761–764 (1968). [CrossRef]

**56. **P. Besl, “Active, optical range imaging sensors”, Machine Vision and Applications **1**, 127–152 (1988). [CrossRef]

**57. **T. E. Gureyev, A. W. Stevenson, D. M. Paganin, T. Weitkamp, A. Snigirev, I. Snigireva, and S. W. Wilkins, “Quantitative analysis of two component samples using in-line hard X-ray images,” J. Synchrotron Rad. **9**, 148–153 (2002). [CrossRef]

**58. **W. R. Brody, G. Butt, A. Hall, and A. Macovski, “A method for selective tissue and bone visualization using dual energy scanned projection radiography,” Med. Phys. **8**, 353–357 (1980). [CrossRef]

**59. **L. A. Shepp and B. F. Logan, “Reconstructing interior head tissue from X-ray transmissions,” IEEE Trans. Nucl. Sci. **21**, 228–236 (1974). [CrossRef]