Fast three-dimensional phase retrieval in propagation-based X-ray tomography

The following article describes a method for 3D reconstruction of multi-material objects based on propagation-based X-ray phase-contrast tomography (PB-CT) with phase retrieval using the homogeneous form of the transport of intensity equation (TIE-Hom). Unlike conventional PB-CT algorithms that perform phase retrieval of individual projections, the described post-reconstruction phase-retrieval method is applied in 3D to a localized region of the CT-reconstructed volume. This work demonstrates, via numerical simulations, the accuracy and noise characteristics of the method under a variety of experimental conditions, comparing it with both conventional absorption tomography and 2D TIE-Hom phase retrieval applied to projection images. The results indicate that the 3D post-reconstruction method generally achieves a modest improvement in noise suppression over existing PB-CT methods. It is also shown that potentially large computational gains over projection-based phase retrieval for multi-material samples are possible. In particular, constraining phase retrieval to a localized 3D region of interest reduces the overall computational cost and eliminates the need for multiple CT reconstructions and global 2D phase retrieval operations for each material within the sample.

1. Introduction

X-ray imaging is a cornerstone of modern medical imaging, with conventional 2D radiography and 3D computed tomography (CT) being common tools both in clinical and research domains (Bushberg et al., 2012). Phase-contrast imaging (PCI) is a specialized modality where image contrast is achieved by exploiting both the refractive and absorption properties of the imaged object. This technique has been widely studied and refined for over fifty years since the pioneering work of Bonse & Hart (1965). PCI has become a powerful approach for improved imaging of soft-tissue samples, which often exhibit poor contrast-to-noise ratios due to weak absorption contrast and photon fluence constraints in conventional absorption-based radiography (Bushberg et al., 2012). With the ongoing refinement of synchrotron and micro-focus X-ray sources, PCI is being adapted for medical use (Bravin et al., 2013; Tromba et al., 2016). Several PCI methods exist utilizing different mechanisms to encode phase information into projection images. In common use are analyser-based imaging (ABI) (Goetz et al., 1979; Gureyev & Wilkins, 1997; Davis et al., 1995; Nesterets et al., 2004) and grating-based imaging (GBI) (Cloetens et al., 1997; Momose et al., 2003; Weitkamp et al., 2008; Nesterets & Wilkins, 2008) which utilize crystals and gratings, respectively, in the experimental setup. In this work, we study the so-called propagation-based imaging (PBI), also known as the in-line method (Snigirev et al., 1995; Wilkins et al., 1996). In contrast to the other PCI methods, exclusive of combined approaches (Pavlov et al., 2004, 2005; Coan et al., 2005), PBI relies on the free-space propagation between the sample and detector for phase-contrast effects to manifest themselves as detectable intensity variations. Given this, PBI is simpler from an experimental perspective than other PCI methods. However, this simplicity is offset by more stringent requirements of the spatial coherence properties of the incident X-ray beam (Nugent, 2010).

PBI has been shown to produce enhanced image contrast in weakly absorbing objects such as biological samples, generally in the form of Fresnel diffraction fringes at the interfaces between different materials. However, material-specific quantitative information cannot be gleaned from PBI intensity images directly and requires the application of phase retrieval methods prior to CT reconstruction in order to recover the complex refractive index distribution n(r) = 1 − δ(r) + iβ(r) within the sample. Several methods of phase retrieval have been developed in PBI, with different restrictions imposed on the object and imaging system. The transport of intensity equation (TIE) based methods for phase retrieval from the work of Teague (1983) and refined by others (Cloetens et al., 1999; Bronnikov, 1999, 2002) are commonly used. To reconstruct the 3D distribution of n(r), most of these methods require multiple X-ray projections (at different propagation distances and/or X-ray energies) to be acquired at each angular position of the object, which can be difficult to achieve under experimental conditions with time and dose constraints.

