X-ray phase-contrast radiography and tomography with a multi-aperture analyser

TL;DR

This paper introduces a multi-aperture analyser setup for X-ray phase-contrast imaging and tomography, enabling quantitative 3D reconstructions and simultaneous retrieval of absorption, phase, and scattering with high sensitivity.

Contribution

The work presents a novel multi-aperture analyser configuration that allows for quantitative phase-contrast imaging and tomography, including simultaneous measurement of absorption, phase, and scattering.

Findings

Successful tomographic reconstruction using the method.

Reliable retrieval of absorption and phase even with strong scatterers.

High sensitivity and dynamic range demonstrated on custom phantoms.

Abstract

We present a multi-aperture analyser set-up for performing X-ray phase contrast imaging in planar and three-dimensional modalities. A multi-slice representation of the sample is used to establish a quantitative relation between projection images and the corresponding three-dimensional distributions, leading to successful tomographic reconstruction. Sample absorption, phase and scattering are retrieved from the measurement of five intensity projections. The method is experimentally tested on custom built phantoms with synchrotron radiation: sample absorption and phase can be reliably retrieved also in combination with strong scatterers, simultaneously offering high sensitivity and dynamic range.

Equations

$$I(\bar{y})$$

$$=$$

$$t$$

$$=$$

$$\Delta y_R$$

$$=$$

$$t$$

$$\Delta y_R$$

$${\sigma }_O^2$$

Full text

Conventional radiography is based on the attenuation of X-rays traversing a sample for generating contrast. In X-ray phase contrast imaging (XPCI) Bravin et al. (2013); Wilkins et al. (2014) additional contrast mechanisms, generating from the phase-shifts imparted to the beam by the sample, contribute to the image formation process. This can improve the visibility of a large variety of details of interest and its applications span across many different fields, encompassing material science, security screening, biology and medicine. Amongst various implementations with which it is possible to obtain X-ray phase-contrast images Bonse and Hart (1965); Goetz et al. (1979); Davis et al. (1995); Ingal and Beliaevskaya (1995); Wilkins et al. (1996); Chapman et al. (1997); Clauser (1998); David et al. (2002); Momose et al. (2003); Mayo and Sexton (2004); Pfeiffer et al. (2006); Wen et al. (2009); Morgan et al. (2012); Miao et al. (2016), we focus here on edge illumination (EI) Olivo et al. (2001). EI can provide quantitative attenuation, phase Munro et al. (2012) and dark-field Endrizzi et al. (2014a) representation of a sample; and can be adapted for use with synchrotron radiation, microfocus tubes and conventional rotating anode sources with extended focal spots Olivo and Speller (2007); Endrizzi et al. (2014b). It exhibits negligible requirements in terms of temporal or spatial coherence Munro et al. (2010); Endrizzi et al. (2015a), provides high sensitivity Marenzana et al. (2012); Diemoz et al. (2013a, b), it is robust against mechanical and thermal instabilities Millard et al. (2013); Endrizzi et al. (2015b) and enables low-dose implementations of XPCI in planar and three-dimensional imaging Olivo et al. (2013); Hagen et al. (2014).

We introduce here a new multi-aperture analyser set-up for performing edge illumination with synchrotron radiation in situations where very wide angular ranges must be explored, for example in the presence of strong scatterers. Through experimental tests on custom built phantoms, we show that refraction can be accurately retrieved independently from the presence of large amount of scattering. Two spatial resolutions are investigated, showing that the retrieved signals are independent from this parameter. A simple model based on a multi-slice representation of the sample enables the reconstruction of the three-dimensional images. The method retrieves absorption, refraction and dark-field images, simultaneously providing high sensitivity and large dynamic range for both planar and three dimensional XPCI applications.

The experimental set-up is sketched in Figure 1(a).

