Introduction

Several imaging modalities are used for the detection and monitoring of abnormalities within the human brain, most commonly Magnetic Resonance Imaging (MRI) and Computed Tomography (CT). While MRI is the preferred clinical modality for imaging the anatomy of the brain, due to its lack of ionising radiation and a higher contrast resolution1, CT plays an essential role in diagnostic imaging due to a much faster acquisition speed, which reduces movement artefacts, and a higher spatial resolution. In addition, CT can be used in the presence of metallic implants and tattoos with metallic ink, for which MRI is contraindicated due to the dangers imposed by the high magnetic fields required. CT is also better suited for the evaluation of penetrating brain injuries and other trauma and acute neurological emergencies2,3, where impromptu access to MRI scanners for time-critical assessment is often unavailable, and the presence of magnetic materials cannot be determined prior to the scan.

X-ray absorption has provided the contrast for computed tomography (CT) since the inception of CT in the 1960s and early 1970s4,5,6, enabling the reconstruction of internal ‘slices’ of the human body and of other objects from a series of external views. While absorption CT is very well-suited to imaging soft tissue and bone, its use for delineating more subtle features within the tissues is more limited. Phase contrast X-ray imaging (PCXI) allows us to exploit the diffraction of X-rays, rather than their attenuation, resolving features in soft tissues and other low-density, low-Z materials for which standard attenuation imaging is insufficient. The refractive index decrement, which describes the phase change of X-rays as they pass through matter, is typically three orders of magnitude smaller than the attenuation coefficient for the case of soft tissue in the diagnostic X-ray regime7,8; this suggests that a significant increase in contrast resolution can potentially be obtained with PCXI. It has recently been shown that phase contrast can enable a radiation dose reduction in CT by factors in the thousands over conventional CT without loss in image quality9.

The highest spatial resolution currently achieved for ex vivo MRI imaging of a brain is around 50–100 μm when using extremely high field strengths on the order of 9–16 T10 (for reference, clinical CT scanners operate up to 3 T). The best in vivo MRI resolution achieved to date is 100 μm using 7.00–11.7 T scanners11,12. To achieve these resolutions, acquisition times must be quite long, typically 1–2 hours. For synchrotron-based phase contrast computed tomography, the highest spatial resolution achieved to date is ~1 μm13. While prior work has visualised the brain in mouse fetuses14, to date no in vivo or in situ phase contrast X-ray CT (PCXI-CT) of a brain within a substantially calcified skull has been published, which is likely due to the significant artefacts that typically arise from the skull upon back-projection.

Recent studies have demonstrated that PCXI-CT is an effective tool for small animal studies of the brain, providing high resolution images of tissue structures and clear delineation between grey and white matter15,16,17,18,19. Beltran et al.17,20. showed a 200-fold increase in signal-to-noise ratio using PCXI-CT over absorption contrast, indicating that in situ PCXI-CT can lead to very large improvements over conventional absorption contrast CT. A similar example from our data collected at the SPring-8 synchrotron in Japan is shown in Fig. 1, showing the signal-to-noise ratio (SNR) and contrast-to-noise ratio (CNR) improvement that can be achieved with phase contrast and phase retrieval.

Figure 1
figure 1

(a) A single propagation-based projection image of an excised close-to-term rabbit kitten brain, suspended in agarose. The bracket indicates the position of the brain in the tube. (b) Phase contrast tomographic slice of the rabbit kitten brain in (a). Note that the visible structures can only be seen due to the phase shifts imparted by propagation and would not be visible at all in absorption contrast CT. (c) Phase-retrieved tomogram of the brain from (a) and (b). CTs were acquired at 24 keV with an object-to-detector distance of 5 m. See methods section for experimental details.

These previous studies have all been limited to ex vivo imaging on brains that have been excised from their skulls. While these results show the clear potential for brain PCXI-CT in preclinical studies, in vivo imaging is much more difficult due to the strong phase gradients between tissue and bone as well as the strong absorption by bone, causing artefacts from the skull which obscure structures that would otherwise be well-resolved. Overcoming these artefacts is important for future in vivo preclinical research using imaging and, ultimately, for adaptation to the clinic. Herein we demonstrate the first visualisation of the brain in situ in a small animal model, performed using propagation-based PCXI-CT.

Background

Propagation-based phase contrast imaging

Propagation-based imaging (PBI) is the simplest phase-contrast method, wherein the phase shifts resulting from refraction within the object are converted to intensity variations via propagation between the object and the detector, with no optical components required along the path. As the wavefront propagates, small differences in phase accumulate between contrasting materials and/or changing thicknesses so that Fresnel fringes become clearly visible at the detector. The experimental setup for PBI differs from conventional X-ray absorption imaging only in the distance between the imaged object and the detector and in the requirement of a source with sufficiently high spatial coherence21,22.