A significant breakthrough was made by Paganin et al. (2002) who developed the so-called `homogeneous' TIE phase retrieval algorithm (TIE-Hom) that accurately reconstructs the complex refractive index of monomorphous objects (Turner et al., 2004). The algorithm requires only a single projection for each view angle and is robust to noise. As such, it has become the de facto standard for phase retrieval in PB-CT. The TIE-Hom algorithm makes use of a spatially uniform factor or delta-to-beta ratio, γ = δ/β, which defines the relative weight of the phase shift and absorption in the material of interest and results in the simplification of the reconstruction of the complex refractive index to

In reality, most samples do not consist of a single material, so use of TIE-Hom phase retrieval requires a compromised choice of γ for one particular material interface. Qualitatively, this choice results in the blurring of edges at the interfaces of materials where γ is overestimated and the retention of phase-contrast fringes for the underestimated case. Quantitatively, there are corresponding errors in the reconstructed distribution of n(r). These errors have been shown to be reduced by the collection of additional projections or utilizing suitable a priori information (Gureyev et al., 2013).

The TIE-Hom method has also been extended to enable quantitatively accurate phase retrieval of images containing non-overlapping projections of two materials (Gureyev et al., 2002) and subsequently of m materials (Beltran et al., 2010). In both cases, a priori information is available for values of γ_m for each material interface. The method proposed by Beltran et al. (2010) demonstrates that a composite 3D reconstruction of n(r) can be produced by m separate applications of TIE-Hom phase retrieval using a different value of γ_m for the projection set, each followed by CT reconstruction from which localized sub-volumes are spliced into the final reconstructed volume. The present work seeks to extend this method for phase retrieval of multi-material samples, differing from the previously described approach whereby the material-specific TIE-Hom phase retrieval step is performed as a post CT-reconstruction filtering operation in 3D. Importantly, this 3D version of TIE-Hom is applied to localized sub-volumes, again using a priori values of γ_m in each. The use of TIE-Hom phase retrieval discussed thus far has been applied to plane projections. This will now be referred to as pre-reconstruction 2D TIE-Hom (PreTIE-Hom2D). Multi-material phase retrieval has also been attempted by the application of a 3D correction filter operation applied to the reconstructed volume in addition to conventional 2D TIE-Hom (Ullherr & Zabler, 2015). Also, other non-TIE based methods of 3D phase retrieval have been developed (Vassholz et al., 2016; Ruhlandt et al., 2014; Maretzke et al., 2016). For clarity, our proposed method will be referred to as post-reconstruction 3D TIE-Hom (PostTIE-Hom3D) for which the derivation is given in Section 2. A numerical simulation framework has been created to simulate, evaluate and compare these methods and is discussed in Section 3.

2. Derivation – post-reconstruction 3D TIE-Hom phase retrieval

Consider the imaging system shown in Fig. 1. Let an object be illuminated by a monochromatic plane X-ray wave with wavelength λ and intensity I_in, I_in^1/2 exp(ikz) with k = 2π/λ. The image of the object is recorded on a position-sensitive detector located at a distance R downstream from the object. In the following we assume that the dimensions of the object are small compared with the source-to-object distance ρ and that ρ ≫ R. Interactions of the X-rays and object matter are described via the spatial distribution of the complex refractive index, n(r) = 1 − δ(r) + iβ(r), r = (x, y, z).

If the projection approximation (Paganin, 2006) is applied to the PB-CT experimental setup depicted in Fig. 1, the following equations allow for the calculation of the transmitted phase and intensity,

where P_θf(x′, y) represents the projection operator defined as

Here δ_D(x) is the Dirac delta function. The 2D Fourier transform operator

combined with the so-called filtered back projection (FBP) operator, (Natterer, 2001),

allows the construction

This inversion equation forms the mathematical basis of CT, permitting the reconstruction of the 3D distribution f(x, y, z) from a set of measured projections, P_θf(x′, y), at θ angles within the interval (0, π).

Figure 1 PB-CT experimental setup. The direction of the incident X-ray wave forms an angle θ ′ with the z axis, −π/2 ≤ θ ′ < π/2 and θ = θ ′ + π/2. The object and detector planes are located at z′ = 0 and z′ = R, respectively.

One can now utilize the FBP operator to obtain expressions to separately reconstruct the real and imaginary parts of the complex refractive index. By substituting rearrangements of equations (2) and (3) into equation (7), we obtain the following pair of equations for β and δ, respectively,

Conventional X-ray radiography and CT are generally concerned with measuring the intensity distribution of transmitted radiation in the object plane and reconstructing the imaginary part of the complex refractive index β, relating the absorption characteristics of the object to the measured intensity. We now seek to utilize the TIE to infer phase information from the visible diffraction fringes created upon propagation through a given distance. The finite-difference form of the TIE is given by

where and is an in-line phase-contrast image in the detector plane. This method generally requires multiple (at least two) intensity measurements (at different R) in order to solve equation (10) for the unknown phase distribution (at a given λ), which may be difficult or problematic in the context of an experimental implementation. Let us assume that the sample object is monomorphous such that a spatially independent (but energy-dependent) proportionality constant, γ = δ/β, holds for the complex refractive index (Paganin et al., 2002; Mayo et al., 2003). This assumption is valid, for example, for objects consisting of a single material and objects composed of light elements (z < 10) when irradiated with high-energy X-rays (60–500 keV). Utilizing this property, one can link phase and intensity by rearranging equations (2) and (3) into

applying to equation (11) to obtain

and then inserting equation (12) into the TIE [equation (10)] to arrive at

From the expression above, one notes that the intensity at z′ = R includes, in addition to the contact intensity, a phase-contrast term proportional to the 2D Laplacian of the contact intensity. A further simplification can be introduced if one considers the case of weak absorption, whereby 2kP_θ β(x′,y) ≪ 1 in equation (2),

Inserting this approximation into equation (13) and rearranging, gives

Moreover, TIE implies weak phase contrast (Gureyev et al., 2004), . This allows the approximation of the in-line contrast function, . Using this approximation and the Fourier space Laplacian identity, , the Fourier transform of equation (15) results in the following equation,

Inserting equation (16) into the FBP operator, equation (7) provides a `single-step' phase retrieval and CT reconstruction expression for monomorphous objects (Gureyev et al., 2006; Arhatari et al., 2007, 2012),

A similar result has been derived by Bronnikov (1999, 2002) for pure-phase objects with negligible absorption which can be seen to correspond to equation (17) with γ → ∞. Similarly, equation (17) reduces to conventional pure absorption CT when R = 0.

Using an explicit form of the FBP operator, equation (5) and applying , where , to both sides of equation (17) and utilizing the identity

such that with some algebraic manipulation the terms in the FBP operator cancel leading to the following expression for ,

thus concluding the derivation.

Moreover, the result in equation (18) can be rewritten in the Fourier domain,

which is a direct counterpart to equation (17) and makes use of the 3D Fourier transform operator