The beam from a bending magnet source is monochromatised with a double Si (111) crystal and then shaped by a narrow aperture, before it impinges on the sample. It is then analysed by a set of apertures positioned immediately before a digital detector. In order to achieve different levels of illumination, the detector and the detector apertures are scanned along the $y$ direction. The data collected with each detector exposure result in a single line of pixels, and a two-dimensional image is built by scanning the sample along the $y$ direction and exposing the detector multiple times.

The illumination function describes how the detected intensity varies depending on the relative displacement $\bar{y}$ between pre-sample and detector apertures. By using a geometrical optics model, the intensity measured at the detector is described as a convolution between the illumination function $L(\bar{y})$ and the object function $O(x,y)$:

[TABLE]

Transmission, refraction and scattering properties of the sample were retrieved by using a multi-Gaussian model for the illumination function Endrizzi et al. (2014a):

[TABLE]

where $L(\bar{y})=$ $\sum_{n}(A_{n}/\sqrt{2\pi\sigma_{n}^{2}})$ $\exp{[-{(\bar{y}-\mu_{n})}^{2}/2\sigma_{n}^{2}]}$, ($n=1\ldots N$); and $O(\bar{y})=$ $\sum_{m}(A_{m}/\sqrt{2\pi\sigma_{m}^{2}})$ $\exp{[-{(\bar{y}-\mu_{m})}^{2}/2\sigma_{m}^{2}]}$, ($m=1\ldots M$). The parameters are defined as follows: $\mu_{mn}=\mu_{m}+\mu_{n}$, $\sigma_{mn}^{2}=\sigma_{m}^{2}+\sigma_{n}^{2}$ and $A_{mn}=A_{m}A_{n}(1/\sqrt{2\pi\sigma_{mn}^{2}})$. We note that beam absorption and refraction can be both included in the object function $O(y)$ as a multiplicative factor $t$ and a shift $\Delta y_{R}$ of the center of the distribution, respectively. In this formulation, a purely absorbing object is represented by a Dirac’s delta function centred in zero and multiplied by a factor $t$, which indicates the transmitted to incident intensity ratio. If the sample is also refracting, the delta function is laterally shifted $O(y)=t\delta(y-\Delta y_{R})$, thus the effect on the illumination function is a reduction in the total transmitted intensity by the factor $t$ plus a $\Delta y_{R}$ lateral shift due to refraction. In the case of a sample exhibiting also dark-field contrast, the object function $O(y)$ has a finite width and the presence of the sample results into a broadened illumination function.

Phase images can be obtained by numerical integration of the refraction images by taking into account that the refraction angle $\alpha=(\lambda/2\pi)\partial_{y}\Phi(x,y)$ is directly proportional to the gradient of the object’s phase shift $\Phi(x,y)$, where $\lambda$ is the wavelength. In the case of a parallel beam geometry, of which a synchrotron set-up is usually a good approximation, the refraction angle can be measured as the ratio between the relative shift of the illumination function and the sample-to-detector distance $\alpha=\Delta y_{R}/z_{sd}$.

A thick sample can be represented with a multi-slice approach, by separately considering the subsequent effect of thin sample sections along the beam axis $z$. We define $O_{k}$ a thin but finite thickness $\Delta_{z}$, section of the object along $z$. The entire object is obtained by the sum of its sections (with $k=1\ldots K$) and its extent along $z$ is given by $Z_{o}=\Delta_{z}K$. The intensity measured at the detector can then be expressed in terms of the contributions from separate layers

[TABLE]

where the more compact notation $*$ was introduced to indicate convolution. Let’s consider, for simplicity’s sake, a single-Gaussian model for the object’s section functions $O_{k}$. In this case, the entire object function $O(y)=(O_{1}*O_{2}*\ldots*O_{K})(y)$ is still a Gaussian, with mean $\mu_{O}=\sum_{1}^{K}\mu_{k}$, variance $\sigma_{O}^{2}=\sum_{1}^{K}\sigma_{k}^{2}$ and amplitude $A_{O}=(1/\sqrt{2\pi\sigma_{O}^{2}})\prod_{1}^{K}A_{k}$; this is now expressed as a function of the contributions from the individual thin layers. The intensity transmitted by the entire object becomes

