1. Introduction

X-ray phase contrast tomography is a technique capable of reconstructing the phase map of an object, often leading to an increased image contrast when compared with conventional, absorption-based X-ray computed tomography (CT), with the added benefit of a potential reduction in radiation dose [1]. This makes it a desirable tool for the investigation of low-absorbing samples, which naturally suffer from weak contrast in conventional CT, and particularly those constrained by dose limits, e.g. biological samples. Many phase contrast methods require a highly coherent X-ray source and are therefore limited to synchrotrons, however, several techniques have been implemented with laboratories sources [2]. In propagation-based imaging (PBI), although the technique relies on spatial coherence to detect interference patterns [3], the development of microfocus sources enabled its implementation outside synchrotron facilities [4,5]. While microfocus sources are mostly associated with reduced flux, therefore leading to long exposures, promising new technologies such as the liquid-metal-jet source can overcome this limitation [6]. In the grating interferometry method, phase sensitivity is achieved by employing two gratings to create and detect the Talbot self-imaging effect [7]. Spatial coherence is therefore required, however, the technique can be made compatible for use with extended laboratory sources, provided that a third, “source” grating is used [8, 9]. Another technique, the single grid method, employs an absorption (or phase) grid to create a reference beam, which is later compared to the detected beam in the presence of a sample, in order to extract information about the sample’s complex refractive index [10]. While the method was first designed for use with synchrotron radiation, a recent study explored the possibility of adapting it for use with microfocus, polychromatic sources to perform dynamical imaging [11]. The edge illumination (EI) technique, considered in this study, is a non-interferometric, achromatic imaging method [12], which can be used with either synchrotron radiation or with polychromatic, extended sources (with spot size up to ∼100 μm) [13–16].

The working principles of EI have been discussed in detail previously [17, 18]; in general, in the laboratory, phase sensitivity is achieved by employing sets of apertures (masks), placed behind and in front of the sample. While the first mask divides the incoming, divergent beam into physically separated beamlets, the second mask is placed in contact with the detector, creating insensitive regions between neighboring pixels. By slightly misaligning the masks such that only part of each beamlet reaches the active area of the pixel, refraction due to the sample can be detected, since the latter will deflect the beamlet towards or away from the pixel, therefore leading to an increase or decrease in the number of detected photons. [18,19]

Until recently, in EI, at least two frames per projection were required to extract an object’s phase information, usually acquired at different positions of the sample-mask [20]. Diemoz et al. have recently developed a phase retrieval algorithm where a single image per projection is sufficient [21]. In a similar way to the algorithm developed by Paganin et al. for PBI [22], Diemoz’s algorithm for EI relies on several assumptions, the main one being that the sample is quasi-homogeneous (i.e. it has a constant ratio between the real and imaginary parts of its complex refractive index). Using projection data acquired with synchrotron radiation, the algorithm has been shown to provide quantitative results when applied to a homogeneous sample, and enhanced image quality when imaging complex biological samples [21].

Translation of this single-image retrieval algorithm to CT scans would have significant impact, as it would considerably reduce scan times and simplify the acquisition sequence, as well as possibly lead to a dose reduction. Indeed, recent developments have made the algorithm suitable for use on CT data acquired in a laboratory environment, with an extended and polychromatic X-ray source [23]. While high image quality was observed for a range of samples which can be considered approximately homogeneous (e.g. rat heart), many biological samples of interest are not quasi-homogeneous (e.g. contain both soft-tissue and bone), and applying the algorithm to such samples will result in under or over-retrieved characteristic features (dark/bright fringes, blurring). For example, in a paper by Diemoz et al., this can be seen as blurring around a chalk detail embedded in a phantom containing various plastic materials [23].

In this paper, we present an extension to the retrieval algorithm developed by Diemoz et al., which enables EI phase retrieval of multi-material samples. Our new algorithm follows the work of Beltran et al. who extended Paganin’s single-image retrieval algorithm for PBI, for the case of a multi-material object [24, 25]. While still requiring only a single image per projection, our retrieval algorithm can correctly reconstruct the interface between any two materials (when one is encased by the other), by tuning the relevant input values. Using a “splicing” method presented here, when a sample is made of various materials, the algorithm can be applied multiple times, once for each pair of materials, and a composite CT slice can then be obtained which will simultaneously display all materials adequately.