Equation (18) is the mathematical basis of our new approach for a 3D phase retrieval method for weakly absorbing monomorphous objects. As is seen, the phase retrieval step is decoupled from the FBP operator, implying that this step can be performed in 3D after the conventional CT reconstruction. A similar derivation for pure-phase objects is described by Baillie et al. (2012). One of the more interesting aspects to the form of equation (18) relates to the ability of performing phase retrieval localized to a region of interest (ROI). This property potentially gives an advantage over pre-reconstruction phase retrieval techniques for samples which may contain a range of different materials of interest whereby it is often difficult to optimize global parameters to achieve optimal contrast across the entire sample (Gureyev et al., 2013). We will illustrate that with localization we can apply a specific value of γ chosen to `focus' on a desired material composition inside a 3D region, thus giving a quantitatively accurate phase retrieval within that region. A spatially localized form of phase retrieval also leads to some potential computational gains over existing methods, such as the ability to divide and parallelize phase retrieval over the sample reconstruction. The use of fast discrete Fourier or GPU based filters for implementing the 3D TIE-Hom filter on sub-regions would make this phase retrieval method computationally efficient, even for large datasets.

3. Numerical experiment setup

3.1. X-ray CT simulation framework

In order to evaluate the accuracy and characteristics of PostTIE-Hom3D, a computational simulation framework was constructed for conventional absorption CT and PB-CT workflows. This framework allows for the definition of an analytical 3D model consisting of multiple simple geometric primitives in space representing a material defined by its complex refractive index at a given X-ray energy. With such a model, the workflow can be used to generate a volume image in addition to contact or PBI projections for a given number of rotation angles with specified resolution and photon statistics. These simulated projections can then be subsequently reconstructed, including phase retrieval as part of the processing pipeline and quantified with a range of metrics. The framework was primarily constructed using the ITK image processing toolkit (Ibáne~z et al., 2005) and uses elements of the Astra-toolbox (van Aarle et al., 2016), X-TRACT (Gureyev et al., 2011) and RTK (Rit et al., 2014) packages for CT-reconstruction routines.

3.2. Multi-material numerical phantom

The numerical phantom model defined for these simulations consists of an air-filled 1024 µm × 1024 µm × 256 µm cuboid containing a central 700 µm-diameter cylinder consisting of breast tissue with the composition described by Hammerstein et al. (1979) and serving as the background material. Embedded within the central cylinder are four smaller 80 µm-diameter cylinders of different materials distributed azimuthally at a fixed distance from the origin. A schematic representation of the phantom and corresponding material properties for the objects at an X-ray energy of 20 keV are shown in Fig. 2 and Table 1. Values obtained for material specific β and δ were calculated using the web-based tool (Gureyev et al., 2011). The relative delta-to-beta ratio versus breast tissue, γ_rel, has been calculated for each material as γ_rel = (δ − δ_breast)/(β − β_breast), where δ_breast and β_breast are the values for breast tissue.

Table 1
Material properties at E = 20 keV

Object	Material	β × 10^10	δ × 10^7	γrel (breast tissue)
1	Air	0.004641	0.00625	1681
2	Breast tissue	3.48	5.85	1
3	Weddellite (CaC2O4·2H2O)	35.17	10.36	141
4	Paraffin (C31H64)	1.93	5.35	322
5	Adipose	2.54	5.36	523
6	Tumour	3.57	5.93	950

Figure 2 Top view of the numerical phantom, illustrating dimensions and materials as defined in Table 1.

The materials modelled within the inner cylinders include organic tumorous and adipose tissue. The calcite weddellite (CaC₂O₄·2H₂O, mass density 1.94 g cm⁻³) represents one of the primary components of type I and type II breast calcifications (Ghammraoui & Popescu, 2017). Additionally, the hydro­carbon paraffin (C₃₁H₆₄, mass density 0.9 g cm⁻³) was included as it can be used to enclose biopsy samples and has X-ray absorption and refractive properties similar to those of the other organic materials.

3.3. Simulations

This research seeks to evaluate the numerical implementation of PostTIE-Hom3D derived in Section 2 and compare its performance with conventional absorption CT and the Beltran et al. (2010) PreTIE-Hom2D method for the simulated multi-material phantom defined in Section 3.2.

3.3.1. Projection simulation

Distributions of transmitted intensity and phase shift from a simulated X-ray source (at 20 keV) were computed for N_p sample rotation angles analytically from the phantom model definition in Section 3.2 using equations (2) and (3). The complex amplitude of the transmitted wave was then calculated over a 2D plane sampled at a given resolution. In the case of the simulations for this article, a relatively fine sampling interval of 0.25 µm was chosen for the generation of initial projections. Due to the invariance of the phantom along the rotation y axis, we can employ the simplification only requiring the generation of a single-row projection for the phantom at each angle step, thus resulting in 1D plane-projection images of 8192 pixels.

For PBI projections, the Fresnel propagation operator (Paganin, 2006) implemented as a 2D Fourier filter was applied to each 2D projection, with the transfer function, FG(ξ′,η) = exp[−iπλR(ξ′²+η²)]. In the case of contact projections, the above propagation step was skipped leading to the final step of projection simulation where a discrete Gaussian smoothing filter was applied (with the variance σ² = 2 µm²) to simulate finite resolution of the detector. This step led to the smearing of sharp edges at the interfaces between objects and to reducing potential aliasing artefacts of subsequent phase retrieval and CT reconstruction. Finally, projections were sampled with a finite square aperture of 2 µm resulting in N_p 1D row intensity projection images of 1024 pixels (32-bit real).

3.3.2. Dose and noise

In this work we seek to model a fixed total radiation dose per scan utilizing a specified exposure time per angle. To quantify the dose, we introduce the simulation parameter `total photons per pixel' (TPP) as the total number of incident photons per pixel per scan. In practice, we simulate noisy projections by applying a Poisson distribution with a known mean number of photons per individual pixel.

3.3.3. Pre-reconstruction 2D TIE-Hom phase retrieval

For the application of the PreTIE-Hom2D method, a set of N_p 2D projections was constructed by vertically stacking copies of the simulated, propagated and binned 1D row projections generated as per Section 3.3.1. Poisson noise was then generated as described in the previous section, followed by 2D TIE-Hom phase retrieval [equation (21)] applied to each 2D projection using the corresponding phantom-material-specific value of γ_m,

