## Abstract

A method of tomographic phase retrieval is developed for multi-material objects whose components each have a distinct complex refractive index. The phase-retrieval algorithm, based on the Transport-of-Intensity equation, utilizes propagation-based X-ray phase contrast images acquired at a single defocus distance for each tomographic projection. The method requires *a priori* knowledge of the complex refractive index for each material present in the sample, together with the total projected thickness of the object at each orientation. The requirement of only a single defocus distance per projection simplifies the experimental setup and imposes no additional dose compared to conventional tomography. The algorithm was implemented using phase contrast data acquired at the SPring-8 Synchrotron facility in Japan. The three-dimensional (3D) complex refractive index distribution of a multi-material test object was quantitatively reconstructed using a single X-ray phase-contrast image per projection. The technique is robust in the presence of noise, compared to conventional absorption based tomography.

©2010 Optical Society of America

## 1. Introduction

X-ray absorption is a powerful tool for imaging objects comprised of multiple materials. In this absorptive imaging regime, the logarithm of the intensity image formed over the detector is proportional to the projection of the object’s X-ray linear attenuation coefficient along a series of straight-line raypaths [1]. However, features that exhibit poor absorption contrast, such as biological soft tissues, can be difficult to resolve with this method.

Numerous techniques have been developed to successfully overcome this difficulty, including K-edge subtraction imaging [2, 3], magnetic resonance imaging [4, 5], positron emission tomography [5], X-ray interferometry [6], analyzer-based phase contrast [7], and X-ray grating methods [8].

The technique utilized in the present paper, which differs from the abovementioned techniques insofar as it employs no optical elements between the illuminating radiation and the detector, is propagation-based X-ray phase contrast imaging (PBI) [9–11]. In this imaging modality, the act of free-space propagation—from the exit surface of a sample illuminated by a spatially coherent source, to the surface of a two-dimensional position-sensitive detector—renders visible the transverse phase shifts imparted by the sample upon the illuminating radiation (see Fig. 1 ).

Although PBI can improve the visualization of weakly absorbing features in a sample, quantitative information cannot be directly obtained from the raw phase contrast X-ray images. Phase-retrieval methods for extracting quantitative information from intensity measurements alone have been developed (see e.g. [12–21]) to obtain the projected phase and absorption information of the object. Notwithstanding their successes, these methods often require multiple intensity measurements, impose strong restrictions on the object under study or apply iterative solution techniques. Acquiring multiple images can prove problematic for correct alignment of images and induces a higher radiation dose, which is important for biomedical applications. Several phase retrieval methods require the object to be ‘weak’ such that they provide little to no absorption contrast with limited phase gradients introduced by the sample [18, 20, 21]. Since most inanimate materials and biological tissues cannot be considered as ‘weak objects’, phase retrieval algorithms developed under these approximations have only limited use for biomedical imaging or materials science applications. Iterative phase retrieval algorithms can also be problematic as convergence to the correct solution cannot be guaranteed and are computationally more intensive than analytic solutions.

For the case of a single-material object illuminated by paraxial coherent X-rays, Paganin *et al.* [22] developed a noise-robust deterministic phase-retrieval algorithm to reconstruct the projected linear attenuation coefficient from a single PBI image, which has been subsequently utilized in a number of tomographic studies [23–25]. These phase-contrast tomography investigations, which incorporate the effects of both absorption and phase contrast, generalize the seminal work for the phase-contrast tomography of pure phase objects by Bronnikov [26,27].

In this paper, we extend the work of Paganin *et al*. [22] and Mayo *et al*. [23] to enable analytic propagation-based phase-retrieval tomography to be performed on a multi-material object in which the spatially-dependent complex refractive index is quantized, *i.e.* it takes one of a series of distinct values. The algorithm makes use of a single PBI image per projection, separately and selectively reconstructing each interface between any given pair of distinct materials. Having separately reconstructed the interfaces between different distinct pairs of materials, a “spliced” image of all materials present in the sample can then be computed. Section 2 of the paper develops the underlying theory for our method, with an experimental implementation in Section. 3.

## 2. Theory

Here we outline the theory for quantitative phase–amplitude tomography of a multi-material object with a quantized distribution of complex refractive indices (“spatially-quantized object”). Section 2.1 reviews an existing theory for two-dimensional (2D) single-image propagation-based phase-retrieval of a one-material object, with Section 2.2 generalizing this to the case of 2D single-image phase retrieval of an object of a given material which is embedded in a matrix of a second material. Section 2.3 generalizes this result to the case of three-dimensional (3D) propagation-based phase-contrast imaging of spatially-quantized objects, given a single image per projection.

#### 2.1 Two-dimensional phase retrieval for a binary object