As incident X-rays pass through a sample, the intensity measured downstream at the detector contains a combination of attenuation and phase information. In Fig. 2, a simple, single-material sample is shown with the resultant image at the detector exhibiting contrast that is mostly due to attenuation in the regions within and surrounding the sample; however the phase effects are clear at the boundaries between materials of differing projected refractive index (in this case, the sample material and the surrounding air). This results in the fringe pattern discussed above, seen on the right side of the figure. This fringe pattern can be used to extract quantitative information about the sample, as described in the next section.

Figure 2
figure 2

Propagation-based X-ray phase contrast imaging (adapted from Kitchen et al.40).

Simplified phase retrieval

The behaviour of an X-ray wavefield, as it passes through an object, is governed by the complex refractive index of the sample, which, for a monochromatic incident wavefield, ignoring polarization, is given by,

$$n({\bf{r}})=1-\delta ({\bf{r}})+i\beta ({\bf{r}}\mathrm{).}$$
(1)

At each point r in the sample, the real part δ is known as the refractive index decrement, describing the phase component, and β = λμ/4π is the attenuation index, describing the absorption component, where λ is the wavelength and μ is the attenuation coefficient of the material.

Under the projection approximation, which assumes that the transverse scattering contribution to the deflection of X-rays through a sample is negligible, the intensity immediately downstream of the sample (i.e. the detector side of the sample) is given by Beer’s Law,

$$I({{\bf{r}}}_{\perp })={I}_{0}{e}^{-\int \mu ({r}_{\perp },z)dz},$$
(2)

where I0 is the intensity of the monochromatic plane waves assumed to be incident on the sample, z is the direction of propagation, and r represents the coordinates perpendicular to propagation. Similarly, the phase is given by,

$$\varphi ({{\bf{r}}}_{\perp })=-\,k\int \delta ({{\bf{r}}}_{\perp },z)dz.$$
(3)

For a single-material object, these become,

$$I({{\bf{r}}}_{\perp })={I}_{0}{e}^{-\mu T({{\bf{r}}}_{\perp })}\,{\rm{and}}\,\varphi ({{\bf{r}}}_{\perp })=-\,k\delta T({{\bf{r}}}_{\perp }),$$
(4)

where T(r) is the projected thickness of the sample. Since the intensity of a propagating wave is I(r) = |ψ(r)|2, the propagating wavefield ψ can be represented as:

$$\psi ({\bf{r}})=\sqrt{I({\bf{r}})}{e}^{i\varphi ({\bf{r}})}.$$
(5)

In the near-field regime (Fresnel number \({N}_{F}\gg 1\)), the transport-of-intensity equation23 holds:

$${\nabla }_{\perp }\cdot [I({\bf{r}}){\nabla }_{\perp }\varphi ({\bf{r}})]=-\,k\frac{\partial I({\bf{r}})}{\partial z}.$$
(6)

Here, denotes the gradient operator in the xy plane perpendicular to the optic axis z.

From this equation and from the intensity and phase given above under the projection approximation, Paganin et al.24 derived the following expression to recover the projected thickness for a single-material sample,

$$T({{\bf{r}}}_{\perp })=-\,\frac{1}{\mu }{\mathrm{log}}_{e}({ {\mathcal F} }^{-1}\{\mu \frac{ {\mathcal F} \{I({{\bf{r}}}_{\perp },z={R}_{2})/{I}_{0}\}}{{R}_{2}\delta {|{{\bf{k}}}_{\perp }|}^{2}+\mu }\}),$$
(7)

where I(r, z = R2) is the intensity at the detector plane and R2 is the sample-to-detector distance. Equation (7) was expanded upon by Beltran et al.20 for the case of a two-material sample, giving:

$${T}_{2}({{\bf{r}}}_{\perp })=-\,\frac{1}{{\mu }_{2}-{\mu }_{1}}\times {\mathrm{log}}_{e}({ {\mathcal F} }^{-1}\{\frac{1}{[{R}_{2}({\delta }_{2}-{\delta }_{1})/({\mu }_{2}-{\mu }_{1})]{{\bf{k}}}_{\perp }^{2}+1} {\mathcal F} \,[\frac{{\rm{I}}({{\bf{r}}}_{\perp },z={R}_{2})}{{I}_{0}{e}^{-{\mu }_{1}A({{\bf{r}}}_{\perp })}}]\}),$$
(8)

where the subscripts 1 and 2 refer to the two different materials, the first embedded within the second, and A(r) is the combined projected thickness, A(r) = T1(r) + T2(r). Determining A(r) is not easy, or may only be possible with limited accuracy, depending on the sample. Instead, we diverge slightly from their method, as outlined below.

Recovery of the exit surface intensity of an object from projection images of one material embedded within another can be achieved by modifying the work of Beltran et al.20. Their Equation (8) is