3.3.4. CT reconstruction

To recover the imaginary component β of the complex refractive index from the thus far simulated projections, the `gold-standard' FBP CT-reconstruction algorithm, equation (7), was applied. Prior to performing the actual reconstructions, standard background (flat-field) correction was performed before the −ln transform. In the case of the PreTIE-Hom2D multi-material method and as discussed in Section 3.3.3, N_p 2D phase-retrieved projections are reconstructed into a 1024³ voxel volume. To construct the final spliced multi-material volume, 2D TIE-Hom phase retrieval followed by FBP CT-reconstruction was performed for each phantom material separately with the composite volume constructed by inserting sub-volumes enclosing each cylinder. For PostTIE-Hom3D, N_p 2D simulated projections were constructed following the approach described in Section 3.3.3, without the phase retrieval step and then reconstructed with FBP producing a 1024³ voxel volume to which post-reconstruction 3D phase retrieval was applied.

3.3.5. Post-reconstruction 3D TIE-Hom phase retrieval

PostTIE-Hom3D phase retrieval is implemented as a 3D Fourier filter [see equation (19)] and applied to sub-volumes corresponding to ROIs contained within the whole CT reconstructed volume. Successful application of PostTIE-Hom3D requires that the selected 3D ROI meets the following two criteria. Firstly, the ROI should fully contain the single-material object under investigation. Ideally, this region should not be `polluted' with the inclusion of other objects or reconstruction artefacts which will lead to further undesirable artefacts in the phase-retrieved sub-volume. Secondly, the ROI should be chosen to consider the width of the 3D TIE-Hom point-spread function (PSF), P_TIE, the Fourier transform of which is defined as

Here, A = πγλR is a positive constant (we restrict our consideration to the case of positive γ). The corresponding real-space expression is

Here, l_TIE = A^1/2/2π and r = (x² + y² + z²)^1/2, with the standard deviation of this distribution equal to . To overcome situations where the criteria are unable to be met with an appropriately large ROI due to neighbourhood constraints, it is possible to artificially enlarge the ROI by zero-padding to the required dimensions. Zero-padding to a power of two is also commonly applied in implementations of fast Fourier transform (FFT) based filters for optimal performance. For the simulations performed in this article, the dimensions and locations of the cylindrical phantom objects are known explicitly, therefore a 128 µm × 128 µm × 128 µm cubic ROI was chosen for each of the four cylinders which adequately extends beyond the 80 µm cylinder diameter and easily satisfies the requirement for P_TIE. For example, given an X-ray energy of 20 keV, adipose tissue with γ ≃ 523 (Table 1) and a propagation distance of R = 60 mm results in the calculated characteristic length scale 1_TIE ≃ 12µm, corresponding to only several pixels at the simulated detector resolution.

3.3.6. Evaluation metrics

To quantify and compare the performance of the evaluated phase retrieval methods, two metrics have been selected, the contrast-to-noise ratio (CNR) (Gureyev et al., 2014) and a universal image quality index (UIQI) introduced by Wang & Bovik (2002). CNR is defined as follows,