[TABLE]

and a similar result is obtained for the refraction

[TABLE]

where $\delta$ is the decrement from unity of the material’s refractive index, $\beta$ its imaginary part and $\Phi_{k}=-(2\pi/\lambda)\delta_{k}\Delta_{z}$ is the phase shift introduced by the $k-$th object’s layer. If the section thickness $\Delta y$ is small enough to allow a transition to the continuous formulation, we obtain

[TABLE]

and similarly for the dark-field signal

[TABLE]

Along with absorption and phase, also the dark-field signal $\sigma_{O}^{2}$ can be cast as an integral along the beam path. By collecting a number of different views while rotating the sample around the $y$ axis, the three dimensional distribution of the width of the object function $\sigma^{2}(x,y,z)$ can be calculated by means of, for example, the Inverse Radon Transform or the Filtered Back Projection algorithms as is routinely done for absorption and phase-based computed tomography. Equivalent results were obtained in the context of analyser based imaging Rigon et al. (2008), grating interferometry Bech et al. (2010) and beam-tracking Vittoria et al. (2015).

The experiment was performed at the SYRMEP beamline (Elettra Sincrotrone Trieste, Italy). The pre-sample slit was $40$ mm along $x$ and its aperture along $y$ was set at $20$ $\mu$m. The detector apertures were $23$ $\mu$m and arranged at a regular period of $79$ $\mu$m. The beam energy was $20$ keV, the pre-sample aperture was at about $22$ m from the source, the sample stage was at $26$ cm from the aperture and the analiser and detector a further $2.5$ m downstream. The detector was a CCD with GadOx scintillation screen (Photonic Science, UK) with a pixel size of 12.5 $\mu$m.

For this experiment, five terms were retained for the illumination function (N=5) and the scattering distribution was assumed to be a single Gaussian (M=1). The three parameters representing absorption, refraction and scattering in the sample were obtained by a pixel-wise non-linear fitting procedure Endrizzi et al. (2015b) that compared the intensities recorded with and without the sample in the beam. For image acquisition, the illumination function was sampled in five positions with $\{\pm 24,\pm 12,0\}$ $\mu$m displacement with respect to the position of maximum intensity in the central slit of the analyser, with an exposure time of $400$ milliseconds each.

The phantom used for planar imaging was composed of a cylinder of acrylic material with density $1.18$ g/cm3 and radius $1.45$ mm and a step wedge made of paper layers. A melamine sponge prism was used to experimentally measure the dependency of the dark-field signal upon thickness. Mono disperse borosilicate micro spheres with diameters of $5$, $10$ and $12$ $\mu$m were embedded into an acrylic support for computed tomography acquisitions. A plastic scaffold was used to test the three-dimensional reconstruction on a phantom with a more complex geometry.

The data were recorded at a pixel size of $12.5$ $\mu$m and subsequently binned in such a way that the integrated intensity going through a single aperture in the detector mask was combined in a single image line. For the radiography image mode, the sample was scanned along $y$ in $15$ $\mu$m steps to build up a two dimensional image. For the computed tomography image mode, $600$ views were acquired with $0.3$ degree angular step.

The illumination function, acquired with a fine sampling step, is compared to the one obtained through a wave optics numerical simulation Vittoria et al. (2013) that incorporates all the experimental parameters. A good match between the intensities experimentally recorded in the positions corresponding to the five apertures at the detector mask and the ones expected from the simulation can be observed in Figure 1(b).