Consider an object with complex refractive index distribution $n(\mathbf{p})=1-\delta (\mathbf{p})+i\beta (\mathbf{p})$, where *δ* and $\beta =\lambda \mu /4\pi $ respectively quantify the refractive and absorptive properties of the object, *λ* is the radiation wavelength and *μ* is the linear attenuation coefficient [28]. The position vector $\mathbf{p}=({p}_{1},{p}_{2},{p}_{3})$ represents the object’s coordinate system as illustrated in Fig. 1. For an object composed of a single material “1” (including voids; herein referred to as a binary object), $n(\mathbf{p})$ takes only the values $1-{\delta}_{1}+i{\beta}_{1}$ or unity. For such objects, and under the assumptions of paraxial monochromatic scalar electromagnetic-wave illumination and optical thinness, the intensity and phase of the wavefield downstream of the object obey the Transport-of-Intensity equation (TIE) [12]:

Here, the intensity and phase of the wavefield are denoted by $I({\mathbf{r}}_{\perp},z)$ and $\phi ({\mathbf{r}}_{\perp},z)$, ${\mathbf{r}}_{\perp}=(x,y)$ is the position vector perpendicular to the optic axis *z*, ${\nabla}_{\perp}=({\partial}_{x},{\partial}_{y})$ is the transverse gradient operator, $\mathbf{k}=2\pi /\lambda $ is the wave number and ${\partial}_{m}$ denotes differentiation with respect to $m=x,y,z$. For normally incident plane-wave illumination of an optically-thin object, the intensity and phase at the contact plane $(z=0)$ is given by the projection approximation (see e.g. [29]), such that $I({\mathbf{r}}_{\perp},z=0)={I}_{0}\mathrm{exp}[-{\mu}_{1}{T}_{1}({\mathbf{r}}_{\perp})]$ and $\phi ({\mathbf{r}}_{\perp},z=0)=-k{\delta}_{1}{T}_{1}({\mathbf{r}}_{\perp})$. The former expression is the Beer–Lambert law of absorption, ${T}_{1}({\mathbf{r}}_{\perp})$ is the projected thickness of the object (projected in the *z* direction) and *I*_{0} is the uniform incident intensity. By substituting $I({\mathbf{r}}_{\perp},z=0)$and $\phi ({\mathbf{r}}_{\perp},z=0)$ into Eq. (1) one can solve for the projected thickness as [22]:

Here, **F** and ${\mathbf{F}}_{}^{-1}$are the two-dimensional forward and inverse Fourier transforms with respect to ${\mathbf{r}}_{\perp}$, ${\mathbf{k}}_{\perp}=({k}_{x},{k}_{y})$ are the Fourier coordinates corresponding to ${\mathbf{r}}_{\perp}$, and *d* is the propagation distance between the exit-surface of the object and the surface of the detector (see Fig. 1). This algorithm assumes the “near field approximation” that *d* be sufficiently small that the Fresnel number ${N}_{F}={\ell}^{2}/[\lambda (D+d)]$ be significantly larger than unity, where *ℓ* is the characteristic transverse length scale associated with the wavefield at the exit surface of the object and *D* is the distance from the smallest feature to the exit surface (see e.g. [29]). For a given wavelength and spatial resolution we can therefore estimate the minimum feature size that can be imaged and still satisfy the near field approximation.

A mathematically identical form of Eq. (2) that allows for the recovery of the projected electron density $\rho ({\mathbf{r}}_{\perp})$, has also been derived via a different approach by Wu *et al.* [30]. In their derivation, the authors make explicit use of the Klein-Nishina formula and assume the energies of the incident X-rays must range from 60 to 500 keV. Their derivation is not restricted to a single material object, but enables the same mathematical simplification as the single material assumption made in deriving Eq. (2) when lower energy X-rays are considered.

#### 2.2 Two-dimensional phase retrieval for a ternary object

Here, the phase-retrieval algorithm in Eq. (2) is generalized to the case of a ternary object (two materials including voids). To this end, consider an object made of material “j” with projected thickness ${T}_{j}({\mathbf{r}}_{\perp})$, embedded in another object made of material “1” with projected thickness ${T}_{1}({\mathbf{r}}_{\perp})$. We seek to recover ${T}_{j}({\mathbf{r}}_{\perp})$ from a single PCI image. Introduce the total projected thickness [25]:

Under the projection approximation the intensity and phase at the contact plane are given by

We can solve for ${T}_{j}({\mathbf{r}}_{\perp})$ upon substitution of Eqs. (3), (4) and (5) into the TIE (Eq. (1)) and making the assumption that the projected thickness of the encasing material, ${T}_{1}({\mathbf{r}}_{\perp})$, varies slowly across the plane ${\mathbf{r}}_{\perp}=(x,y)$. Under this approximation, spatial derivatives of the encasing material are neglected. As such, the term ${\delta}_{1}{I}_{0}\mathrm{exp}[-{\mu}_{1}A({\mathbf{r}}_{\perp})]{\nabla}_{\perp}\cdot \{\mathrm{exp}[-({\mu}_{j}-{\mu}_{1}){T}_{j}({\mathbf{r}}_{\perp})]{\nabla}_{\perp}A({\mathbf{r}}_{\perp})\}$ appearing in the derivation is ignored, which leads to the following non-linear differential equation:

Fourier transform both sides of Eq. (8), use the Fourier derivative theorem, and inverse transform to obtain the following expression for ${T}_{j}({\mathbf{r}}_{\perp})$:

To employ this algorithm one requires *a priori* knowledge of the total projected thickness $A({\mathbf{r}}_{\perp})$ in addition to the values *δ*
_{j}, *δ*
_{1} and *μ*
_{j}, *μ*
_{1} corresponding to the different materials in the object. If an object is of complicated shape and is known to contain no internal voids, then $A({\mathbf{r}}_{\perp})$ can be found using techniques such as laser profilometry [25]. In the case of tomography, a more convenient and practical approach to obtain $A({\mathbf{r}}_{\perp})$ for each projection angle can be used. By utilizing Eq. (2) the encasing material can first be correctly reconstructed tomographically. Internal and external voids can then be located by computationally searching for a predefined threshold in each slice of the reconstructed volume. $A({\mathbf{r}}_{\perp})$ may then be calculated for each projection angle. On an additional note, for point source illumination it is necessary to account for image magnification in both Eq. (2) and Eq. (9). Both equations will be altered by identical factors [31].

#### 2.3 Three-dimensional phase retrieval: Quantitative phase–amplitude tomography of multi-material objects with quantized refractive-index distribution, from a single phase-contrast image per projection

For both absorption and phase contrast tomography we utilize a conventional filtered backprojection (FBP) algorithm to reconstruct a two-dimensional slice of the imaged object [1, 24]. Fourier-transform-based FBP utilizing projection phase-contrast images processed using Eq. (2), allows one to “focus in” on an interface between air (cavities) and the encasing object “1”, enabling any voids within the encasing material to be quantitatively reconstructed. Here *δ*
_{1} and *μ*
_{1} take the values corresponding to the encasing object “1”. For objects with spatially quantized refractive indices, such as shown in Fig. 2
, the “correct” interfaces will be sharply reconstructed.

Interfaces between embedded objects “j” (where $j=2,3,\mathrm{...}$) and the encasing object “1”, will be incorrectly reconstructed using this procedure. However, such incorrectly reconstructed interfaces will bear a characteristic local signature of either (i) residual phase-contrast fringes, or (ii) blurred interfaces. Moreover, such “incorrect” interfaces (e.g. between materials “j” and “1”) may be quantitatively (sharply) reconstructed by repeating the FBP analysis using inputs derived by inserting ${(\Delta \delta )}_{j1}\equiv {\delta}_{j}-{\delta}_{1}$ and ${(\Delta \mu )}_{j1}\equiv {\mu}_{j}-{\mu}_{1}$ into Eq. (9), provided the material of interest is not in contact with or in the immediate vicinity of either another embedded object or a cavity. From a series of FBP analyses for each distinct pair of material interfaces present in the sample, one can build up a quantitative 3D map of all materials. These can be spliced together into a composite tomographic reconstruction, in which all interfaces between each distinct pair of materials is sharply reconstructed. This is the core idea underpinning our procedure for phase–amplitude tomography of a multi-material object with spatially quantized complex refractive index distribution, given a single propagation-based phase contrast image per projection.

Crucial to the success of the splicing procedure described above is the idea that “incorrect” interfaces are only *locally* polluted (in three dimensions) by locally incorrect choices for *δ* and *μ* [23]. How may one quantify the smoothing artefacts for interfaces which are incurred as a result of a locally incorrect filtering parameter, $\alpha =d\delta /\mu $, being utilized in Eq. (2) at a given interface? In regions where over-smoothing occurs (*i.e.* where *α* is underestimated), the amount of blurring (“bleed width”) of an “incorrect” interface can be estimated by analyzing the low-pass Fourier filter in Eq. (2), which has the form of a Lorentzian function. This Lorentzian can be rewritten as

In the case of a medium “j” embedded in an encasing medium “1”, $\epsilon =d({\delta}_{j}-{\delta}_{1})/({\mu}_{j}-{\mu}_{1})$ is the ratio that will correctly reconstruct this particular interface. Equation (10) allows us to approximately interpret the operation of Eq. (2) (using *α* rather than *ε*) on this interface as the correct Lorentzian function multiplied by the ratio of quadratic functions on the right side of Eq. (10). This quotient is a Fourier filter that locally causes either an over- or under-smoothing artefact in the vicinity of the “incorrect” interface. Incidentally, in the case of over-smoothing, when $\epsilon /\alpha <1$, this ratio of quadratic functions is a Lorentzian function that has been vertically scaled by a factor $(\alpha -\epsilon )/\epsilon $ and shifted positively along the vertical axis by $(\epsilon /\alpha )$, with a maximum value of unity. In the opposite case, when $\epsilon /\alpha >1$, we have an inverted Lorentzian, giving rise to a high-pass filter that leaves a characteristic fringe representing under-compensated phase contrast. In either case, one has a clear signature associated with “incorrect” interfaces, which allows them to be unambiguously identified and thereby enabling one to “focus in” upon a particular material interface of interest.