where V is the volume of the object ROI (in voxels), 〈β_o〉 and σ_o² are the mean and variance of β within the object ROI and 〈β_b〉 and σ_b² are the mean and variance for a `background' ROI devoid of any phantom objects. In these simulations, this background region corresponds to a similarly sized region as the object ROI and is located outside all object ROIs in the reconstructed volume. CNR is calculated to compare the associated statistical gain between absorption only contact CT and the two phase retrieval methods which will highlight the respective noise suppression properties. CNR gain, G_CNR, is calculated as the ratio of the CNR in the CT reconstructions using phase-retrieved and contact projection data, for same-object and background ROIs,

The second metric UIQI, where

requires a reference image for comparison, which is provided by the same ROI extracted from the numerical phantom model, with 〈β_r〉, σ_r² and σ_or representing the mean, variance and covariance between the reference and sample ROIs, respectively. The UIQI produces a single numerical index that combines three factors: loss of correlation, luminance distortion and contrast distortion which its authors suggest permit the measurement of information loss as opposed to the quantification of error with other metrics (Wang & Bovik, 2002).

3.4. Results

Simulation parameters have been selected to be consistent with an experimental PB-CT imaging scenario representative of a small biological biopsy. For the simulations presented and discussed here, the following series of fixed parameters have been selected: an incident monochromatic X-ray beam with E = 20 keV and pixel size h = 2 µm (in both directions). A total of 900 projections were generated, N_p = 900 corresponding to an angular step of 0.2°. This number conforms to the Nyquist sampling condition N_p > (π/2)d/h, for estimating the lower-bound for the number of required projections to avoid sub-sampling, with d being the reconstruction diameter (Hsieh, 2009). A propagation distance of R = 60 mm was chosen to be just below the limit for the TIE-Hom validity condition of λR/h² = 1, such that good phase contrast should be obtained in the simulated propagated projections. A range of photon statistics were applied to the simulated projections as described in Section 3.3.4. An exponentially increasing series of values in terms of the TPP and relative noise as a percentage were used, 1 × 10⁴ (1.0%), 1 × 10⁵ (∼0.32%), 1 × 10⁶ (0.1%), 1 × 10⁷ (∼0.032%) and 1 × 10⁸ (0.01%). Notably, noise-free and contact projections were also simulated for comparative purposes.

The grid of 16 images presented in Fig. 3 comparatively illustrates PB-CT and the application of phase retrieval to FBP reconstructions of the simulated phantom for the imaging scenario described above in the case of 0.1% noise. Each image displays the central 2D slice through a 128 µm³ cubic ROI of each material cylinder. Grid rows represent the materials: weddellite, paraffin, adipose and tumorous tissue; whilst columns display contact CT (R = 0 mm), PB-CT (no phase retrieval), PreTIE-Hom2D and PostTIE-Hom3D for the propagation distance R = 60 mm. Visually inspecting the reconstructions for weddellite [Figs. 3(a)–3(d)] some object contrast is achieved in the contact reconstruction [Fig. 3(a)], albeit with a reasonably high level of noise. Phase-contrast fringes are clearly visible in the reconstruction of the propagated case without phase retrieval [Fig. 3(b)] and similar high-quality results are shown for both TIE-Hom phase retrieval methods in Figs. 3(c) and 3(d). The results for paraffin [Figs. 3(e)–3(h)] and adipose tissue [Figs. 3(i)–3(l)] are visually very similar, in both cases no apparent object contrast is perceptible in contact reconstructions. Again, noted is the presence of phase-contrast fringes in the non-phase retrieved images in Figs. 3(f) and 3(j), although they are less defined and contain significantly more visible noise than the weddellite case. Both variants of phase retrieval produced visually similar results for the two materials and achieve relatively good contrast and reduced levels of noise such that the object is clearly distinguishable from the background. In the case of tumorous tissue [Figs. 3(m)–3(p)] the results are far less revealing. It is worth noting that, at E = 20 keV, tumorous tissue has absorption and refractive properties more similar to the background breast tissue (Table 1) than the other materials in the phantom, making the task of distinguishing the two regions after reconstruction significantly more difficult. This difficulty is clear in the non-phase retrieved images; both contact and propagated images in Figs. 3(m) and 3(n) exhibit only noise. Notably, the success in the application of phase retrieval for tumorous tissue is not convincing when viewing a single reconstructed slice as in the case in Figs. 3(o) (PreTIE-Hom2D) and 3(p) (PostTIE-Hom3D), where little apparent contrast or structure is present for one to identify the reconstructed object.

Figure 3 Comparison of FBP reconstructions for the simulated phantom across the four materials with and without phase retrieval. Each image displays the central reconstructed slice in the volume ROI for the corresponding material with 0.1% noise. Rows display results for weddellite, paraffin, adipose and tumorous tissue. Columns display reconstructions for contact-CT (R = 0 mm), PB-CT (no phase retrieval), PreTIE-Hom2D and PostTIE-Hom3D (R = 60 mm).

3.4.1. Reconstruction line profiles

To gain insight into the relative quantitative merits of the two-phase retrieval methods, line profiles of the reconstructed central slice in each of the four materials are presented in Fig. 4. Each plot presents the line profile of β values over a central line extending 40 µm horizontally across the ROI. Separate plots for the reference computational model and reconstructed values after phase retrieval with PreTIE-Hom2D and PostTIE-Hom3D for R = 60 mm and 0.1% noise are plotted.

Figure 4 Line-profile plots comparing the imaginary part of the complex refractive index, β, across the central reconstructed slice between the phantom model and PreTIE-Hom2D and PostTIE-Hom3D phase retrieval for R = 60 mm and 0.1% noise. (a) Weddellite, (b) paraffin, (c) adipose and (d) tumorous tissue.

Inspection of Fig. 4(a) (weddellite) shows a similar reconstruction result with phase retrieval in both cases resulting in the elimination of diffraction fringes which are clearly visible in the contact reconstruction [Fig. 3(c)]. Good contrast and suppression of noise have been achieved, giving rise to values of β in agreement with the model. Notably, there is a similar degree of smoothing of the edges for both methods which may be attributed to the finite resolution of the imaging system and interpolation of projection data during CT reconstruction. The noise-suppressing properties of PostTIE-Hom3D compared with PreTIE-Hom2D are illustrated for the less-absorbing materials paraffin [Fig. 4(b)] and adipose tissue [Fig. 4(c)]. Here, both profiles and images are similar, both materials exhibit comparable absorption and refractive properties at the simulated X-ray energy as shown by the values in Table 1. Again, for these simulation parameters, PostTIE-Hom3D displays qualitatively improved noise reduction for these materials. For the final material, tumorous tissue [Fig. 4(d)], the profiles indicate only limited success of both PreTIE-Hom2D and PostTIE-Hom3D in improving object contrast in this case. The object is barely perceptible after phase retrieval although it may be suggested that PostTIE-Hom3D is marginally less noisy.

3.4.2. Evaluation with respect to noise

The evaluation of the noise suppression, contrast enhancing properties and quantitative accuracy of both methods as displayed in the previous section needs to be viewed with respect to varying levels of Poisson noise. Fig. 5 presents plots of G_CNR and UIQI (Section 3.3.6) for each material over a range of values of Poisson noise on a logarithmic scale. In the case of UIQI, an additional plot for `contact' reconstructions without any phase retrieval is displayed as a point of reference. The relative noise suppression characteristics between PreTIE-Hom2D and PostTIE-Hom3D are quantified by the calculated values of G_CNR and UIQI, respectively. Almost universally, both methods exhibit improved CNR at all noise levels for all four materials with PostTIE-Hom3D generally superior. However, some variation in the overall trends are observed across different materials. For example, in Figs. 5(a) (weddellite) and 5(b) (paraffin) we see that the reported values of G_CNR are both maximal at the highest simulated noise level, TPP = 1 × 10⁴ (1.0%), and decrease monotonically approaching 1 (no gain) at the lowest simulated noise level, TPP = 1 × 10⁸ (0.01%). For adipose tissue [Fig. 5(c)], one notes quite different behaviour with both methods displaying an initial increase in gain, peaking at TPP = 1 × 10⁵ (∼0.32%) and then decreasing, approaching 1 at 0.01% noise. In this case, it is also noted that PostTIE-Hom3D significantly outperforms PreTIE-Hom2D with respect to G_CNR until their values align at around TPP = 1 × 10⁷ (∼0.032%). For tumorous tissue [Fig. 6(d)], one sees some similarity with that of adipose with a general increase in G_CNR with a peak shifted towards the lower end of the simulated noise spectrum at TPP = 1 × 10⁷ (∼0.032%) and then falling away. Again, PostTIE-Hom3D exhibits a near uniform improvement in G_CNR over PreTIE-Hom2D at noise levels approaching the peak, after which both methods perform similarly, as previously noted.