2. Theory

In this section we will derive the theory for EI single-image phase retrieval for samples made of multiple materials. A schematic of a laboratory EI system is shown in figure 1. A divergent X-ray beam is divided into beamlets by a mask made of alternating absorbing and transmitting lines, placed immediately before the sample. A second mask is placed before the detector, with a corresponding and magnified period and aperture size. When both masks are misaligned with respect to each other, refraction of the beamlets caused by the sample in the direction perpendicular to the masks’ apertures is converted into detectable intensity changes. These intensity changes are proportional to the first derivative of the phase, ϕ [20]. Image formation in the direction parallel to the masks’ apertures can be treated as PBI, where the signal recorded is proportional to the second derivative of the phase [26, 27]. It should be noted however that in our case, the PBI signal is very weak due to the use of an extended laboratory source. Assuming a near-field regime, and that the sample’s attenuation and phase are varying slowly, the signal can be expressed by the transport-of-intensity equation [26,28].

 

Fig. 1 A top-view schematic of the laboratory implementation of an EI setup. The sample’s rotation axis is aligned with the y-direction.

Download Full Size | PDF

In their work, Diemoz et al. have shown that the normalized signal recorded by an EI setup is given by [21]:

When considering a homogeneous sample which is sufficiently thin as to satisfy the projection approximation [29], the transmission and phase at the sample’s exit surface can be expressed as a function of the sample’s thickness, t(x, y):

Diemoz et al. have shown that by substituting Eqs. (2) and (3) into Eq. (1) and defining J_EI = z₂C′(y_e)C⁻¹(y_e), an expression for the sample’s projected thickness can be obtained:

Therefore, when a homogeneous sample of known material is imaged with monochromatic radiation, the algorithm will retrieve quantitative values for the object’s thickness. While for the case of inhomogeneous samples the retrieved results would not be quantitative, often high image quality can still be achieved by tuning the input δ and μ values following an initial guess. For simplicity, here we will derive the multi-material retrieval algorithm assuming monochromatic radiation.

Here, we adapt the derivation first presented for PBI by Beltran et al. for a ternary object (composed of two materials and voids) [24], to the case of EI. Consider an object made of two materials denoted by 1 & j, where material j is fully encased by material 1. The corresponding refractive indices are n₁ = 1 − δ₁ + iβ₁ and n_j = 1 − δ_j + iβ_j, with β₁ = μ₁/2k and β_j = μ_j/2k. For monochromatic radiation, Eq. (1) describes the detected normalized signal. Assuming the sample satisfies the projection approximation, the transmission and phase shift at the sample’s exit plane are given by:

The third term on the RHS is expanded as:

By inserting these terms into Eq. (1) and rearranging, we get:

By taking the two-dimensional Fourier transform of Eq. (13) and making use of the Fourier derivative theorem, we get the following expression:

An expression for t_j can then be obtained by inverse transforming both sides of Eq. (14) and rearranging:

For a given experimental setup, Eq. (15) can be readily implemented to obtain quantitative results of the projected thickness of material j, provided that the total projected thickness, A, and both materials’ refractive indices are known. The total projected thickness can be easily estimated in common situations where the sample is placed in a cylindrical, full container (as is the case in both examples provided below). However, in other cases where there are voids to consider, A(x, y) can be estimated by forward-projecting a CT slice reconstructed by using Eq. (4) for material 1, and setting a threshold to separate object from void, as was previously suggested by Beltran et al. [24].

For samples containing multiple materials j = 2, 3, ...N, Eq. (15) can be applied N − 1 times, each time adjusting Δμ and γ such that they correspond to the pair of materials of interest. As long as each material j is fully encased by material 1, the interface between the two materials will be correctly reconstructed. Independent N CT slices can be obtained by applying a conventional filtered-backprojection (FBP) algorithm to sinograms of t₁ and t_j=2,3,...N. However, each of these slices will provide an accurate reconstruction of the corresponding pairs of materials, while other interfaces will be either under or over-retrieved (appearing either as residual fringes or as blurred interfaces). A composite CT slice, with all interfaces sharply reconstructed, can then be obtained by splicing together these CT slices corresponding to all different materials. The splicing procedure consists of extracting from each slice the area containing the pair of materials focused on in the retrieval, and digitally inserting these into one composite slice, after the adjustment of background offsets.