The “bleed width” for over-smoothing can be estimated via the optical uncertainty principle:

where Δ**x**is the uncertainty in real space (

*i.e.*, the “bleed width”) and Δ

**k**is the uncertainty in reciprocal space. Take Δ

**k**as the half-width at half-maximum (HWHM) of the second quotient on the right-hand side of Eq. (10); this HWHM is taken relative to the asymptotic baseline of $\epsilon /\alpha $. Subsequently, using Eq. (11) we obtain a simple lower bound for the associated bleed width that is independent of

*ε*:

We close this section by noting that this “bleed width”, associated with interfaces *other* than those selected via a given choice for *α*, is constructive insofar as it “tags” particular over-smoothed interfaces as not corresponding to the particular interface of interest, but rather to an interface with $\epsilon <\alpha $. Under-smoothed interfaces are also tagged, as corresponding to $\epsilon >\alpha $, via the signature of incompletely-compensated phase contrast fringes.

## 3. Experimental results and discussion

Here we give an X-ray synchrotron-based experimental implementation of the theory developed above. Section 3.1 discusses the experimental setup. Section 3.2 reports 2D phase–amplitude retrieval of multi-material spatially quantized objects with embedded features that are non-overlapping in projection (cf. Sections 2.1 and 2.2). This requirement for non-overlap in projection is relaxed is Section 3.3, which implements 3D phase–amplitude retrieval of multi-material spatially quantized objects, given a single propagation-based X-ray phase contrast image per projection (cf. Section 2.3).

#### 3.1 Experimental setup for X-ray phase–amplitude imaging of spatially quantized objects

X-ray propagation-based phase–amplitude computed tomography experiments were performed in Hutch 3 of beamline 20B2 at the SPring-8 synchrotron radiation source, Japan (proposal 2009A1882). The beamline employs a bending magnet source with a Si(111) double-crystal monochromator [32]. An X-ray energy of 24 keV provided acceptable phase and attenuation contrast for materials of relatively low atomic number. This energy is also commensurate with the diagnostic X-ray energies used in mammography [33]. Phase and absorption contrast images were recorded with the detector positioned at distances of $d=1\text{\hspace{0.17em}}m$ and $d<3\text{\hspace{0.17em}}cm$from the object, respectively, as shown in Fig. 1. An optically coupled Hamamatsu CCD camera (C4880-41S) with a 10 μm thick gadolinium oxysulfide (Gd_{2}O_{2}S) scintillator was used to acquire each data set. The images had a window size of 3000 x 1500 pixels and an effective pixel size of 5.9 *μ*m, which gave a region of interest (ROI) of 17.70 mm (H) × 8.85 mm (V). For each tomographic data set, 1800 projections were collected over 180° of rotation, each acquired with an exposure time of 2.5 s. Flat field images (no object in the beam) were recorded every 43 projections. Frequent flat field recordings were required due to instabilities in beam position with time. Dark field images also were acquired with the main beam shutter closed at the beginning and end of each scan to correct for the detector dark current offset. Prior to quantitative analysis, all phase-contrast images were flat-field and dark-field corrected, to compensate for both illumination non-uniformities and CCD dark current.

The imaged object was a PMMA (Polymethyl-methacrylate (C_{5}H_{8}O_{2}); commonly known as Perspex) cylinder 10 mm in height and 12.75 ± 0.05 mm in diameter. This contained four cavities each with diameter 1.02 ± 0.05 mm. An Aluminium (Al) and PTFE (Polytetrafluoroethylene (H_{2}F_{4}); commonly known as Teflon) pin of 1.00 ± 0.05 mm diameter were inserted into two of the cavities to create a quaternary object (three materials plus voids) with three distinct interfaces, these being air/PMMA, Aluminium/PMMA and PTFE/PMMA (see Fig. 2). The *δ* and *μ* values for the materials in the object, for 24 keV X-rays, are listed in Table 1
. The Aluminium/PMMA interface will be denoted as $j=2$ and the PTFE/PMMA interface as $j=3$.

#### 3.2 Projection imaging

Here we use the theory of Section 2.2 to implement 2D phase–amplitude retrieval of a multi-material spatially quantized object with non-overlapping projections of embedded features composed of different materials, noting that the assumption of non-overlapping embedded features will be dropped when moving to a tomographic analysis in Section. 3.3.

Absorption and phase contrast images of our test object (see Fig. 2) are shown in Fig. 3(a) and (b) , for an orientation in which all embedded features (corresponding to different materials) are non-overlapping in projection. In the absorption contrast image, all interfaces except Al/PMMA exhibit weak contrast, since Al is significantly more attenuating than the other materials in the sample (see Table 1). Conversely, the phase contrast image displays a significantly greater signal-to-noise ratio (SNR), dramatically enhancing the PTFE/PMMA and air/PMMA interfaces and improving visualization of the object overall [10]. Line profiles from left-to-right of the centre of the images in Fig. 3(a) and (b) are shown in Fig. 3(f) and (g), respectively. The profile in Fig. 3(g) clearly reveals strong propagation-based phase contrast fringes formed at interface boundaries [10].