Figure 5 Plots of GCNR and UIQI for R = 60 mm and varying levels of noise. (a) Weddellite, (b) paraffin, (c) adipose tissue and (d) tumorous tissue.

Figure 6 Comparison of FBP reconstructions of the simulated phantom within the ROI containing tumorous tissue with and without phase retrieval and varying noise levels. Each image displays the central reconstructed slice in the volume ROI. Rows display results for noise levels of 0.1%, ∼0.03% and 0.01%. Columns display contact-CT (R = 0 mm), PB-CT (no phase retrieval), PreTIE-Hom2D and PostTIE-Hom3D (R = 60 mm) reconstructions.

Turning now to the plots for UIQI, as mentioned in the definition (Section 3.3.6), this metric attempts to quantify `image quality' from a broader perspective than the purely statistical nature of CNR. Its use of a reference image in the form of the numerically accurate model ROI implies that it may provide values more representative of the quantitative deviation of the subject image from the reference. For two of the less absorbing materials, paraffin [Fig. 5(b)] and adipose [Fig. 5(c)], one observes very similar results with a monotonically increasing value of UIQI as TPP increases, with PostTIE-Hom3D giving a slightly higher quality result for higher noise levels until converging with PreTIE-Hom2D at around 0.01% noise. For these materials, one also notes that UIQI for `contact' reconstructions remains relatively uniform until around 0.032% noise, from where the quality metric increases with a corresponding decrease in noise. With tumorous tissue [Fig. 5(d)] the UIQI for both phase retrieval methods dips below that of contact reconstructions before increasing in line with decreasing noise as with the other materials. This initial degradation in UIQI at higher noise levels for tumorous tissue can be attributed to the relatively little object information available to `retrieve' relative to noise such that phase retrieval introduces a greater deviation from the model. The more highly absorbing weddellite sample [Fig. 5(a)] illustrates different behaviour from the others in that UIQI for both phase retrieval methods remains relatively uniform and similar irrespective of the noise level. In contrast, UIQI for contact reconstructions of weddellite increases monotonically (as expected) with decreasing noise, exceeding the corresponding values for both PostTIE-Hom3D and PreTIE-Hom2D for TPP > 1 × 10⁷.