$$\frac{{\rm{I}}({{\bf{r}}}_{\perp },z={R}_{2})}{{I}_{0}\exp [-{\mu }_{1}A({{\bf{r}}}_{\perp })]}=[-\frac{{R}_{2}({\delta }_{2}-{\delta }_{1}){\nabla }_{\perp }^{2}}{({\mu }_{2}-{\mu }_{1})}+1]{\rm{e}}{\rm{x}}{\rm{p}}[-({\mu }_{2}-{\mu }_{1}){T}_{2}({{\bf{r}}}_{\perp })]\mathrm{.}$$
(9)

Here

$$A({{\bf{r}}}_{\perp })={T}_{1}({{\bf{r}}}_{\perp })+{T}_{2}({{\bf{r}}}_{\perp }\mathrm{).}$$
(10)

Since A(r) is assumed to be a slowly varying function of r, the denominator of the left-hand side of Equation (9) can be approximated as a constant when a spatial derivative such as \({\nabla }_{\perp }^{2}\) is applied to it. When both sides of Equation (9) are multiplied by this term, it can therefore be grouped with the exponential term on the right-hand side to give

$${I}_{0}\exp [-{\mu }_{1}A({{\bf{r}}}_{\perp })]\exp [-({\mu }_{2}-{\mu }_{1}){T}_{2}({{\bf{r}}}_{\perp })]={I}_{0}\exp [-{\mu }_{1}{T}_{1}({{\bf{r}}}_{\perp })-{\mu }_{2}{T}_{2}({{\bf{r}}}_{\perp })]={I}_{0}{\rm{I}}({{\bf{r}}}_{\perp },z=\mathrm{0).}$$
(11)

Therefore, we can rewrite Equation (9) as

$$\frac{{\rm{I}}({{\bf{r}}}_{\perp },z={R}_{2})}{{I}_{0}}=[-\frac{{R}_{2}({\delta }_{2}-{\delta }_{1}){\nabla }_{\perp }^{2}}{({\mu }_{2}-{\mu }_{1})}+1]{\rm{I}}({{\bf{r}}}_{\perp },z=\mathrm{0).}$$
(12)

An alternative approach is to consider that, for a two-material sample under the projection approximation, Equations (4) become,

$$I({r}_{\perp })={I}_{0}{\rm{e}}{\rm{x}}{\rm{p}}[-{\mu }_{1}{T}_{1}({{\rm{r}}}_{\perp })-{\mu }_{2}{T}_{2}({{\rm{r}}}_{\perp })]\,{\rm{and}}\,\varphi ({{\bf{r}}}_{\perp })=-\,k\delta {T}_{1}({{\bf{r}}}_{\perp })-k\delta {T}_{2}({{\bf{r}}}_{\perp }\mathrm{).}$$
(13)

Substituting for T1(r) using Equation (10) gives,

$$I({{\bf{r}}}_{\perp })={I}_{0}\exp [-{\mu }_{1}A({{\bf{r}}}_{\perp })]exp[-({\mu }_{2}-{\mu }_{1}){T}_{2}({{\bf{r}}}_{\perp })]\,\,{\rm{and}}\,\,\varphi ({{\bf{r}}}_{\perp })=-\,k{\delta }_{1}A({{\bf{r}}}_{\perp })-({\delta }_{2}-{\delta }_{1}){T}_{2}({{\bf{r}}}_{\perp }\mathrm{).}$$
(14)

Assuming, again, that A(r) can be approximated as a constant, this results in a multiplicative change in the measured intensity and an additive shift in the phase, which disappears on differentiation in Equation (6); hence, this becomes like a single-material object where δ = δ2 − δ1 and μ = μ2 − μ1 in Equation (7).

Following Paganin et al.24, Beltran et al.20 and Gureyev et al.25, the “absorption” contrast image can be recovered using the Fourier derivative theorem:

$$I({{\rm{r}}}_{\perp },z=0)={ {\mathcal F} }^{-1}\{\frac{1}{[{R}_{2}({\delta }_{2}-{\delta }_{1})/({\mu }_{2}-{\mu }_{1})]{{\bf{k}}}_{\perp }^{2}+1} {\mathcal F} [\frac{{\rm{I}}({{\bf{r}}}_{\perp },z={R}_{2})}{{I}_{0}}]\}.$$
(15)

A tomographic tilt series of such projections can then be used to quantitatively reconstruct the objects embedded within the encasing material. A similar approach can be used to solve for the encasing material itself by following the algorithm of Paganin et al. We see that our revised equation does not require knowledge of the total projected thickness A(r), unlike the algorithm of Beltran et al. (Equation (8)). We also note that, for tomography, we need only recover the absorption contrast image before using a tomographic reconstruction algorithm. In the slightly different context of phase retrieval of several sharp boundaries in tomography, Haggmark et al.26 recently came to the same equation, also utilising the assumption of a slowly-varying thickness in addition to assuming a linear relationship between δ and μ for different materials at a given energy.

Single image phase retrieval of a brain in situ