Figure 3(c) shows an image of the projected thickness of the object computed via Eq. (2), under the assumption that the object was composed entirely of PMMA. Here we have used the known values of *δ*
_{1} and *μ*
_{1} for PMMA listed in Table 1. The line profile in Fig. 3(h) shows the distribution of the recovered projected thickness from the phase contrast image in Fig. 3(b). From the profile it can be seen that the maximum projected thickness of the PMMA is close to the expected values (12.75 mm) and the boundaries of the air/PMMA interfaces are sharp. Moreover, due to reasons mentioned in section 2.3, the Al/PMMA and PTFE/PMMA interfaces suffer from an “over smoothing” and the associated projected thicknesses are also overestimated.

To recover the projected thicknesses of the Al and PTFE, that is ${T}_{2}({\mathbf{r}}_{\perp})$ and ${T}_{3}({\mathbf{r}}_{\perp})$, we “focus” on each of these individually using Eq. (9). Since implementation of this algorithm requires knowledge of the total projected thickness, the function $A({\mathbf{r}}_{\perp})$ was generated with the assumption that the object was a cylinder with no internal voids, which results in negative values arising in the projected thickness in the presence of internal voids (see Fig. 3(i) and (j)). To recover ${T}_{2}({\mathbf{r}}_{\perp})$, the values used for *δ*
_{2} and *μ*
_{2} correspond to Aluminium and *δ*
_{1} and *μ*
_{1} to PMMA (see Table 1). To recover ${T}_{3}({\mathbf{r}}_{\perp})$, the values used for *δ*
_{3} and *μ*
_{3} correspond to PTFE. The phase-retrieved images are shown in Fig. 3(d) and (e). The images show how the “selected” interfaces are sharpened and how the attenuated intensity due to PMMA is effectively removed. Line profiles in Fig. 3(i) and (j) illustrate this more clearly and reveal that the recovered projected thicknesses at the “selected” interfaces ${T}_{2}({\mathbf{r}}_{\perp})$ and ${T}_{3}({\mathbf{r}}_{\perp})$ are close to their expected thickness (~1 mm). Additionally, we see that the algorithm breaks down in regions where the encasing material exhibits strong phase contrast; however, this occurs at the boundary of the encasing material far from the interfaces of interest. These artefacts are expected since Eq. (9) assumes that phase gradients due to the encasing material are negligible. Also, we now see that the interfaces where the density is less than that of Aluminium are now under-smoothed and slight residual fringing is still visible (see for example Fig. 3(i), interfaces 1, 3 and 4; cf. Section. 2.3). Note that the upper threshold in (e) has been reduced to enhance the PTFE, as a linear palette obscures this feature and the Al/PMMA interface appears saturated as a consequence.

#### 3.3 Tomography of spatially quantized objects, from a single phase-contrast image per projection

Section 3.3.1 presents our results for interface-specific phase-retrieval tomography, whereby a single phase-contrast image per projection may be used to selectively reconstruct a given interface between two given materials, in a spatially-quantized object. In Section 3.3.2, we experimentally test our uncertainty-principle estimate (see Eqs. (11) and (12)) for the “bleed width” associated with material interfaces other than that which is selectively reconstructed. In Section 3.3.3 we show how a complete set of material interfaces in a given spatially-quantized object may be spliced together, into a composite tomogram, given a single phase-contrast view per projection.

### 3.3.1 Interface-specific phase-retrieval tomography

For absorption-contrast FBP tomography, using the parameters described in Section 3.1, the reconstruction of a single slice of the sample is shown in Fig. 4(a) . Qualitatively the reconstruction shows good detail of the object features, particularly at the Al/PMMA interface where there is a large difference in absorption. In Fig. 4(f) line profiles from left-to-right of the centre of each interface of Fig. 4(a) are shown in a single plot. This displays the distribution of the attenuation coefficient at each interface where it can be seen that, on average, the attenuation coefficients are similar to the theoretical values in Table 1. However, the reconstruction also contains a substantial amount of high frequency noise making it difficult to define the boundaries at the interfaces where absorption contrast is poor (air/PMMA and PTFE/PMMA). Such noise is typically reduced in tomography by use of an additional low-pass filter during the filtered backprojection, but must be traded against a reduction in spatial resolution [1].

Figure 4(b) shows a slice of a tomographic reconstruction obtained by applying FBP directly to the raw phase contrast images (cf. [34]), without any phase–amplitude retrieval. Like the absorption contrast result it also reveals good detail of the object features. The reconstruction allows us to see how the boundaries of each interface are enhanced as a result of the high SNR provided by phase contrast. Although phase contrast provides better feature visibility, the line profile in Fig. 4(g) shows that non-physical negative values and sharp spikes arise in the attenuation coefficient map without the phase retrieval step.