The relative improvements in image quality due to phase retrieval in comparison with contact CT as presented quantitatively via G_CNR and UIQI as noise levels are varied are shown in Fig. 6. Here, in a similar form to Fig. 3, central 2D FBP reconstructed slices for the ROI containing tumorous tissue are displayed for contact and PB-CT with phase retrieval for three different levels of noise (rows) corresponding to 0.1%, ∼0.03% and 0.01%. Columns display contact-CT, PB-CT (no phase retrieval), PreTIE-Hom2D and PostTIE-Hom3D for the propagation distance R = 60 mm. Figs. 6(a)–6(d) show the same images as Figs. 3(m)–3(p), corresponding to 0.1% noise. Reviewing the images for contact-CT in the first column for each noise level [Figs. 6(a), 6(e) and 6(j)] one observes no perceivable contrast; all reconstructed slices appear to display uncorrelated noise. For the second column, PB-CT without phase retrieval, one notes at ∼0.03% [Figs. 6(f)] and 0.01% noise [Fig. 6(j)] the emergence of a diffraction fringe at the boundary interface between the tumorous and breast tissue in addition to reconstruction `streak' artefacts which are most likely due to sharp edges of the neighbouring weddellite object and are similar in appearance to those for adipose tissue and paraffin at the higher noise level of 0.1% as seen earlier in Figs. 3(f) and 3(j). Turning now to the phase retrieval results in columns three and four, one sees that both PreTIE-Hom2D and PostTIE-Hom3D are able to produce results with enough contrast that the circular region of tumorous tissue is clearly distinguishable from the background breast tissue for ∼0.03% noise [Figs. 6(g)–6(h)] and slightly higher perceivable contrast at 0.01% noise [Figs. 6(k) and 6(l)]. Again, noted is the appearance of streak artefacts and a marginally visually enhanced result in the case of PostTIE-Hom3D than PreTIE-Hom2D.

Overall, the results in this section illustrate the ability of TIE-Hom-based phase retrieval in conjunction with PB-CT to reveal useable image contrast in the presence of noise where conventional contact-CT fails. Moreover, this research demonstrates that the level of contrast achieved is a function of the material, noise statistics and intrinsic properties of the imaging system and is consistent with theoretical derivations for the effect of phase retrieval in PB-CT (Nesterets & Gureyev, 2014).

3.4.3. Evaluation with respect to propagation distance

With some insight gained as to the characteristics of the phase retrieval schemes in the presence of varying levels of simulated noise and hence dose, the focus now becomes an evaluation of the response of varying propagation distance, which serves as the key experimental parameter in PB-CT. Several studies (Nesterets & Gureyev, 2014; Kitchen et al., 2017) have demonstrated that increasing propagation distance in combination with TIE-based phase retrieval results in significantly improved signal-to-noise ratio in resulting images. Such improvements permit the use of a smaller radiation dose, producing images of similar quality to those produced without PB-CT methods at higher doses.

For the simulations previously described at a fixed propagation distance, R = 60 mm was used and corresponded to just below the distance for maximum phase contrast in the TIE regime defined by λR/h² = 1. For the analysis in this section the noise level is fixed at 0.1%, TPP = 1 × 10⁶ and the propagation distance varied with 10 mm ≤ R ≤ 70 mm at 10 mm increments. Fig. 7 presents plots of G_CNR and UIQI across this range of propagation distances. With respect to G_CNR, one sees relatively consistent trends across the four materials with gain generally increasing as propagation distance increases with PostTIE-Hom3D, recording greater levels than PreTIE-Hom2D. The magnitude of the gain does vary with material, with weddellite [Fig. 7(a)] achieving modest levels of the order 2.3 ≤ G_CNR ≤ 2.5 for PreTIE-Hom2D and 2.4 ≤ G_CNR ≤ 2.8 for PostTIE-Hom3D. The other less-absorbing materials exhibit greater levels of G_CNR, with paraffin [Fig. 7(b)] exhibiting around 10 ≤ G_CNR ≤ 32 (PreTIE-Hom2D) and 15 ≤ G_CNR ≤ 36 (PostTIE-Hom3D). Adipose tissue [Fig. 7(c)] reports slightly higher levels with 15 ≤ G_CNR ≤ 40 for PreTIE-Hom2D and 20 ≤ G_CNR ≤ 50 for PostTIE-Hom3D. In this particular case, one also notes that there is a slight reduction in the magnitude of G_CNR for propagation distances where R > 50 mm, which is the likely result of reconstruction artefacts beginning to dominate over the random noise. The results for tumorous tissue [Fig. 7(d)] differ again from previous cases with G_CNR for both phase retrieval methods showing an initial rise then dip in the vicinity of 10 mm ≤ R ≤ 40 mm followed by a sustained increase for R > 40 mm where the magnitude for PreTIE-Hom2D is 15 ≤ G_CNR ≤ 60 and 45 ≤ G_CNR ≤ 80 for PostTIE-Hom3D. The graphs of G_CNR over all four materials reveal that PostTIE-Hom3D consistently produces greater contrast over PreTIE-Hom2D by a relatively fixed, non-converging magnitude as the propagation distance increases.

Figure 7 Plots of GCNR and UIQI for 10 mm ≤ R ≤ 70 mm and TPP = 1 × 10^6. (a) Weddellite, (b) paraffin, (c) adipose tissue and (d) tumorous tissue.

Examination of the graphs for UIQI shows that again both paraffin [Fig. 7(b)] and adipose tissue [Fig. 7(c)] trend similarly with increasing propagation distance with nearly identical values in the range 0.2 ≤ UIQI ≤ 0.5. As evident with G_CNR, PostTIE-Hom3D tends to outperform PreTIE-Hom2D by a fixed magnitude with UIQI levelling off from around R ≥ 60 mm, after which the theoretical validity conditions of TIE-Hom are exceeded, leading to a loss of fidelity in the phase retrieval process. Tumorous tissue [Fig. 7(d)] shows some inconsistency in the evaluated UIQI over both methods as was the case with G_CNR. This is likely due to the lack of available phase contrast in the projections to adequately reconstruct the object as seen visibly in the line profiles in Fig. 5, and is quantified by the relatively low values of UIQI compared with the other materials. Despite this, PostTIE-Hom3D produces a UIQI nearly double that of PreTIE-Hom2D.

The plot in Fig. 7(a) for weddellite shows distinctly different behaviour from the other materials. In this case one observes a monotonic decrease in UIQI for both methods as the propagation distance increases, opposite in behaviour to the others. Additionally, one notes that PreTIE-Hom2D outperforms PostTIE-Hom3D from around R > 20 mm where the latter decreases more rapidly as the propagation distance increases. Fig. 8 illustrates line profiles plotting PostTIE-Hom3D β values for weddellite in the interval 10 mm ≤ R ≤ 70 mm with TPP = 1 × 10⁶ in addition to the model and contact profiles. Evident from these plots is the ever-increasing degree of smoothing produced around the edges of the object as the propagation distance increases, resulting in an overall increased deviation from the reference model profile and degrading UIQI relative to the model as shown in Fig. 7(a). Such over-smoothing can be attributed to the validity conditions of the TIE being exceeded for weddellite at increasing propagation distances.

Figure 8 Line profiles of β values for the phantom model, contact and PostTIE-Hom3D phase retrieval for weddellite over propagation distances of 10 mm ≤ R ≤ 70 mm.

Overall, the results show that the PB-CT in conjunction with the phase retrieval methods investigated result in significant improvements to image quality over contact-CT imaging. The results are consistent with those reported by Gureyev et al. (2014), illustrating improvements in contrast-to-noise of several orders of magnitude with increasing propagation distances. Moreover, such improvements correspond with an even greater relative reduction in the radiation dose required to obtain the equivalent quality to conventional contact imaging.

3.5. Computational analysis

Up to this point, this article has considered the merits of the two phase retrieval methods discussed purely in terms of their imaging characteristics. This section will now examine their performance from a computational perspective, which is increasingly relevant in the context of real-life applications of PB-CT. Often, a realistic experimental scenario involves the use of synchrotrons or laboratory-based micro-CT systems leading to large datasets requiring dedicated high-performance computing (HPC) infrastructure to process within feasible timeframes. Such conditions lead to the need to consider the computation costs of the methods employed. Computationally, the application of the multi-material method by Beltran et al. (2010) for a sample consisting of N_m distinct materials over N_p projections requires N_m × N_p 2D phase retrieval operations in addition to N_m FBP CT reconstructions. The computational complexity of the FBP algorithm is well studied (Natterer, 2001) and is dominated by the back-projection step. If one assumes optimal sampling conditions, the total work for reconstructing an N_w² × N_h volume is proportional to N_hN_w³. TIE-Hom phase retrieval in 2D projection space is generally implemented as a 2D Fourier filter with total computational cost for N_p projections proportional to N_pN_wN_h ln(N_wN_h). From these two components it is clearly seen that the dominant computational element of the Beltran and coworkers multi-material method is the N_m FBP CT reconstructions, thus computationally bound by the number of materials to be examined and the sample dimensions, N_mN_hN_w³. In contrast, the PostTIE-Hom3D method only requires a single FBP CT reconstruction with the application of phase retrieval implemented as a set of N_o 3D Fourier filtering operations, where N_o is the number of localized object ROIs examined. The complexity of each 3D Fourier filter is proportional to N³_ROI ln(N³_ROI), where N³_ROI is assumed to be a cube of uniform dimensions representing the ROI. Therefore, the PostTIE-Hom3D method is computationally bounded by the number of objects examined and their size, N_o and N_ROI, respectively. It is evident from these observations that PostTIE-Hom3D requires significantly less computation when both the number of objects and their size relative to overall sample under investigation is constrained. Conversely, the multi-material PreTIE-Hom2D method becomes computationally expensive as the number of materials examined increases. In terms of memory use, PostTIE-Hom3D requirements correspond to that for computing the forward and inverse 3D FFTs for each ROI sub-volume. If an in situ FFT algorithm is considered, then the memory required mirrors that of the memory used to store the volume itself. Moreover, the Fourier filtering operation that computes the phase retrieval can also be performed in-place, thus conserving the memory footprint of the total operation. For example, a 128³ pixel ROI sub-volume would require a total of 128³ memory elements to compute the in-place 3D FFT for the phase retrieval step, corresponding to around 8 KB if using 4-byte floating point data.

A common experimental scenario is the case of a single-material and single-phase retrieval ROI that extends over the full tomography volume. This example represents the pathological `worst case' for the PostTIE-Hom3D method whose computational gains are derived from the ability to perform phase retrieval on one or more sub-volumes in parallel with dimensions smaller than the full tomographic volume. As such, in this case PostTIE-Hom3D would offer no benefit over conventional pre-CT phase retrieval. The method is specifically targeted at applications where the full volume may contain isolated features of different materials that need to be accurately segmented or classified as part of a CT workflow and as such require material-specific phase retrieval.