For materials of similar composition (and hence refractive index), such as the grey and white matter of the brain, the single-material phase retrieval algorithm applied with respect to the tissue/air interface works quite well to resolve the boundaries between these materials. Structures that are unresolved or poorly resolved with attenuation contrast alone become visible upon phase retrieval, since the phase retrieval filter suppresses noise, thereby increasing the SNR and CNR. Ideally, phase retrieval performed with respect to the grey/white matter boundary provides the best contrast resolution of brain structures, however when imaging the brain in situ, the inaccurate phase retrieval at the bone/tissue interface causes excessive blurring that overwhelms the features within the soft tissue of the brain. Figure 3 shows the in situ analogue to Fig. 1. In panel 3c, phase retrieval has been performed with respect to bone/air interface, resulting in an over-blurring at both the grey/white matter and bone/tissue boundaries. Nevertheless, the brain structure is more clearly resolved in Fig. 3c than in Fig. 3b, despite the incorrect phase retrieval at the bone/tissue interface and the associated artefacts caused by the the highly attenuating bone. In panel 3d, phase retrieval has been performed with respect to the tissue/air interface, again resulting in a highly-blurred reconstruction. In this case, brain features are very well resolved (panel 3e) but are dominated by the bone artefacts.

Figure 3
figure 3

Panels (a–c) are the in situ analogue to the panels in Fig. 1(a–c) for a full head of a dead close-to-term rabbit kitten, though the region of the brain is not the same as Fig. 1. CTs were acquired at 24 keV with a propagation distance of 5 m. (a) A single propagation-based projection image of a whole rabbit kitten head suspended in agarose. (b) Phase contrast tomogram of the rabbit kitten head in (a,c) Phase retrieved tomogram of the head from (a) and (b) created using the single-material algorithm with respect to a bone/air interface. d) Phase retrieved tomogram using the single-material algorithm with respect to a tissue/air interface. (e) A close-up of the square region outlined in (d,f) Phase retrieved tomogram using the two-material algorithm (Equation (15)) with respect to the bone/tissue interface. (g and h) Profiles are plotted for each phase retrieval method for two regions on the left and right side of the images, marked ‘L’ and ‘R’, respectively. The profiles in (g) run across a section of the skull that has not yet fully calcified and highlights the accuracy of each method, showing the peaks blurring together with inaccurate phase retrieval.

When phase retrieval is performed with respect to the bone/tissue interface using Equation (15) and μ and δ values derived from the NIST physical reference databases27,28δμ = (δ2δ1)/(μ2μ1) = 5.66 × 10−10), both interfaces are better resolved with a more consistent resolution across the image than with the single-material algorithm (δ/μ = 1.54 × 10−9 and δ/μ = 8.60 × 10−9 for bone and brain tissue, respectively). These results are shown in panel 3 f, where it can be seen that the artefacts due to the highly absorbing bone are largely abated. The remaining artefacts fall into two types - ring artefacts and streak artefacts - and are ones which should ideally be corrected for prior to the phase-retrieval step. Streak artefacts across high-contrast edges are caused by multiple factors, and determining the dominant contributing factors is important for improving the SNR and CNR for in situ brain imaging.

Streak artefacts in CT