Figure 4(c) shows a slice of a tomographic reconstruction obtained using Eq. (2) to yield ${T}_{1}({\mathbf{r}}_{\perp},\theta )$ for each projection angle, giving a series of two-dimensional projected-thickness maps that were then tomographically reconstructed using FBP [24]. The reconstruction shows the distribution of ${\delta}_{1}({\mathbf{p}}_{\perp})$ along the plane ${\mathbf{p}}_{\perp}=({p}_{1},{p}_{3})$ (see Fig. 1), which in this case is PMMA $(j=1)$. It can be seen that the noise has been substantially suppressed while preserving the sharpness of the air/PMMA interface, as the line profile in Fig. 4(h) indicates. Conversely, over-smoothing (cf. Section 2.3) is clearly apparent at the remaining interfaces. Despite this localized blurring of these regions the result shows that Eq. (2) can be applied to multi-material objects and accurately reconstruct the voids as long as they are not near the vicinity of the over-smoothed interfaces. The line profile in Fig. 4(h) shows that at the air/PMMA interface, the distribution of ${\delta}_{1}({\mathbf{p}}_{\perp})$ has an average value of 4.2 × 10^{−7}, which is within 10% of the theoretical value in Table 1. This is acceptable as it is known that attenuation coefficients have discrepancies of up to 10% between theoretical values [35].

Figure 4(d) and (e) show reconstructions of the same slice from phase retrieved images using Eq. (9) to obtain ${T}_{2}({\mathbf{r}}_{\perp},\theta )$ and ${T}_{3}({\mathbf{r}}_{\perp},\theta )$ before using FBP to respectively recover the distributions of ${\delta}_{2}({\mathbf{p}}_{\perp})$ (Aluminium; $j=2$) and ${\delta}_{3}({\mathbf{p}}_{\perp})$ (PTFE; $j=3$). Figure 4(d) focuses on the Al/PMMA interface and (e) focuses on the PTFE/PMMA interface. Due to the rotational symmetry of the surface of the cylinder the same $A({\mathbf{r}}_{\perp})$ fitted function was used for

every projection. The distribution of ${\delta}_{2}({\mathbf{p}}_{\perp})$and ${\delta}_{3}({\mathbf{p}}_{\perp})$are shown in the line profiles in Fig. 4(i) and (j). On average the values of ${\delta}_{2}({\mathbf{p}}_{\perp})$ and ${\delta}_{3}({\mathbf{p}}_{\perp})$ are 7.9 × 10^{−7} and 6.6 × 10^{−7} respectively, and are within 15% of the theoretical values in Table 1.

### 3.3.2 Quantification of over-smoothing

Experimentally, we find that the “bleed width” blurring typically extends between 3 and 5 times the lower bound in Eq. (12). Hence the following rule-of-thumb is used to estimate the bleed width $\Delta \mathbf{x}$:

This formula gives a rough indication as to how far an over-smoothed interface will be blurred by the phase retrieval (cf. Section 2.3). We tested this rule by measuring the blurring at the Al/PMMA and PTFE/PMMA interfaces in Fig. 4 (h). The blurring was just visible up to a distance of 443 *μ*m from the Al/PMMA interface and 454 *μ*m from the PTFE/PMMA interface with an uncertainty of ±30 *μ*m. Evaluating *α* for PMMA gives a value of 106 *μ*m. The measured bleed widths are thus 4.17×*α* for Al/PMMA and 4.28× *α* for PTFE/PMMA, which agrees with the inequality in Eq. (13). We also note that both interfaces bleed almost the same amount, providing empirical evidence that the parameter *α* is the only contributor to the smoothing artefact.

### 3.3.3 Spliced tomographic reconstruction

Having quantitatively reconstructed all three interfaces and established a rule-of-thumb to estimate the blurring widths (see Eq. (13)), we can now combine these images to form a single “spliced” image of our quaternary (or higher order) object (cf. Section. 2.3). The image was constructed by digitally inserting the individually reconstructed interfaces in Fig. 4(d) and (e) into the appropriate regions of the encasing material in Fig. 4(c). The size of the region was chosen by considering the amount of blurring at these interfaces. The refractive index of these segments were offset so that the background zeroes match that of the encasing material and the resulting amplitude was rescaled to maintain the original value. The spliced image is shown in Fig. 5(a) . Line profiles of the different interfaces of the image are shown in Fig. 5(b). We now have a quantitative tomographic reconstruction of a multi-material object, at each point of which the refractive index takes one of N distinct values, using a single PCI image per projection.