This illumination function is taken as the reference one, and subsequently used for the retrieval of the sample images. The approximate rocking curve of a Si (111) crystal at the energy of $20$ keV is also plotted in Figure 1(b), for a visual comparison of the sensitivity curves of the EI and the analyser crystal based imaging methods. It is interesting to note that the EI configuration reported here offers comparable curve width to that of the crystal analyser technique, at least when the $(111)$ diffraction is exploited in the latter.

In order to investigate the dependence of the retrieved signals on the spatial resolution, data were analysed two times. In the first case (high resolution) the pixel size was $50$ $\mu$m $\times$ $60$ $\mu$m while in the second case (low resolution) the pixel size was $500$ $\mu$m $\times$ $150$ $\mu$m. These two sets of data were obtained by binning the raw intensity acquisitions, after which the same analysis procedure is applied to both datasets. The results of this analysis are summarised in Figure 2. The effect of lowering the spatial resolution can be observed by comparing panels (a)-(c) to panels (d)-(f) in Figure 2, where the high- and low-resolution images are reported, respectively. For a quantitative comparison, line profiles are shown in the (g) and (h) panels of Figure 2. They were extracted along the directions highlighted in the panels (c), (f) for the plots in (g) and in the panels (b), (e) for the plots in (h). They represent the dark-field signal measured along the paper step wedge and the phase signal measured across the plastic cylinder. As it can be seen, the values obtained in the high- and low-resolution configurations are quantitatively very close to each other. Moreover, the theoretical phase shift shows a very good match with the experimentally measured one ($\delta=6.61\cdot 10^{-7}$ and $\beta=3.34\times 10^{-10}$ CSIRO CSS eToolbox ), even though the experimental profiles in Figure 2(h) were extracted from the portion of the cylinder placed behind the part of the wedge with strongest scattering. The dispersion of the retrieved refraction angles, measured in an empty background region (around zero), can be used to estimate the sensitivity of the imaging system Modregger et al. (2011); Diemoz et al. (2013a, b). By following this procedure, the high and low spatial resolution configurations gave a standard deviation of $17$ nrad and $5.9$ nrad, respectively. The behaviour of transmission and dark-field signal as a function of thickness were investigated by means of the melamine sponge prism. Two line plots are reported in panels (a) and (b) of Figure 3. It can be clearly seen that the dark-field signal grows linearly with the sample thickness, and that the expected exponential relationship holds for the transmission. The dark-field CT result for the micro spheres sample is shown in the inset of Figure 3(b), along with the standard absorption CT image again showed as an inset in Figure 3(a). The values $\{0.95\pm 0.06$, $1.10\pm 0.05$, $1.19\pm 0.08\}$ $\cdot 10^{3}$ $\mu$rad2 mm*-1* were measured for the $\{5$, $10$, $12\}$ $\mu$m diameter spheres details, respectively.

Finally, as an example for a phantom with a more complex geometry, the image obtained by fusing absorption- and dark-field- CT reconstructions of a plastic scaffold is shown in Figure 4. The two colour channels emerge from different sample details, meaning that they represent different properties of the sample and can offer a complementary visualisation of its characteristics.

In conclusion, we have presented a multi-aperture analiser set-up for radiography and computed tomography X-ray phase-contrast imaging. The linearity of the dark-field signal with the sample thickness was experimentally verified and quantitative computed tomography was performed on a custom phantom made of micro spheres embedded into an acrylic support. The method is based on the edge illumination method and it simultaneously provides absorption, refraction and dark-field contrast images. Very good agreement was observed between theory and experiment for the phase-shift image of a test sample. Notably, the phase shift was quantitatively measured to a high degree of accuracy also in the presence of strong scatterers in the beam. Finally, as an example of a complex geometry sample, the images of a scaffold were presented.

Acknowledgements

This project was supported by the UK Engineering and Physical Sciences Research Council Grant EP/I021884/1. ME was supported by the Royal Academy of Engineering under the RAEng Research Fellowships scheme. We thank Elettra Sincrotrone Trieste for access to SYRMEP beamline (proposal 20140147) that contributed to the results presented here.