One scenario where PostTIE-Hom3D offers a significant advantage over pre-reconstruction projection-based phase retrieval is that where the method is applied to a small ROI and varying the γ value with the aim of obtaining a subjectively optimal value of γ for a certain imaging task (for example, for feature detectability). This calculated value can be subsequently used in a global or localized CT-reconstruction workflow. Obtaining γ in this way eliminates the need for projection-wide phase retrieval followed by full volume CT reconstruction as in pre-CT phase retrieval, thus reducing the computational cost substantially.

Synchrotron CT for mammography is another real-life case whereby successful quantitative detection and differentiation between glandular, adipose, tumorous tissue and micro-calcifications are required. In such a workflow, an initial classification step is required to define the ROIs and material-specific gamma value for phase-retrieval; this may be a manual or automated process. In the mammography usage case described, the number of ROIs detected may be of the order of tens or hundreds for which PostTIE-Hom3D can be computed in parallel. In the case of HPC systems, this may be achieved in a distributed approach across many compute nodes and/or cores etc. Moreover, the PostTIE-Hom3D algorithm can be fully implemented on GPU architectures such that both the forward/inverse FFTs in addition to the filter can be performed as a completely GPU resident operation subject to available GPU memory to contain the ROI volume. For larger volumes, a more sophisticated distributed multi-GPU approach is theoretically possible. Such a GPU implementation would achieve real-time phase retrieval.

4. Conclusions

This research has derived a variant of TIE-Hom phase retrieval that can be applied directly to localized 3D regions of interest, consisting of isolated single-material objects within a greater reconstructed volume of the distribution of the complex refractive index. Each ROI requires a priori information relating to the absorption and refractive properties of the contained material. This method allows for efficient and accurate reconstruction of multi-material samples within the TIE regime, with several marked benefits over alternative approaches. A simple numerical framework for PB-CT X-ray imaging has been constructed allowing for the method to be simulated and compared with the approach developed by Beltran et al. (2010) on a synthetic phantom for a range of experimental parameters.

It is shown numerically using the contrast-to-noise and universal image-quality index metrics that this new method achieves improved noise suppression, contrast enhancement and overall quantitatively correct results compared with the contact absorption-only CT and TIE-Hom phase retrieval methods applied to projections prior to CT reconstruction.

It is also shown that the proposed 3D TIE-Hom method offers significant computational efficiencies for localized and multi-material samples where performing post-reconstruction phase retrieval eliminates the need for multiple CT reconstructions, representing the most computationally expensive element. The post-reconstruction nature of the method also offers several practical benefits for tomographic workflows. Namely, by having the phase retrieval step after reconstruction, it allows for the retention of existing absorption-tomography experimental configurations, where the insertion of extra processing steps between acquisition and reconstruction may be difficult. In contrast, 3D TIE-Hom phase retrieval may be performed optionally on a previously acquired and reconstructed dataset without having to repeat the full CT-reconstruction workflow. Finally, a practical application of this approach has been suggested by Paganin (2015), whereby PostTIE-Hom3D phase retrieval is applied to a conventional CT-reconstruction of phase-contrast projections, interactively adjusting the TIE-Hom γ parameter and observing the result in order to obtain subjective `focus' for a specific material within the ROI. It is envisaged that such a method, if implemented in a software application using modern GPU hardware, could allow real-time combined reconstruction and visualization of multi-material objects.