To quantitatively compare our spliced reconstruction with the absorption-based reconstruction, the SNR was calculated for each medium. The Al, PTFE and PMMA gave SNRs of 312, 309 and 98.9, respectively, where the signal is the average value inside the medium. This significantly improves the respective SNR values of 18.7, 5.19 and 1.17, which were calculated at the same regions in the absorption contrast reconstruction in Fig. 4(a). The ordering of these SNRs is consistent with the relative amount of smoothing associated with the filters in Eqs. (2) and (9). The increase in SNR is largest (~85×) for the encasing material as the Fourier-space damping of high spatial frequencies is largest for the PMMA/air interface. Each of the SNRs were calculated over a 90×90 pixel region of interest. The mathematical origin of the stability with respect to noise is due to the regularizing presence of a non-zero denominator in the Fourier filter of Eqs. (2) and (9), as the ratio $d\delta /\mu $ is always greater than zero and the ratio $d({\delta}_{\text{j}}-{\delta}_{1})/({\mu}_{\text{j}}-{\mu}_{1})$ is unlikely to be less than zero; also, since *μ*is never zero and ${(\Delta \mu )}_{\text{j}1}$ is generally a non-zero quantity. This avoids instability problems of the ‘division by zero’ type that would arise if the denominator were zero [22].

## 4. Conclusion

We have developed a numerically efficient phase retrieval algorithm to reconstruct the complex refractive index distribution of known materials embedded in a second medium from a single X-ray phase contrast image per projection. This interface-specific phase-retrieval tomography algorithm requires *a priori* knowledge of each material’s complex refractive index and the total projected thickness of the sample at each orientation. The algorithm was successfully applied to experimental data collected using X-ray synchrotron radiation. We have used our method to quantitatively reconstruct a quaternary object by reconstructing all interfaces separately. A complete tomographic reconstruction of all interfaces in the object was produced by splicing the individual reconstructions together. When it comes to imaging multi-material objects, our technique proves superior to conventional absorption contrast in terms of the signal-to-noise ratio (SNR). For the test sample used here, the improvement in the SNR was between 17 and 85 fold.

Given the success of our results, combined with the practicality of the technique, we anticipate that it will have applications in a number of different fields including medicine, biomedical science, geology and materials science. Our methodology may also be adapted to electron and visible light microscopy.

## Acknowledgements

MAB acknowledges funding from the Monash University Dean’s Scholarship Scheme. MJK and DMP acknowledge funding from the Australian Research Council, via the Discovery-Projects Programme. DMP acknowledges Timur Gureyev (Commonwealth Scientific and Industrial Research Organization, Clayton, Australia), Glenn Myers (Australian National University, Canberra) and Peter Cloetens (European Synchrotron Radiation Facility, Grenoble, France), for stimulating discussions.

## References and links

**1. **A. C. Kak, and M. Slaney, *Principles of Computerized Tomographic Imaging* (IEEE Press, New York, 1988).

**2. **S. Bayat, G. Le Duc, L. Porra, G. Berruyer, C. Nemoz, S. Monfraix, S. Fiedler, W. Thomlinson, P. Suortti, C. G. Standertskjöld-Nordenstam, and A. R. A. Sovijärvi, “Quantitative functional lung imaging with synchrotron radiation using inhaled xenon as contrast agent,” Phys. Med. Biol. **46**(12), 3287–3299 (2001). [CrossRef]

**3. **S. Monfraix, S. Bayat, L. Porra, G. Berruyer, C. Nemoz, W. Thomlinson, P. Suortti, and A. R. A. Sovijärvi, “Quantitative measurement of regional lung gas volume by synchrotron radiation computed tomography,” Phys. Med. Biol. **50**(1), 1–11 (2005). [CrossRef] [PubMed]

**4. **B. Driehuys and L. W. Hedlund, “Imaging techniques for small animal models of pulmonary disease: MR microscopy,” Toxicol. Pathol. **35**(1), 49–58 (2007). [CrossRef] [PubMed]

**5. **D. P. Schuster, A. Kovacs, J. Garbow, and D. Piwnica-Worms, “Recent advances in imaging the lungs of intact small animals,” Am. J. Respir. Cell Mol. Biol. **30**(2), 129–138 (2004). [CrossRef] [PubMed]

**6. **U. Bonse and M. Hart, “An x-ray interferometer,” Appl. Phys. Lett. **6**(8), 155–156 (1965). [CrossRef]

**7. **E. Förster, K. Goetz, and P. Zaumseil, “Double crystal diffractometry for the characterization of targets for laser fusion experiments,” Krist. Tech. **15**(8), 937–945 (1980). [CrossRef]

**8. **F. Pfeiffer, T. Weitkamp, O. Bunk, and O. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nature **2**, 258–261 (2006).

**9. **A. Snigirev, I. Snigireva, V. Kohn, S. Kuznetsov, and I. Schelokov, “On the possibilities of x-ray phase contrast microimaging by coherent high-energy synchrotron radiation,” Rev. Sci. Instrum. **66**(12), 5486–5492 (1995). [CrossRef]

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

**11. **P. Cloetens, R. Barrett, J. Baruchel, J.-P. Guigay, and M. Schlenker, “Phase objects in synchrotron radiation hard X-ray imaging,” J. Phys. D Appl. Phys. **29**(1), 133–146 (1996). [CrossRef]

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