Note that while Eq. (15) was derived for the case of monochromatic radiation, an analogous expression could be developed for polychromatic beams. In a similar way to Diemoz’s work [23], since both δ and β vary with energy, a polychromatic extension would involve expressing the measured detector signal as a weighted sum of all its monochromatic components, taking into account the source spectrum and the detector’s energy response. Therefore, as long as there is no significant beam hardening by the sample, Eq. (15) can also be applied to data acquired with a polychromatic beam, if the input values for δ and β are calculated at their effective energies [30].

3. Quantitative validation

The novelty of the proposed algorithm lies in its ability to enable the imaging of samples containing materials with significantly different refractive index properties. In order to test its validity, a simulation study was performed, followed by the collection of experimental data for direct comparison. As a test sample, a water filled cylinder containing aluminium and low-density polyethylene (LDPE) rods was chosen.

A wave-optics simulation, previously developed to replicate the experimental EI setup [31], was used with a photon energy of E = 17.5 keV (corresponding to the k-alpha line of the molybdenum target in the used laboratory X-ray tube). A numerical phantom was created, which consisted of a vertically-aligned cylinder of water (diameter = 1.85 cm) containing an aluminium rod (diameter = 2 mm) and a LDPE rod (diameter = 3.9 mm). Corresponding refractive index values were taken from the ICRU 44 database [32] and are reported in table 1.

Table 1. Input Refractive Index Values - Phantom, Simulation

View Table | View all tables in this article

Mask parameters were chosen to match those used in the laboratory, where each mask is 150 μm thick and is made of gold. The masks’ periods and aperture sizes were 48 μm and 12 μm, respectively, for the sample-mask, and 62 μm and 15 μm for the detector-mask. The masks’ apertures were aligned with the vertical direction. The detector pixel size was 62 × 62 μm². The source-to-sample distance was z₁ = 1.6 m and the sample-to-detector distance was z₂ = 0.4 m. For simplicity, in the simulation we assumed that there is no gap between the sample-mask and sample, and similarly no gap between the detector-mask and detector; in truth, a small (∼1–2 cm) gap is present which is not expected to have any noticeable effect on the images. An illumination curve was simulated to obtain C(y_e) and C⁻¹(y_e). A CT scan was then simulated, consisting of 2400 views over 360 degrees (i.e. angular step = 0.15 degrees). The relative displacement of the sample-mask with respect to the detector-mask was set as y_e = 10 μm. To increase the spatial sampling rate, at each angular view 4 images were taken at different sub-pixel positions of the sample and were later recombined to form a high-resolution image. This process is often referred to as “dithering” in EI, and is not required unless an increase in spatial resolution over that determined by the detector pixel size is sought [33]. The phase retrieval part comprised of three steps. First, Eq. (4) was applied to the raw projections using the refractive index values of water. The second and third steps consisted of applying Eq. (15) to raw projections, once for the water-aluminium interface, and then for the water-LDPE interface, using the values listed in table 1. In this case, A(x, y) was calculated as the projection of a circle with a diameter corresponding to that of the water cylinder. Three separate sinograms were created, corresponding to the different retrieval steps. CT slices of each of these phase-retrieved sinograms were reconstructed by means of a FBP algorithm.