Numerous phenomena contribute to the formation of streak artefacts in CT reconstruction, particularly along high contrast edges. The three main causes are (1) insufficient data within the dynamic range of the detector, (2) noise, (3) physical effects that reduce contrast and resolution. These all create strong artefacts when there are large attenuation gradients within the sample, particularly when trying to detect subtle variation in soft-tissue contrast (e.g. grey/white matter). For a brain in situ, artefacts from the bone that are not usually problematic become important, since they can overwhelm the underlying, relatively low-contrast brain structures. Even a sub-pixel offset in the centre of rotation correction can create obvious tuning-fork artefacts29 from the bone that are not apparent when viewing an intensity palette that is scaled to include higher-density structures and hence a larger contrast range. The following phenomena are the primary contributors to streak artefacts in CT:

  • Photon starvation: This phenomenon occurs when an insufficient number of photons reach the detector through a highly-attenuating region of a sample, such as a metal implant or bone. When this occurs, the signal at that part of the detector is dominated by noise, leading to streak artefacts on reconstruction. There are a number of methods employed to correct for these artefacts, many of which involve thresholding and interpolation of the sinogram. These methods work quite well for compact regions where the region can easily be removed by thresholding in the sinogram without interfering significantly with neighbouring regions (for a discussion of these techniques, see Mouton et al.30).

  • Beam Hardening: The relationship between the logarithm of the X-ray intensity transmitted through a sample and the sample thickness is linear at any given energy (Beer’s Law) under the projection approximation; however, this is only strictly true for a monochromatic source. Lower-energy photons are absorbed more readily than higher-energy photons, leading to a deviation from this linearity for a polychromatic source (i.e. beam hardening). This change in the attenuation profile from the monochromatic case results in cupping artefacts through individual materials and streak artefacts along edges between different materials in the tomographic reconstruction.

  • Energy Harmonics: For synchrotron sources, monochromators are used to limit transmission of the source to a narrow band of energies. Energies outside the desired value can cause a form of beam hardening if the spectral bandwidth is too broad or if higher-order harmonics of the monochromator crystal are allowed to pass through. In the former case, a standard beam-hardening correction can be applied, while the latter requires a correction that accounts for the specific energies associated with the higher order harmonics.

  • Compton scattering: Compton scattering occurs as X-rays pass through the sample and can contribute to streak artefacts by increasing the background intensity, hence decreasing the contrast across edges. The scattering cross-section is lobed in the forward direction but amounts to a relatively uniform contribution for the energies and fields of view typically used for small animal imaging. This effect diminishes on propagation and is expected to have a small effect for PBI imaging at object-to-detector propagation distances over 2 m.

  • Poisson Noise: Noise in an imaging system has a similar effect to Compton scattering, with the exception that noise exhibits substantial spatial variations while Compton scattering is typically spatially smooth. Like Compton scattering, it will vary the background signal, resulting in a change in the attenuation gradient across high-contrast boundaries. This leads to an over- or under-estimation of the attenuation in the regions tangential to these edges in the CT acquisition plane.

  • Point Spread Function (PSF): Imperfect detector response plays a significant role in streak artefacts, since it results in a decrease in image resolution, leading to blurring across edges. Deconvolution of the PSF sharpens the edges, but since it also amplifies noise, there is a tradeoff in doing so. This is covered in more detail in the results section under ‘simulation and experiment’.

  • Edge Effects: The discrete nature of detector pixels means that there is a spatial averaging of the signal at each pixel when a sharp edge lies within the boundaries of a pixel31. The relevance of this effect increases with increasing pixel size, so it can be minimised with the use of high-resolution detectors. Edge effects are also tied in with the detector PSF and can therefore be accounted for, to some degree, in the process of deconvolution.

  • Centre of rotation offset: When acquiring a parallel-beam CT, it is nearly impossible to position the sample such that the centre of rotation corresponds exactly to the centre or edge of a pixel, resulting in an offset that must be accounted for when reconstructing. This is relatively easy to do, and there exist a number of different methods to determine the offset (for discussion of some of these methods, see Vo et al.32). In most cases, it is sufficient to determine this offset to the nearest pixel. When it is not sufficient, the effect can be minimised to a local blurring, rather than streaks, by acquiring projections across 360° rather than 180°. Alternatively, we find that by mirroring the projections of a 180° CT and reconstructing the volume as if it were acquired over 360°, the same effect can be achieved.

Compton scattering can be ruled out as a significant contributor to the streak artefacts in our experiments, since the scattered photon density is both proportional to the sample size and inversely proportional to the object-to-detector distance33. We expect the scatter-to-primary ratio to be less than a few percent due to the very small beam size (3 × 3 cm) and large propagation distance (5 m) at relatively low energy (24 keV). Several of the other effects described above are discussed further in the results section.

Ring artefacts

Synchrotron CT is also subject to ring artefacts, which are due to temporal variations in the intensity of the beam and possibly also non-linearities in the detector response, both of which prevent proper correction using a ‘flat-field’ image. Minimising ring artefacts is particularly challenging for in situ brain imaging, since many ring removal methods exploit their circular symmetry (see Münch et al.34 and Prell et al.35 for examples); since the skull is also largely circularly symmetric, this makes it difficult to decouple signal from artefact, resulting in an over-correction of some skull regions, where sections of bone are interpreted as artefacts, and an associated under-correction of regions that are diametrically opposite. While the development of effective ring artefact correction algorithms is essential for preclinical studies of the brain, it is omitted here as beyond the scope of this work; however it should be noted that a better option may be to prevent these artefacts from forming in the first place by shifting the sample in regular manner with respect to the detector (e.g. by translating the sample or detector vertically) during acquisition of the CT. This would prevent the variations in pixel sensitivity that are responsible for the ring artefacts from being consistently present in a single reconstruction plane. Another possibility is to address the problem where it initiates by more accurately characterising the detector components, as well as any time fluctuations in the X-ray source, in order to properly account for them.

Results

Streak artefacts from limited spatial resolution: simulation and experiment

To explore the effects that the detector PSF has on CT imaging, an absorption contrast CT of an aluminium rod in air was simulated, assuming an X-ray source energy of 26 keV and using the corresponding attenuation coefficient for aluminium, μ = 445.23 m−1. The projection images were then convolved with a 2D Gaussian (σ = 30.2 μm) before back projection. The simulation results are shown in Fig. 4, where panel 4a is the initial ideal object and panel 4b is the resulting reconstruction after PSF blurring. The latter is still a reasonable representation of the ideal object when the intensity palette is scaled to display the full intensity range; however, when the scaling is adjusted to highlight the subtle variations in the background region (4c), dark streak artefacts are clearly evident emanating from the high-contrast edges, with larger effects along longer edges.

Figure 4
figure 4