**13. **T. E. Gureyev and K. A. Nugent, “Phase retrieval with the transport-of-intensity equation. II. Orthogonal series solution for nonuniform illumination,” J. Opt. Soc. Am. A **13**(8), 1670–1682 (1996). [CrossRef]

**14. **T. E. Gureyev and K. A. Nugent, “Rapid quantitative phase imaging using the transport of intensity equation,” Opt. Commun. **133**(1-6), 339–346 (1997). [CrossRef]

**15. **D. Paganin and K. A. Nugent, “Noninterferometric phase imaging with partially coherent light,” Phys. Rev. Lett. **80**(12), 2586–2589 (1998). [CrossRef]

**16. **O. M. Bucci, L. Crocco, M. D'Urso, and T. Isernia, “Inverse scattering from phaseless measurements of the total field open lines,” J. Opt. Soc. Am. A **23**(10), 2566–2577 (2006). [CrossRef]

**17. **L. Crocco, M. D’Urso, and T. Isernia, “Faithful non-linear imaging from only-amplitude measurements of incident and total fields,” Opt. Express **15**(7), 3804–3815 (2007). [CrossRef] [PubMed]

**18. **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**(12), 1617–1619 (2007). [CrossRef] [PubMed]

**19. **M. D'Urso, K. Belkerbir, L. Crocco, T. Isernia, and A. Litman, “Phaseless imaging with experimantal data: facts and challenges,” J. Opt. Soc. Am. A **25**(1), 271–281 (2008). [CrossRef]

**20. **A. Pogany, D. Gao, and S. W. Wilkins, “Contrast and resolution in imaging with a microfocus x-ray source,” Rev. Sci. Instrum. **68**(7), 2774–2782 (1997). [CrossRef]

**21. **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**(19), 2912–2914 (1999). [CrossRef]

**22. **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**(1), 33–40 (2002). [CrossRef] [PubMed]

**23. **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**(19), 2289–2302 (2003). [CrossRef] [PubMed]

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

**25. **G. R. Myers, T. E. Gureyev, D. M. Paganin, and S. C. Mayo, “The binary dissector: phase contrast tomography of two- and three-material objects from few projections,” Opt. Express **16**(14), 10736–10749 (2008). [CrossRef] [PubMed]

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

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

**28. **J. Als-Nielsen, and D. McMorrow, *Elements of Modern X-ray Physics* (John Wiley & Sons, New York, 2001).

**29. **D. M. Paganin, *Coherent X-Ray Optics* (Oxford University Press, New York, 2006).

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

**31. **M. J. Kitchen, R. A. Lewis, M. J. Morgan, M. J. Wallace, M. L. Siew, K. K. W. Siu, A. Habib, A. Fouras, N. Yagi, K. Uesugi, and S. B. Hooper, “Dynamic measures of regional lung air volume using phase contrast x-ray imaging,” Phys. Med. Biol. **53**(21), 6065–6077 (2008). [CrossRef] [PubMed]

**32. **S. Goto, K. Takeshita, Y. Suzuki, H. Ohashi, Y. Asano, H. Kimura, T. Matsushita, N. Yagi, M. Isshiki, H. Yamazaki, Y. Yoneda, K. Umetani, and T. Ishikawa, “Construction and commissioning of a 215m-long beamline at SPring-8,” Nucl. Instrum. Meth. A. **467-468**, 682–685 (2001). [CrossRef]

**33. **A. Abrami, F. Arfelli, R. C. Barroso, A. Bergamaschi, F. Bille, P. Bregant, F. Brizzi, K. Casarin, E. Castelli, V. Chenda, L. Dalla Palma, D. Dreossi, C. Fava, R. Longo, L. Mancini, R.-H. Menk, F. Montanari, A. Olivo, S. Pani, A. Pillon, E. Quai, S. R. Kaiser, L. Rigon, T. Rokvic, M. Tonutti, G. Tromba, A. Vascotto, C. Venanzi, F. Zanconati, A. Zanetti, and F. Zanini, “Medical applications of synchrotron radiation at the SYRMEP beamline of ELETTRA,” Nucl. Instrum. Meth. A **548**(1-2), 221–227 (2005). [CrossRef]

**34. **P. Cloetens, W. Ludwig, J. Baruchel, J.-P. Guigay, P. Pernot-Rejmánková, M. Salomé-Pateyron, M. Schlenker, J.-Y. Buffière, E. Malre, and G. Peix, “Hard X-ray phase imaging using simple propagation of a coherent synchrotron radiation beam,” J. Phys. D Appl. Phys. **32**(10A), 145–151 (1999). [CrossRef]

**35. **C. T. Chantler, C. Q. Tran, Z. Barnea, D. Paterson, D. J. Cookson, and D. X. Balaic, “Measurement of the x-ray mass attenuation coefficient of copper using 8.85-20 keV synchrotron radiation,” Phys. Rev. A **64**(6), 062506 (2001). [CrossRef]