For experimental data collection, a phantom was created which matched the simulated one (apart from the fact that a hollow plastic cylinder of unknown material, with inner and external diameters of 1.75 cm and 1.85 cm, respectively, was filled with water for practical reasons). As a source, the Rigaku MicroMax 007 HF rotating anode (molybdenum) X-ray tube (Rigaku Corporation, Japan) with a focal spot of approximately 70 μm was used. The source was operated at 25 mA and 40 kVp, with a 30 μm Molybdenum filter. A Pixirad single photon counting, energy-resolving detector [34,35], with a pixel size of 62 × 62 μm² was used and placed 2.07 m away from the source. The sample and detector masks were placed at 1.6 m and 1.97 m downstream of the source, respectively, following studies from the method’s early days in which distances between the source, masks and detector were optimized [14,36]. Both masks were fabricated by electroplating gold strips onto a graphite substrate (Creatv Microtech Inc., Potomac, MD, USA). The masks’ period and aperture sizes were the same as used in the simulation, apart from the period of the detector-mask which was 59 μm, to compensate for the limited magnification arising from the small gap between the detector-mask and detector. An illumination curve scan was performed, followed by a CT scan of the phantom with the same parameters reported for the simulation. Flat-field images (i.e. images in the absence of the sample) were acquired at each angular view, and were later used to normalize the sample images. The exposure time per projection was 4 s.

For phase retrieval, Eqs. (4) and (15) were applied to the normalized data. Here, since the source spectrum is polychromatic and the complex refractive index varies with energy, the refractive index values used in the retrieval were estimated. In general, these can be estimated using the concept of effective energies if the source spectrum, detector energy response and sample composition and geometry are known [30]. In the case discussed in this paper, one can consider that an optimal CT reconstruction has been achieved when two conditions are satisfied. The first one is the maximization of the sharpness of the interface between the pair of materials under question. The second is that the reconstructed material j has a mean value of 1 in the CT slice. In order to simultaneously fulfill these conditions, CT reconstruction is carried out multiple times while varying the refractive index values. The chosen values used in the retrieval algorithm of the experimental data are presented in table 2, along with the effective energies these correspond to for each material. These effective energies are reasonable considering the source spectrum produced by the molybdenum target. The difference between the values obtained for water and LDPE compared to aluminium can be explained by taking into account the beam hardening caused by the higher absorption of the aluminium rod. Phase retrieved sinograms were created, followed by a CT reconstruction using FBP.

Table 2. Input Refractive Index Values - Phantom, Experimental

View Table | View all tables in this article

For a quantitative evaluation of the proposed algorithm, simulated and experimental projected thicknesses of each material were plotted together and are shown in Fig. 2. The plots demonstrate the quantitativeness of the algorithm: all simulated profiles retrieve the nominal thickness values accurately (retrieving diameters of 1.85 cm, 2 mm and 3.9 mm for water, aluminium and LDPE, respectively). As expected, in each case the profiles are quantitative only for the interface between the two materials targeted by the algorithm through selection of the relative input parameters. A good agreement can also be noted between the experimental and simulated profiles for each material of interest. Discrepancies in e.g. the retrieved thickness of water, in the location where the aluminium rod is embedded, arise from the difference in (effective) energies between the monochromatic simulation and the experimental polychromatic spectrum. Another inconsistency in retrieved values can be observed between simulated and experimental profiles of both aluminium and LDPE, just near the edges of the cylinder. This inconsistency arises from the fact that while a water cylinder was implemented in the simulation, the experimental phantom consisted of a plastic cylinder filled with water. Finally, by looking at the retrieved projected thickness of both aluminium and LDPE, it can be seen that the algorithm breaks down at the edges of the outer cylinder; this can be expected since in that case our assumption of a slowly varying encasing material thickness is violated. However, this affects the retrieved values only locally, at the edges of the encasing material, and values at the interface of interest are not polluted.

 

Fig. 2 Retrieved projected thicknesses for water, aluminium and LDPE, using values reported in tables 1 and 2. The red arrows point at the material of focus in each case. The experimental profiles were averaged over 5 pixel rows.

Download Full Size | PDF

In order to present all interfaces correctly in a single CT slice, the three different slices were spliced together. This step was performed using the open source software ImageJ, for both the experimental and simulated slices. While the segmentation part of the splicing procedure could be automated using methods such as Otsu thresholding [37, 38] or by developing expressions for the “bleed-width” in a similar fashion to Beltran’s work [24], for this proof-of-concept study the segmentation was done manually, by visually estimating the blurring bleed-width.