(a) A simulated absorption contrast CT of an aluminium rod surrounded by air. (b) The same rod from (a) reconstructed from projections that were blurred with a 2D Gaussian. (c) The same blurred rod from (c) with the palette scaled to focus on features in the background. The simulation was designed to have similar projected thicknesses and length to the aluminium phantom.

Figure 5 shows a reconstruction of the aluminium phantom, measuring approximately 5 mm in length in the CT acquisition plane. An energy of 26 keV was used for imaging to ensure sufficient transmission through the relatively dense aluminium. Panel 5a shows the image scaled to include the full greyscale palette in the same way as the simulation in panel 4b, clearly delineating the aluminium strip with only minimal artefacts. When scaled to the background palette (panel 5b), the extent of the artefacts is more apparent. In the inset image, dark streak artefacts can be seen that resemble those in the PSF simulation of panel 4c. The detector line spread function was measured to sub-pixel accuracy in two orthogonal directions using a 250 μm thick tungsten edge and a Pearson type VII function36:

$$y=c{[1+\frac{{(x-{x}_{0})}^{2}}{m{a}^{2}}]}^{-m}.$$
(16)
Figure 5
figure 5

(a) A reconstructed slice through a phantom consisting of a 0.3 mm thick aluminium strip suspended in agarose, scaled to the full greyscale palette. (b) The same phantom as in (a), with the palette scaled to the background region containing artefacts. (c) The same phantom after Wiener deconvolution of the measured point spread function (see the simulation and experiment section for PSF details). The phantom CTs were acquired at 26 keV at a propagation distance of 12 cm.

The PSF was then approximated by a 2D elliptical Pearson VII function:

$$y=c{(1+\frac{1}{m}[\frac{{(x-{x}_{0})}^{2}}{{a}_{h}^{2}}+\frac{{(y-{y}_{0})}^{2}}{{a}_{v}^{2}}])}^{-m},$$
(17)

where (x0, y0) is the centre of the PSF image, m is the larger of the two parameters m from the orthogonal line spread function fits, and ah and av are the amplitudes of the horizontal and vertical fits, respectively. The Pearson type VII function was chosen to represent the PSF, because it ensures that the denominator in the deconvolution is always greater than zero. The exact parameters used for this fit are c = 0.0701, m = 3.222, ah = 1.443, and av = 1.477. The PSF was found to be much broader than expected, with visible effects across sharp edges spanning ~80 pixels. The full-width at half-maximum (FWHM) is ~2.6 pixels. When the PSF is deconvolved using Wiener deconvolution, the effect of the streak artefacts is improved, but there is clearly still a large component remaining (panel 5c), indicating that other effects are involved.

These results show substantial artefacts that dominate the background of the image in the regions immediately surrounding the metal phantom. In many ways, the skull is more forgiving in that its roughly circular symmetry results in an averaging out of many of the streak artefacts in addition to restricting the bulk of them to the region outside its boundary; however, the skull is also the biggest obstacle in correcting those artefacts as it completely surrounds the brain. This fact rules out the use of many of the most common metal artefact reduction techniques available, which often involve some form of thresholding and in-filling of the sinogram, that might otherwise be applied to reduce bone artefacts. Even with as precise as possible thresholding and adaptation of a normalisation method similar to that used by Meyer et al.37, we find that the skull is too broadly distributed across the sinogram to differentiate it from soft tissue for in-filling. Iterative reconstruction methods, however, may prove to be useful for artefact reduction30.

Harmonic contamination measurement

Finally, a harmonic contamination test was performed to determine the significance of the monochromator harmonics. Images were acquired with aluminium sheets of increasing thickness placed in front of the detector until the measured signal fell to that of the dark current. The negative natural logarithm of the mean normalised intensity, \(-{\mathrm{log}}_{{\rm{e}}}(I/{I}_{0})\), was plotted as a function of aluminium thickness and fit, using least squares minimisation, as:

$$-{\mathrm{log}}_{e}(I/{I}_{0})=-{\mathrm{log}}_{e}[\mathrm{(1}-a)\exp (\,-\,{\mu }_{F}T)+a\exp (\,-\,\,{\mu }_{H}T)],$$
(18)

where μF and μH are the attenuation coefficients corresponding to the fundamental and third harmonics of the Si (1 1 1) monochromator crystal, T is the thickness of the aluminium, and a is the fractional contribution of the fundamental harmonic to the mean intensity at the detector. The fit can be seen in Fig. 6, with a value of a < 1%. The X-ray conversion efficiency decreases substantially with increasing energy38; as expected, we find the contribution of the third harmonic to be an insignificant factor for the streak artefacts created on back projection.

Figure 6
figure 6