Individual simulated CT slices, each displaying reconstructions of the water-air (a), aluminium-water (b) and LDPE-water (c) interfaces are shown in Fig. 3, along with the simulated and experimental spliced slices. Each slice in panels (a–c) demonstrates a sharp reconstruction of the material pair of interest, while residual phase-contrast fringes and blurred interfaces of other materials are highlighted by arrows. In order to achieve a quantitative spliced slice preserving across the phantom the δ values used in the retrieval, each slice reconstructed from t_j had to be manipulated according to: ${\mathit{slice}}_{t_j}^{*}=\left[{\mathit{slice}}_{t_j}\times ({\delta }_j-{\delta }_1)\right]+{\delta }_1$, before digitally inserting its segmented area of interest into the corresponding region of ${\mathit{slice}}_{t_1}^{*}={\mathit{slice}}_{t_1}\times {\delta }_1$. The same splicing procedure was applied to both simulation and experimental slices, however, the δ values used for quantitative splicing were different. For each case, the δ values applied during splicing were the same as the ones previously used in the retrieval algorithm, hence some differences exist (most notably for aluminium, see tables 1 and 2). This can be appreciated in Fig. 3(f) where profiles across both spliced slices are plotted together for comparison. Ignoring the discrepancies which have already been discussed above (due to the plastic container and the polychromaticity of the spectrum), both profiles demonstrate a sharp, quantitative reconstruction of all parts of the phantom. These results confirm that the algorithm is indeed quantitative for the case of monochromatic radiation, and can be applied when a polychromatic spectrum is used, if the concept of effective energies is adopted [30].

 

Fig. 3 Simulated CT slices focusing on the water-air (a), aluminium-water (b) and LDPE-water interfaces, using E=17.5 keV. The white arrows indicate artefacts arising from locally incorrect choices of refractive index values. Simulated (d) and experimental (e) spliced slices along with a plot of profiles through them (f), in the position indicated by the red line in (d). The difference in the retrieved values of the aluminium rod in (f, see arrow) is expected and is due to the difference in effective energies (see tables 1 and 2).

Download Full Size | PDF

4. Qualitative evaluation

Following validation by a phantom study, the proposed retrieval algorithm and splicing method were applied to experimental data of a complex biological sample. The imaged sample was a chicken bone surrounded by soft tissue, placed in a plastic cylinder of approximately 8 mm diameter. The sample was “fresh”, i.e. no additional sample preparation was applied (e.g. formalin fixation). The same source used for the phantom scan was used, operated at 40 kVp and 25 mA, with no additional filters. The detector used for this scan was a Hamamatsu C9732DK flat panel detector with a passive-pixel CMOS sensor (Hamamatsu, Japan), and was positioned 2 m away from the source. The detector pixel size was 50 × 50 μm², however the effective pixel size in the x-direction was 100 μm due to the line-skipping design of the detector-mask, where every second detector pixel column is completely covered by the detector-mask. This mask design is used to reduce the negative effect of pixel cross-talk on the EI signal [39]. The sample and detector masks were placed at 1.6 m and 1.96 m downstream the source, respectively. The sample-mask period and aperture size were 79 μm and 10 μm, respectively, while the detector-mask period was 98 μm with a 17 μm aperture size. Following an illumination curve scan, a CT scan was performed with 720 views over 360 degrees (i.e. 0.5 degrees step). To increase the spatial resolution, 8 projections were taken at different sub-pixel positions of the sample, at each angular view. The exposure time per projection was 3 s and the relative position of the sample-mask with respect to the detector-mask was set as y_e = 8 μm. Raw projections were corrected for dark current and flat-field inhomogeneities prior to the phase retrieval step, which used Eq. (4) once (for the cylinder-air interface) and Eq. (15) twice (where the interfaces between soft tissue types, e.g. fat-muscle, and between soft tissue and bone were considered). For the retrieval of the bone-soft tissue interface, refractive index values of bone and water (commonly used as a tissue-equivalent material) at an estimated effective energy of 22 keV were used [32]. Here, the effective energy was estimated according to the theoretical source spectrum and a predicted linear detector energy response, since the detector used in this setup was an integrating detector with an energy response function which has yet to be fully characterized (unlike in the phantom scan reported in section 3, where a photon-counting detector was employed). Moreover, for the retrieval of the interface between soft tissue types, since the specific tissue types were not known, and the complex sample geometry in this case did not strictly satisfy the algorithm’s condition for a material fully encased by another, a purely qualitative evaluation was performed. Therefore, to maximize the contrast between the different soft tissues, an iterative approach was employed, where Δμ and γ were varied until the contrast was maximized. The reconstructed slices are shown in Fig. 4. On the left hand side, the three separate reconstructions are shown, each rendering a specific interface while exhibiting the expected blurring\fringes at other interfaces. A spliced, composite CT slice is presented in the larger image on the right, in which all different materials are adequately reconstructed. As previously noted, although the algorithm was derived for the case of a material j fully encased within material 1, this is not strictly the case for the chicken bone sample. However, as already observed by Beltran, a locally incorrect choice of refractive index values will only affect the part of the interface for which the choice is not ideal [24]. Since quantitative results were not expected, the composite slice was created by windowing each retrieved slice differently, in order to achieve the best visualization. As before, the splicing procedure was performed using ImageJ and manual segmentation. The spliced image shown in Fig. 4 demonstrates that this method could be used to obtain high-quality, albeit not quantitative, images of complex samples with unknown composition, using a single image per projection.

 

Fig. 4 Axial slices of a chicken bone sample reconstructed using values optimized for the cylinder-air interface (a), intra-soft tissue contrast (b) and bone-soft tissue interface (c). Panel (d) shows a composite slice obtained by splicing slices (a–c).

Download Full Size | PDF

5. Discussion and conclusion

We have derived an algorithm capable of quantitatively retrieving the projected thickness of multi-material samples illuminated by monochromatic radiation, using a single input image acquired with an EI setup. A CT phantom simulation study validated the algorithm’s accuracy, while results from an experimental scan provided good agreement, within the limits of polychromaticity. A “splicing” method was then presented as means for constructing a composite CT slice where all parts of a sample are adequately reconstructed, thereby eliminating artefacts previously arising from locally incorrect choices of refractive index values for multi-material samples, during the phase retrieval step. For strictly quantitative results, one material must be fully encased by the other, and the refractive indices of both materials and the total projected thickness must be known. Therefore, for quantitative imaging, certain applications (e.g. in the field of materials science) could benefit from using the proposed method on data acquired using monochromatic synchrotron radiation. Although we have yet to apply the algorithm to synchrotron radiation data, we expect it to provide more accurate quantitative CT reconstructions, as was demonstrated via simulation. However, we have shown that enhanced image quality and visualization can be achieved for complex samples of unknown materials, imaged in a laboratory environment.

It should be noted that, if the requirement for quantitativeness is relinquished and when splicing is used, similar results can be obtained by using Eq. (4) multiple times, and replacing δ and μ with Δδ = δ_j − δ₁ and Δμ = μ_j − μ₁. This approach might be more practical since it does not require the total projected thickness as input; however, the retrieved projected thickness profiles will not be quantitative as they are not “normalized” by the encasing material. During the splicing step, the CT slices are scaled and can be windowed to achieve a similar result to the ones obtained by using Eq. (15).

The EI single-image retrieval algorithm was recently developed and has shown a huge potential for reducing scan times, particularly for CT scans, by simplifying the acquisition sequence [23]. Previously, the position of the sample-mask had to be varied at least once per projection, leading to increased scan times with a potential for introducing instabilities to the imaging system. In this proof-of-principle study, dithering was performed in both reported scans, meaning that multiple images of the sample at different sub-pixel positions were acquired at each view, thereby increasing the delivered dose and scan duration. However, dithering is not strictly necessary and is only used to increase the spatial resolution, as previous studies have shown that quantitative information can be extracted from EI CT when no dithering is performed [33]. Although not exploited in this study, when dithering is not required, the single-image retrieval algorithm enables a continuous rotation of the sample, which can in turn lead to significantly reduced scan times, as was recently reported [23]. We therefore expect that the development presented in this paper will enable obtaining similar results for complex samples containing materials with a wide range of refractive index values. Furthermore, since the algorithm requires one input image as opposed to the previously needed minimum of two, it will be interesting to investigate its potential for dose reduction strategies.