The negative natural logarithm of the normalised intensity as a function of aluminium thickness for a nominally 26 keV spectrum of synchrotron X-rays filtered by the (1 1 1) reflection from a Si monochromator. For a purely monochromatic source, this is expected to be a straight line as per Beer’s Law. This deviation from linearity is due to contribution from the third harmonic of the monochromator crystal, which becomes dominant at larger thicknesses. The slope of the curve below ~20 mm corresponds to the attenuation coefficient of Al at the fundamental energy, and at larger thicknesses it corresponds to that of the third harmonic.

The phenomena outlined in the previous few sections are the most common causes of streak artefacts in CT. These simulations and experiments have shown that the artefacts cannot be easily explained by one individual source. The phenomenon that clearly does contribute - the PSF - is not sufficient to account for the full extent of the artefacts. We conclude that there are likely a number of different phenomena contributing that, while occurring individually might only have a small affect, together amount to a much larger cumulative effect. Isolation and correction of these effects remains the subject of future work.

Rabbit kitten brain CT

The full volume of the in situ rabbit kitten data set was reconstructed using filtered back projection (FBP) and rotated to create axial, sagittal, and coronal views for both absorption contrast and phase retrieved PBI CT. Absorption contrast CTs were acquired with the detector as close to the sample as physically possible, while phase contrast was achieved through propagation of the wavefield. Phase retrieval was performed using μ and δ values derived from the NIST physical reference databases27,28δμ = 5.66 × 10−10). Several slices of the absorption contrast CT volume can be seen in Fig. 7. The views are denoted in blue, red, and green for axial, sagittal, and coronal, respectively, and the crosshairs on each panel correspond to the locations from which the other two views in that row are cut. Note that all of the images in Fig. 7 are conspicuously featureless apart from some bone. However, in the corresponding phase retrieved views in Fig. 8, grey and white matter boundaries are resolved, and several specific brain features are clearly delineated. The overall SNR gain with phase contrast brain CT was found to be 19.7 ± 1.5. This was determined using six of the flattest regions in each image, where the SNR is the ratio of the mean to the standard deviation of each region. This number may seem surprisingly small, given the marked improvement between Figs 7 and 8; however, since the radiation dose required is inversely proportional to the square of the gain9, we can see that ~400 times more dose would be required to obtain the same result with conventional absorption contrast CT.

Figure 7
figure 7

Several axial (a,d and g), sagittal (b,e and h) and coronal (c,f and i) absorption contrast tomograms from in situ brains from dead rabbit kittens. The crosshairs on each image denote the locations of the slices from the other orientations in each row, with axial, sagittal, and coronal slices marked in blue, red, and green, respectively.

Figure 8
figure 8

The same slices from Fig. 7, now with phase contrast and phase retrieval. As with Fig. 7, axial, sagittal, and coronal views are marked in blue, red, and green, respectively. Circularly symmetric ring artefacts can be seen as a white blurring toward the bottom of panels (c) and (i). These also manifest as a diffuse white band in panel (b), running from the lower centre of the image toward the cerebellum (Cb). Note that these artefacts are not visible in the absorption contrast images in Fig. 7, since there they are below the level of the noise. Reference labels delineated in coronal sections are observed at the level of the frontal cortex (FCx) in images (a) and (g), and level of the parietal cortex (PCx) in (d), showing the caudate nucleus (CN), hippocampus (Hip), thalamus (Th) and corpus callosum (CC).

Remarkably, the brain structures visible in these images are not obscured by streak or ring artefacts. This is due to the volume being rotated with respect to the CT acquisition plane in which the artefacts are created; hence the artefacts are minimised. In panel 8a, the frontal lobe, frontal cortex, and striatum can be seen. Panel 8d shows the parietal cortex, hippocampus and thalamus, and in panel 8 g the frontal cortex, corpus callosum, and caudate nucleus are all clearly resolved. The strongest streak artefacts can be seen in the axial images. This is because the axial orientation corresponds to only slightly acute angles with respect to the acquisition plane, while the angles of the other orientations are much larger.

In general, streak artefacts pose a particular problem for brain PCXI-CT due to the very low contrast between structures within the brain. Streaks that are not distinguishable above the noise in absorption contrast CT can dominate brain structures once that noise has been suppressed on phase retrieval, since they often display similar or even higher contrast. The many possible causes of these artefacts can also be very difficult to decouple. Nevertheless, we find that by paying particular care to the orientation of the sample in the CT acquisition plane, we can minimise these effects. Some residual streak artefacts still persist, and further investigation is required to determine the most appropriate method by which to correct for these effects. For preclinical studies (e.g. tissues from deceased animals or brief terminal experiments under anaesthesia) where dose is not an issue, or for studies involving non-biological samples, it should be possible to eliminate all of these artefacts by acquiring two or three CTs in orthogonal orientations.

Ring artefacts can also be seen in the lower part of the coronal images of Fig. 8. As with the streak artefacts, these are minimised by sectioning slices at an angle with respect to the CT acquisition plane. For this volume, the coronal orientation is offset by 40° from the acquisition plane. This means that fewer consecutive image pixels contain ring artefacts originating from the same detector elements, reducing structure in the artefacts. The artefacts, however, persist to some degree along the readout direction of the detector. This results in a diffuse band that can be seen across panel 8b from the lower centre of the image, running diagonally upward and to the right, through the cerebellum. This effect is most prominent at the centre of rotation of the sample and becomes less so at larger radii.

It should also be noted that it is of particular importance for phase contrast brain imaging to accurately account for all of the phenomena that cause each of the different types of artefact. There are many correction methods that work effectively for absorption contrast CT imaging that are insufficient for the low-contrast edges that are enhanced using phase contrast imaging, as the latter exposes effects of these artefacts that, while present in absorption contrast imaging, are not generally problematic as they are below the noise threshold. With further development of artefact correction methods (e.g. PSF deconvolution, improved detector characterisation, modelling source variations, etc.), we anticipate that the SNR and CNR of phase contrast brain CT can be even further improved.

Conclusions

Propagation-based PCXI-CT is shown here to be an effective tool for visualising the brain in situ for preclinical animal studies. While the surrounding skull, temporal fluctuations in the intensity of the source, and detector imperfections provide distinct challenges with respect to reconstruction artefacts, we find that there are ways to work around these limitations to see brain structures that might otherwise be obscured. We present a two-material phase retrieval algorithm for tomography, which was shown to be highly effective at delineating soft tissue from bone. In addition, we find that by changing the ‘sectioning angle’ of the 3D volume, we are able to significantly reduce contamination from streak and ring artefacts. We have identified the most problematic causes of these artefacts, and while further work will be required to address these phenomena, it is clear that it is already possible to identify structures that have previously been unresolved with conventional X-ray CT. The substantial SNR gain achieved using PCXI-CT has no requirement for contrast agents and allows for the visualisation of features that would otherwise require a ~400-fold increase in the radiation dose required to obtain the equivalent results with conventional absorption contrast CT.

Methods

This experiment used rabbit kittens that had been used in experiments conducted with approval from the SPring-8 Animal Care (Japan) and Monash University (Australia) Animal Ethics Committees. All experiments were performed in accordance with relevant guidelines and regulations. The kittens were humanely killed in line with approved guidelines and the carcasses scavenged for this experiment.

To examine CT streak artefacts from strongly-absorbing samples, simulations were performed of a CT of an aluminium rod with a length and thickness designed to mimic those of the rabbit kitten skulls in the CTs discussed below. For comparison, a CT was experimentally acquired of an aluminium phantom at 26 keV at the shortest feasible propagation distance of 12 cm. The phantom consisted of a 0.3 mm thick strip of 99% pure aluminium sheeting suspended in agarose and was specifically designed to aid determination of the primary cause of the streak artefacts seen in in situ brain imaging by mimicking the projected thickness of the highly attenuating skull. The phantom CT was acquired using a 4000 × 2672 pixel Hamamatsu CCD camera (C9300-124) with a tapered fiber optic bonded between the CCD chip and the 20 μm thick gadolinium oxysulfide (Gd2O2S; P43) phosphor, with an effective pixel size of 16.2 μm. Each CT consisted of 1800 projections spanning 180° of rotation, with an exposure time of 80 ms per projection.

To visualise the brain, CTs were acquired of a scavenged New Zealand White rabbit kitten head and an excised rabbit kitten brain, both at 30 days gestational age (GA; term ~32 days), suspended in agarose. They were acquired at an energy of 24 keV and at a 5 m sample-to-detector propagation distance using a 2048 × 2048 Hamamatsu digital sCMOS camera (C11440-22C) with a 25 μm thick gadolinium oxysulfide scintillator and a pixel size of 15.1 μm. 1800 projections were acquired over 180° for both CTs, with an exposure time of 200 ms and a dose rate of 22.4 mGy/s. Further CTs were acquired at a higher resolution in order to test the effects of detector resolution on streak artefacts. These consisted of a rabbit kitten brain in situ, excised from a scavenged New Zealand White rabbit kitten at 29 days GA, suspended in agarose, at propagation distances of 5 m and 11 cm. Both were captured at an exposure time of 83.3 ms using a second 2048 × 2048 Hamamatsu digital sCMOS camera (C11440-22C) with a straight fibre optic and a 15 μm thick gadolinium oxysulfide scintillator with a pixel size of 6.49 μm. Due to the projected sample size being larger than the detector field of view, 7200 projections were taken through 360° of rotation, which were later stitched together with linear blending to create a single dataset of 3600 projections spanning 180°, with a dose rate of 33.9 mGy/s.

All CTs were acquired at a source-to-object distance of 210 m on beamline BL20B2 at the SPring-8 synchrotron in Japan and were reconstructed using FBP. Reconstructions were performed on the MASSIVE supercomputer in Melbourne, Australia using the ASTRA Toolbox CUDA accelerated FBP algorithm39.

Data Availability

The datasets supporting the findings of this study are available from the corresponding author on reasonable request.