Introduction

X-ray imaging, through non-invasive single projection and Computed Tomography (CT), can provide excellent structural detail of samples and is a crucial tool in medical diagnosis, but it exposes the sample to harmful ionising radiation in the process. Although significant steps have been made to reduce the radiation dose, such as using iterative reconstruction methods1, spiral scans, and selective slice scans2, the incorporation of phase contrast (PC) techniques can be used to significantly reduce dose even further3. While several phase contrast techniques exist, we focus here on propagation-based phase-contrast imaging (PBI), taking advantage of the simplicity in optical design, which only requires a sufficiently coherent beam at the sample position and some propagation distance between the object and detector. Phase contrast appears as intensity fringes on the detector and arises from the interference of X-rays that have incurred different phase shifts after passing through different materials. This effect enables PBI to be used effectively to enhance contrast between low-Z materials. To reconstruct the sample properties from the images requires application of an appropriate phase retrieval algorithm. Since we cannot directly measure the phase of an X-ray wavefield, it must be derived from intensity images alone. Phase retrieval is therefore an ill-posed problem and few algorithms provide an accurate and robust solution. The phase retrieval algorithm of Paganin et al.4 provides such a solution for the case of homogenous objects, that is, those comprised of a single material, and converts the phase contrast effects to a signal-to-noise-ratio (SNR) boosted version of the absorption contrast image. This enables low-dose distinction between low-Z materials3,5. The boost in SNR comes from the phase retrieval applying a customised blurring filter to the images, returning the phase-contrast sharpened boundaries to more directly match the object shape while removing high frequency noise in the process6. The most common algorithm for this is a Lorentzian Fourier filter, based on the transport of intensity equation (TIE), which is referred to herein as the Paganin Method (PM).

The PM relies on the paraxial approximation, requires a homogeneous object, and a small propagation distance between the sample and detector \(\Delta\) so that the resulting image is in the near-field regime. While material inhomogeneity leads to erroneous thickness determinations in projection-based images, the single-material requirement was later loosened by Beltran et al.7 to allow for dual-material samples. See also a simplification when combined with CT in Croton et al.8, which allows phase retrieval to be performed on interfaces between any two known materials, regardless of how many other materials are present in the sample. In its original form, the Paganin method of phase retrieval reconstructs the contact intensity as

$$\begin{aligned} I_{\text{PM}}(x, y)&= {\mathscr{F}}^{-1}\left\{ \frac{{\mathscr{F}}[I(x, y, z = \Delta )/I_0]}{1 + \frac{\delta \Delta }{\mu } {\textbf{k}}_{\perp {\text{PM}}}^2}\right\} , \end{aligned}$$
(1)

where \(I_0\) is the incident wavefield intensity, \(\delta\) is the decrement component of the complex refractive index \(n = 1 - \delta + i \beta\), and \(\mu\) is the linear attenuation coefficient, which is related to the imaginary component of the complex refractive index \(\beta\) and X-ray wavelength \(\lambda\) by \(\mu = 4\pi \beta / \lambda\). The vector \({\textbf{k}}_{\perp {\text{PM}}}^2\) represents the perpendicular components of the wavevector, defined as

$$\begin{aligned} {\textbf{k}}_{\perp {\text{PM}}}^2 = k_{x}^2 + k_{y}^2, \end{aligned}$$
(2)

with (\(k_{x}\), \(k_{y}\)) the Fourier spatial frequencies corresponding to the image coordinates (x, y).

In conjunction with CT, Eq. (1) has been applied in a range of settings, including investigations of self-healing thermoplastics9, studies of sandstone10 and coal micro-structure11 as well as in animal imaging studies such as detecting iron oxide nanoparticles in mouse brains12,13, and dose optimization of lung microtomography3,14. Of particular interest are the studies that report improvements to spatial resolution after incorporating sharpening filters to the algorithm, effectively suppressing high spatial frequencies to a lessor degree than implemented by the Fourier Lorentzian filter in Eq. (1). These studies included the ANKAphase15 version 2.1 software package implementation of the PM incorporating a deconvolution filter, and the addition of an unsharp mask to the pyHST2 implementation16. Alternative methods have also included adding high spatial frequency information from the phase contrast image back into the retrieved images17, intended as a compromise between phase retrieval and phase contrast. This motivated Paganin et al.18 to revisit the derivation and find a first principles justification for the success of these approaches18. Previous, first-principles-supported methods for increasing the spatial resolution of the PM have broadened the algorithm’s scope, such as by reducing the filter strength to account for inherent blurring by the system point spread function (PSF)19. However, with an algorithm based on fundamental wave optics, namely the new Generalised Paganin Method (GPM) in18, we can expect a correction that more accurately restores high spatial frequencies than post-hoc sharpening filters.

The over-suppression of high spatial frequencies by the PM is a consequence of applying Eq. (1), derived using continuous Fourier transform integrals, to discrete pixel-based imaging systems. Note that, although Eq. (1) includes Fourier transforms, this algorithm is implemented using discrete Fourier transforms, since digital images have discrete sampling. In the GPM derivation presented by Paganin et al.18, this is addressed by employing a 5-point approximation of the transverse Laplacian20,21 into which the discrete representation of the Fourier transform is substituted. This representation is then used, in place of the Fourier derivative theorem, to expand the Laplacian that appears in the derivation of the original PM. The contact plane intensity in the GPM is given as

$$\begin{aligned} I_{\text{GPM}}(x, y)&= {\text{DFT}}^{-1}\left[ \frac{{\text{DFT}}[I(x, y, z = \Delta )/I_0]}{1 + \frac{\delta \Delta }{\mu } {\textbf{k}}_{\perp {\text{GPM}}}^2}\right] , \end{aligned}$$
(3)

where \({\text{DFT}}\) represents the discrete Fourier transform mapping from the image coordinates (x, y) to the Fourier space coordinates \((k_{x},k_{y})\), and \({\text{DFT}}^{-1}\) is the inverse transform. Equation (3) explicitly accounts for the discrete image sampling at coordinates (x, y), while possessing a similar form as the PM in Eq. (1) differing only in the quantification of the spatial frequencies, now represented by

$$\begin{aligned} {\textbf{k}}_{\perp {\text{GPM}}}^2 = \frac{-2}{W^2}[\cos (W k_{x}) + \cos (W k_{y}) - 2], \end{aligned}$$
(4)

where W is the pixel size of the detector. As demonstrated in Paganin et al.18, a Taylor series expansion of Eq. (4) shows convergence to Eq. (2) for spatial frequencies close to the origin of Fourier space, but this can differ greatly when approaching the Nyquist frequency of the Fourier transform. Paganinet al.18 demonstrated that the GPM filter always suppresses high spatial frequencies to a lesser extent than the PM, providing some first principles justification for the spatial resolution improvement found by users applying sharpening masks alongside the PM algorithm18. Paganin et al.18 provided a comparison between the PM and GPM spatial frequency filters by varying \(\delta \Delta / \mu W\) over several orders of magnitude. They also demonstrate a difference between images reconstructed with the GPM and PM algorithms by subtracting one from the other, but did not consider how the PSF may affect the utility of the GPM. Our study aims to quantify the improvement in spatial resolution of the GPM by looking at the PSF of the entire imaging and reconstruction system, and identifying the experimental conditions under which the GPM provides noticeably improved spatial resolution compared to the PM.

We begin in “Analysis of system transfer functions” section by analyzing the respective transfer functions of each algorithm, comparing their spatial frequency filters before incorporating the system PSF to describe how well the imaging and retrieval system captures the various spatial frequencies present in the sample. “Spatial resolution improvement in projection images” section then provides a brief comparison of the resolution each can achieve on a simulated phantom projection. “Indirect X-ray detectors” section explores improvements achieved in CT and compares the differences between direct and indirect X-ray imaging detectors. Having determined the best experimental configuration to exploit the utility of the GPM, we demonstrate in “Rederivation of additional phase retrieval algorithms using discrete mathematics” section that the discrete Fourier transform notation introduced in the GPM can also be used to prevent over-blurring in other phase retrieval algorithms that use DFTs. Specifically, we experimentally demonstrate more accurate image reconstruction, via improved spatial resolution, for the two-material phase retrieval algorithm of Beltran et al.7 (“Two-material phase retrieval algorithm” section) and spectral decomposition algorithm of Schaff et al.22 (“Dual-energy material decomposition” section).

Simulated analysis in projection

We begin with preliminary comparisons of the PM and GPM methods by looking at the respective spatial filters. We incorporate additional filtering to mimic other stages of imaging, to better emulate how spatial frequencies in the sample are captured in the raw image and then appear in the final retrieved sample image. These two methods are also applied to simulated projection images to measure the resulting spatial resolution of the imaging system.

Analysis of system transfer functions

The PM, Eq. (1), is ultimately a tailored spatial frequency filter, derived under the transport of intensity equation, used to blur an image such that the phase contrast seen at material boundaries is spread to reconstruct the sample thickness. Given that the only difference between the PM and GPM methods is the shape of this spatial frequency filter, this becomes our first point of comparison. From Eqs. (1) and (3), we can define transfer functions for the application of each filter to raw images as

$$\begin{aligned} {\mathscr{H}}_{\text{PM}}(k_{x}, k_{y})&= \frac{1}{1 + \frac{\delta \Delta }{\mu }(k_{x}^2 + k_{y}^2)}, \end{aligned}$$
(5)
$$\begin{aligned} {\mathscr{H}}_{\text{GPM}}(k_{x}, k_{y})&= \frac{1}{1 - \frac{2 \delta \Delta }{\mu W^2}[\cos (W k_{x}) + \cos (W k_{y}) - 2]}, \end{aligned}$$
(6)

given as Eqs. (20) and (19) in Paganin et al.18, where \({\mathscr{H}}_{\text{PM}}(k_{x}, k_{y})\) represents the amplification applied to each spatial frequency amplitude by the PM, and \({\mathscr{H}}_{\text{GPM}}(k_{x}, k_{y})\) the amplification applied by the GPM. From here, a simple comparison between Eqs. (5) and (6) can be performed by taking their ratio (Eq. (21) in2), \(R(k_x,k_y) = {\mathscr{H}}_{\text{GPM}}(k_{x}, k_{y})/{\mathscr{H}}_{\text{PM}}(k_{x}, k_{y})\) and plotting the result on a normalized axis. Figure 1a does this for a range of values of the combined parameter \(\delta \Delta /\mu W^2\), displaying the fractional difference between the GPM and PM as a shaded region bounded by the \(|k_y| \cup |k_x|= 0\) (upper bound) and \(|k_x| = |k_y|\) (lower bound) lines. For the \(|k_x| = |k_y|\) lines, seen as the lower bound of the Fig. 1a plots, we produce the one-dimensional spatial frequency axis through \(k_r = W\sqrt{k_{x}^2 + k_{y}^2}\), leading to values above \(\pi\) where \((k_x, k_y)\) extend into the corners of a square image. The shaded regions help to demonstrate the new-found asymmetry of the GPM filter, correcting the PM algorithm to deliver a more uniform treatment of all edges in real-space, regardless of their orientation. A more helpful way to reveal quantitative differences between each algorithm is through the fractional difference, defined as

$$\begin{aligned} {\mathscr{D}}(k_{x}, k_{y})&= \frac{{\mathscr{H}}_{\text{GPM}}(k_{x}, k_{y}) - {\mathscr{H}}_{\text{PM}}(k_{x}, k_{y})}{{\mathscr{H}}_{\text{PM}}(k_{x}, k_{y})}, \end{aligned}$$
(7)

which converges toward zero when the PM and GPM filters match, as opposed to the ratio of Eqs. (5) and (6) which converges to one. Figure 1b plots Eq. (7) across the \(k_x\) or \(k_y\) axis line (\(|k_y| \cup |k_x| = 0\)). We see again that the larger frequencies experience the greatest difference between the two algorithms, leading to increased spatial resolution in the GPM due to the greater proportion of high spatial frequencies, and that they are effectively indistinguishable at the lower spatial frequencies. However, Eqs. (5)–(7) only represent the transfer function of the post-image processing, whereas to better reflect experimental conditions we must also account for the blurring effect of the detector imaging system and optical system (i.e. all blurring effects aside from the propagation and phase retrieval). To do this we introduce the contrast transfer function \({\mathscr{G}}(k_{x}, k_{y}, \Gamma )\). We describe the real-space detector PSF as an azimuthally symmetric Gaussian, so that \({\mathscr{G}}(k_{x}, k_{y}, \Gamma )\) is given by the Fourier transform

$$\begin{aligned} {\mathscr{G}}(k_{x}, k_{y}, \Gamma )&= {\mathscr{F}}\left[ \exp \left( -\frac{x^2+y^2}{\Gamma ^2/2.355^2}\right) \right] \end{aligned}$$
(8)
$$\begin{aligned} = \exp (-\Gamma ^2(4\pi ^2k_{x}^2 + 4\pi ^2k_{y}^2)/2.355^2), \end{aligned}$$
(9)

where \(\Gamma\) is the a full width at half maximum (FWHM) in real space, measured in pixels, and the factors of \(4\pi ^2\) in the exponential reflect the DFT normalization convention. Combining the imaging system transfer function, Eq. (9), with the phase retrieval transfer functions, Eqs. (5)–(6), gives the complete transfer functions \(\bar{\mathscr{H}}\) as

$$\begin{aligned} \bar{\mathscr{H}}_{\text{PM}}(k_{x}, k_{y})&= {\mathscr{H}}_{\text{PM}}(k_{x}, k_{y}){\mathscr{G}}(k_{x}, k_{y}, \Gamma ), \end{aligned}$$
(10)
$$\begin{aligned} \bar{\mathscr{H}}_{\text{GPM}}(k_{x}, k_{y})&= {\mathscr{H}}_{\text{GPM}}(k_{x}, k_{y}){\mathscr{G}}(k_{x}, k_{y}, \Gamma ), \end{aligned}$$
(11)

allowing us to similarly define a new fractional difference as

$$\begin{aligned} \bar{\mathscr{D}}(k_{x}, k_{y})&= \frac{\bar{\mathscr{H}}_{\text{GPM}}(k_{x}, k_{y}) - \bar{\mathscr{H}}_{\text{PM}}(k_{x}, k_{y})}{{\mathscr{H}}_{\text{PM}}(k_{x}, k_{y})}. \end{aligned}$$
(12)

By incorporating the blurring effects of the imaging and imaging system into Eq. (7), we can expand on the filter comparisons in Paganin et al.18 and demonstrate how imaging PSFs may limit the relative spatial resolution improvement between the algorithms. Figure 1c plots Eq. (12) for the same \(\delta \Delta /\mu W^2\) values as in Fig. 1b, now including a Gaussian PSF with FWHM of 3 pixels. We see that the amplitudes of all spatial frequencies are reduced relative to panels (a) and (b), particularly at the higher spatial frequencies, and the difference between the GPM and PM is reduced overall. This predicts that there may only be a very small difference in spatial resolution between the PM and GPM methods when the PSF \(\Gamma\) is a few pixels wide, as is typical of most indirect X-ray detectors. We also see that the biggest difference is now shifted to medium spatial frequencies in this example. Figure 1d plots Eq. (12) for \(\delta \Delta /\mu W^2=10\), while varying the PSF size in pixels (\(\Gamma\)). The ‘no convolution’ trend displays the fractional difference without incorporating a PSF, \({\mathscr{D}}(k_{x}, k_{y})\), as a reference point. We observe that increasing the PSF size decreases the amplitudes of high spatial frequencies, hence likely decreasing the potential improvement to resolution available via the GPM. However, for PSF widths around 1 pixel, typical of direct X-ray detectors, such as photon counting detectors, the high spatial frequency amplitudes are still 20% increased under the GPM algorithm. This leads us to suspect that direct detectors, such as photon-counting detectors, will be best suited to benefit from the GPM algorithm, while for indirect X-ray detectors, which can possess PSF widths of two or more pixels, the improvement may be relatively minor.

Figure 1
figure 1

Comparative analysis of the PM and GPM transfer functions. (a) Displays the ratio of the phase retrieval transfer functions, Eqs. (6) and (5), for various values of \(\delta \Delta /\mu W^2\), using a horizontal and diagonal slice of the 2D filter to create a bounded region presenting the asymmetry of the GPM filter. (b) Plots the horizontal, (\(k_x, k_y = 0\)), line of the fractional difference in transfer functions described by Eq. (7), used as a comparison to (c) the imaging system transfer function, Eq. (12), which incorporates a \({\mathscr{G}}(k_{x}, k_{y}, \Gamma )\) PSF, set as \(\Gamma =3.0\) FWHM which reduces the plot’s vertical scaling, as well as the phase retrieval transfer function. Finally, (d) directly displays the effect of varying the PSF width, \(\Gamma\), on the imaging system transfer function, for the case \(\delta \Delta /\mu W^2 = 10\).

Full size image

Spatial resolution improvement in projection images

The previous section showed how the shape of the image transfer function can vary between the PM and GPM phase retrieval algorithms as a result of adding realistic PSFs to the system. Here we quantify what effect the PSF has on spatial resolution in combination with phase retrieval via numerical simulation. We performed the comparison by simulating the propagation of a wavefield through cylindrical phantoms using the TIE (details below) until phase contrast fringes were produced. Next, we applied each phase retrieval algorithm and created an azimuthally-averaged profile of the phantom edge, which was differentiated to create a line spread function (LSF) that can be measured to evaluate the spatial resolution of the phase-retrieved image. We then use the measured resolution in the propagation-based phase contrast image, pre-phase retrieval, as a basis for our resolution comparison. Note that, while LSFs are one dimensional, and PSFs are two dimensional, both are distinct measurements of spatial resolution entities, and we will use the terms interchangeably throughout the paper. A typical LSF is well approximated by a Pearson VII function,

$$\begin{aligned} P(x)&= A\left[ 1 + \frac{4(x - x_0)^2}{\Gamma ^2}\left( 2^{\frac{1}{m}} - 1\right) \right] ^{-m}, \end{aligned}$$
(13)

where \(x_0\) is the peak position, \(\Gamma\) is the FWHM and m is the exponent that sets the position on the spectrum between Lorentzian (\(m=1\)) and Gaussian (approximated by \(m>10\)) behaviour. Finally, we use the FWHM values of the imaging system PSF measured in the phase-retrieved images from the PM (\(\Gamma _{\text{PM}}\)) and GPM (\(\Gamma _{\text{GPM}}\)) to calculate a fractional improvement in spatial resolution between the two algorithms,

$$\begin{aligned} \Gamma _{\text{I}}&= \frac{\Gamma _{\text{PM}} - \Gamma _{\text{GPM}}}{\Gamma _{\text{PM}}}. \end{aligned}$$
(14)

Our simulations used end-on cylindrical phantoms composed of water (\({0.998}\,{\text{g}}\,{\text{cm}}^{{-3}}\)) and polymethyl methacrylate (PMMA, \({1.19}\,{\text{g}}\,{\text{cm}}^{{-3}}\)) with a radius of 900.5 pixels, created on a \(2048\times 2048\) pixel array with pixel size \({\text{W}} = {25}\,\upmu {\text{m}}\). The wavefield directly after transmission through the phantom was constructed using the projection approximation23 on a \(\times 5\) up-sampled grid, assuming an object thickness of 6 mm, and propagated with the TIE until a single phase contrast fringe became visible in the wavefield intensity (4 mm). A small Gaussian blurring filter (\(\Gamma = 1.0\) pixel) was applied to the thickness map pre-propagation to suppress artefacts arising from the pixelated boundaries of the circular phantom, before the second blurring filter was applied post-propagation to simulate a detectors with varying PSFs. Phase retrieval was then performed using either the PM or GPM methods. Figure 2a shows example imaging system PSF measurements, before and after phase retrieval, incorporating a \({\mathscr{G}}(k_{x}, k_{y}, \Gamma )\) component, simulated through a Gaussian blurring filter, applied after the TIE propagation24.

Figure 2
figure 2

(a) Azimuthally averaged imaging system Line Spread Functions (LSF)s of the circular phantom image showing the effect of phase retrieval on spatial resolution from the phase contrast (PC) and phase retrieved images using the PM and GPM algorithms for a sample composed of water. Underlying dashed curves represent Pearson VII fits used to measure the LSF width. The phase contrast PSF was rescaled vertically by a factor of 8 for plotting. (b) Plots the percentage improvement in resolution, according to Eq. (14) of the GPM, plotted against the initial resolution of the simulated object.

Full size image

From Fig. 2b we see, for these low-Z materials, a \(\sim {6}\%\) improvement in the resolution when the phase contrast PSF FHWM is equal to the pixel size. This benefit reduces with increasing detector PSF width. At 2 pixels wide, a \(\sim {2}{\%}\) improvement is seen, and only a \(\sim {1}\%\) improvement is seen at 3 pixels PSF FWHM. This reinforces that the benefit of the GPM method is heavily dependent on the detector PSF and will show the greatest improvement over the PM when the detector PSF width is equivalent to a single pixel or smaller.

Experimental analysis in CT

Phase retrieval is often combined with CT to provide three-dimensional separation of materials with very high signal-to-noise ratios due to its three-dimensional smoothing operation3. We recorded tomographic scans of various phantoms to compare the PM and GPM using direct and indirect pixelated detectors to experimentally verify the effects of different detector PSFs on the reconstructions. We used synchrotron radiation, which provides high coherence and a low divergence beam, which are ideal for phase contrast imaging. To reconstruct the CT data, we first performed the standard flat field and dark current correction followed by phase retrieval, then used parallel beam filtered back projection with a ramp filter.

Indirect X-ray detectors

Indirect X-ray detectors typically have a PSF of several pixels in width, introduced during the conversion of X-rays to visible light, where the divergent visible light from a single X-ray photon spreads across several pixels on the detector chip. Our investigation, which focused on indirect X-ray detectors, used a Hamamatsu Orca Flash 4.0 detector, with a \({10}\,\upmu {\text{m}}\) thick Gadox (P43) phosphor directly coupled to the sCMOS sensor, with \({6.5}\,\upmu {\text{m}}\) pixels in a \(2048\times 2048\) geometry. Combined with the imaging system at beamline 20B2 of the SPring-8 Synchrotron in Japan, with a propagation distance of \(\Delta = {2}\,{\text{m}}\) and energy of \({24}\,{\text{keV}}\), this produced a system PSF with FHWM \(\Gamma = {15.6}\,\upmu {\text{m}}\) or 2.4 pixels. The filter shape analysis in “Simulated analysis in projection” section suggests that such conditions would likely lead to negligible improvement in spatial resolution under the GPM. To reproduce this effect in experimental data, here we explore pixel rebinning as a method to probe different pixel sizes for a given detector PSF, \({\mathscr{G}}(k_{x}, k_{y}, \Gamma )\), and hence allow the GPM to demonstrate improvement as the pixel size approaches the PSF width. Rebinning is frequently used to increase the signal-to-noise ratio in circumstances where a longer exposure time is not desired, or to achieve fast data transfer times for high-speed imaging, so is important to consider in the context of this kind of experiment. Figure 3a provides transfer function fractional differences, evaluated using Eq. (12), for the detector settings listed above, with increasing levels of rebinning. We observe that rebinning leads to larger differences between the PM and GPM at higher spatial frequencies, with the GPM potentially providing increased relative resolution at the new pixel size. Figure 3b demonstrates the effect of rebinning on experimentally recorded CT slices created from 1800 projection images acquired over \({180}^{\circ }\) rotation. As expected, analysis without rebinning produced only a negligible improvement in the width of the PSF after applying the GPM filter instead of the PM; approximately \({0.9(5)}{\%}\). Rebinning by factors of 2, 4 and 8 along each axis show further improvements in the GPM spatial resolution compared to the PM. This improvement is quantified as a change in PSF width of \({4(1)}{\%}\) when rebinning by a factor of 2, \({12(2)}{\%}\) at a factor of 4, and \({16(2)}\%\) at a factor of 8. Each level of rebinning shows improvement; however, the difference between rebinning by 2 and 4 is far greater than that between rebinning by 4 and 8, as expected from the higher proportion of high spatial frequency information allowed by the transfer function (Fig. 3a).

Figure 3
figure 3

The effect of rebinning the raw phase contrast image (\(\times n\)) upon the spatial resolution of the phase-retrieved images, using the PM and GPM methods. (a) Plots the transfer function fractional difference, Eq. (12), along the \((k_x,k_y=0)\) line, for various levels of rebinning. Here the filter conditions were chosen to initially match the detector PSF and the conditions used in (b) which shows the Pearson VII fit reconstructions of imaging system PSF measurements conducted at the outer edge of a cylindrical PMMA phantom, shown in inset.

Full size image

We note that the spatial resolution does not vary significantly when initially rebinning by a factor of two along each axis. This implies that data recording could be performed using a detector rebin setting, increasing data transfer rates and read time, without compromising resolution.

Direct X-ray detectors

Figure 1d shows that the proportion of high spatial frequencies in the optical system transfer function, which includes phase retrieval, increases as the width of the PSF decreases. From this, we predict that photon counting detectors will provide the best improvement to spatial resolution when using the GPM instead of the PM, since the detector PSF for these systems is inherently limited to the pixel dimensions. To show this, we imaged the same PMMA phantom as was shown in Fig. 3 with an Advacam Modupix photon-counting detector at a 2 m propagation distance using \({24}\,{\text{keV}}\) synchrotron radiation (SPring-8, BL20B2). The detector comprised a Timepix chip with a \({55}\,\upmu {\text{m}}\) pixel size and silicon sensor, with the threshold set to count X-rays above \({12}\,{\text{keV}}\). Figure 4a provides azimuthally averaged PSF fits at the outer boundary of the PMMA phantom using both the PM and GPM algorithms, showing a percentage improvement in spatial resolution of \({17(2)}\%\) according to Eq. (14) when using the GPM, alongside a drop in SNR from \(92\pm 14\) to \(81\pm 16\) (\(12\pm 35\%\)). In a second experiment, we used an X-Spectrum LAMBDA (Large Area Medipix Based Detector Array) 350K photon-counting detector with a \({1000}\,\upmu {\text{m}}\) thick cadmium telluride (CdTe) scintillator and a \({55}\,\upmu {\text{m}}\) pixel size to image a larger PMMA phantom. The results are shown in Fig. 4b and were recorded at the Imaging and Medical Beamline (IMBL) of the Australian Synchrotron. Although the larger phantom radius was intended to provide more points in the radial averages, we found it was slightly asymmetric and poorly polished, leading us to use only part of the radial arc in our analysis. Defining a radial origin to be from the phantom centre to the top of the image, our radial analysis used an arc from 235\(^\circ\) to 270\(^\circ\), where the shape of the phantom was consistent. These conditions provided a percentage improvement in the PSF FWHM of \({14(1)}{\%}\) when using the GPM, accompanied by a drop in SNR from \(27\pm 6\) to \(25\pm 5\).

Next, we show the application of the GPM to a complex sample of biomedical relevance using a direct detector. We collected a CT scan of the thorax of a recently-deceased juvenile rat. The lung inflation was fixed and the rat set in a plastic tube with agarose to prevent motion during imaging. Imaging was performed, again using the LAMBDA photon counting detector, at a \({2}\,{\text{m}}\) propagation distance and a beam energy of \({26}\,{\text{keV}}\) at the IMBL. The energy threshold of the detector was set to \({6}\,{\text{keV}}\) and phase retrieval was performed on the lung-air interface, using water as an analogous material. Figure 5b shows an interleaved CT slice where the red pixel columns are from a GPM-reconstructed slice, Fig. 5c, and the blue columns are from a PM-reconstructed slice, Fig. 5a. At small scales, the difference between each resolution is subtle; however, expanding the inset in Fig. 5c as in Fig. 5d, shows a much clearer difference between the two methods. Where the PM columns look blurred, the GPM columns appear sharper, with higher definition between the airways (including the alveoli) and the soft tissues. We reiterate that this straightforward improvement to spatial resolution achievable by the GPM when using direct detectors is a consequence of the single-pixel PSFs they possess and hence would equivalently improve data recorded with indirect detectors of similar PSF. Similarly, methods allowing sub-pixel resolution through localisation of charge clouds caused by X-ray interaction with the detector, performed with indirect25,26,27 and direct detectors28,29, would likely benefit further from the GPM rederivation.

Figure 4
figure 4

PSFs measured by azimuthally averaging the outer edge PMMA cylinders in CT recorded at 2 m propagation distance with \({24}\,{\text{keV}}\) beam energy using (a) a Advacam Modupix and (b) a LAMBDA photon-counting detector. Voxel size = \({55}\,\upmu {\text{m}}\). The PMMA cyclinder used in (a) contains an off-centre circular cavity while the cylinder in (b) is solid PMMA nearing the detector width in size.

Full size image
Figure 5
figure 5

CT slices of a rat lung reconstructed from phase retrieved projections using either the (a) PM or (c) GPM algorithms. (b) Shows an interleaved image of the two methods, with a blue band along the bottom in columns from the PM method and a red band along the top in columns from the GPM, with an inset region in (b) bounded by dashed green lines magnified in (d) for direct comparison of the two methods.

Full size image

Rederivation of additional phase retrieval algorithms using discrete mathematics

The GPM rederivation is achieved through a modified version of the Fourier derivative theorem. In a continuous domain, the Fourier derivative theorem is expressed as

$$\begin{aligned} \nabla {\mathscr{F}}[f(x,y)]&= - i {\textbf{k}}_{\perp {\text{PM}}} {\mathscr{F}}[f(x,y)], \end{aligned}$$
(15)

but in a discrete domain we need to incorporate the discrete representation of the Laplacian, described in18 through the 5 point approximation, giving

$$\begin{aligned} \nabla {\text{DFT}}[f(x,y)]&= - i {\textbf{k}}_{\perp {\text{GPM}}} {\mathrm{DFT}}[f(x,y)] . \end{aligned}$$
(16)

Given the improvement to spatial resolution comes from the modified Fourier derivative theorem, we infer the same benefit can be replicated in other phase retrieval algorithms that use Fourier transforms. Here, we consider two algorithms used commonly in propagation-based phase contrast imaging that also utilize the discrete Fourier transform. First, we look at the two-material phase retrieval algorithm of Beltran et al.7, followed by the dual-energy algorithm for decomposing images into their photoelectric and Compton scatter components30. We provide experimental validation of improvements in spatial resolution, achieved through implementation of the discrete Fourier derivative theorem in Eq. (16). We use these to demonstrate the benefit of applying Eq. (16) to derivations in place of Eq. (15). Note that this is not an exhaustive list of all the possible extensions, and we welcome further applications of this approach.

Two-material phase retrieval algorithm

We focus on the highly-stable two-material propagation-based phase retrieval algorithm developed by Beltran et al.7. This algorithm was designed used to correctly reconstruct the projected thickness of one material embedded within another, rather than just focussing on a single material within air. This method requires knowledge of the complex refractive indices of both materials and, for quantitative measures, their total projected thickness

$$\begin{aligned} A(x,y)&= T_1(x,y) + T_2(x,y), \end{aligned}$$
(17)

with \(T_1\) and \(T_2\) representing the spatially variant thickness of each material. However, when combined with CT, it is not necessary to isolate materials in projection as they are spatially separated by backprojection. Croton et al.8 showed that it is therefore not necessary to know A(x, y) for use with CT. This simplifies the algorithm to a low-pass filter that specifically matches the fringe enhancement effects at the boundary between the two materials, creating a noise-reduced representation of the absorption contrast image, given by

$$\begin{aligned} I_{\text{B-PM}}(x, y, z = 0)&= {\mathscr{F}}^{-1}\left[ \frac{{\mathscr{F}}[I(x, y, z = \Delta )/I_0]}{1 + \frac{(\delta _2 - \delta _1) \Delta }{\mu _2 - \mu _1}\, {\textbf{k}}_{\perp {\text{PM}}}^2}\right] . \end{aligned}$$
(18)

Modifying Eq. (18) to use a discrete representation of the Fourier derivative theorem, Eq. (16), gives

$$\begin{aligned} I_{\text{B-GPM}}(x, y, z = 0)&= {\text{DFT}}^{-1}\left[ \frac{{\text{DFT}}[I(x, y, z = \Delta )/I_0]}{1 + \frac{(\delta _2 - \delta _1) \Delta }{\mu _2 - \mu _1}\, {\textbf{k}}_{\perp {\text{GPM}}}^2}\right] . \end{aligned}$$
(19)

We imaged a cylindrical PMMA phantom with three cylindrical cavities; one filled with an aluminium rod and the other two with air. Figure 6 shows a radial PSF fit around the Aluminium inset, recorded using the LAMBDA photon-counting detector at a propagation distance of \({2}\,{\text{m}}\) and an X-ray energy of \({40}\,{\text{keV}}\) at the IMBL. Here the GPM rederivation shows an \({11(3)}\%\) improvement in spatial resolution at the aluminium-PMMA boundary compared to the original method.

Figure 6
figure 6

Pearson VII fit comparison of spatial resolution achieved through a PM and GPM application of the Beltran two-material phase retrieval algorithm, featuring a PMMA phantom with aluminium inset. Data was recorded using the LAMBDA detector at a propagation distance of \({2}\,{\text{m}}\), and spatial resolution analysis was performed using the aluminium-PMMA material interface.

Full size image

Dual-energy material decomposition

Phase retrieval algorithms for PBI have recently been developed to use the phase and attenuation components to enable material decomposition31,32 or to separate attenuation from phase effects30. The methods rely on the acquisition of PBI images at a minimum of two different X-ray wavelengths. Here we further generalise the phase retrieval algorithm of Gursoy and Das30, which uses the Alvarez-Macovski (AM) model of X-ray attenuation33, and begins from the linearised TIE

$$\begin{aligned} -{\mathscr{F}}\ln I(E)&= {\mathscr{F}}\left[ \int \mu (E)dz\right] + \Delta k_{\perp {\text{PM}}}^2 {\mathscr{F}}\left[ \int \delta (E)dz\right] , \end{aligned}$$
(20)

where I(E), \(\mu (E)\) and \(\delta (E)\) are discrete image maps of the spatially variant contact intensity, linear attenuation coefficient, and refractive index decrement and are all functions of energy, E. Integrations are performed along the optical axis, labelled as z. To change basis, \(\mu (E)\) is represented as a linear combination of photoelectric absorption and Compton scattering. Photoelectric absorption is proportional to \(E^{-3}\) while Compton scattering is proportional to both the electron density, \(\rho _e\), and the Klein–Nishina cross section, \(\sigma _{\mathrm{KN}}(E)\). Since \(\delta (E)\) is also a function of electron density \(\delta (E)=(\rho _e h^2 c^2 r_e)/(2\pi E^2)\), where h is Planck’s constant, c is the speed of light and \(r_e\) is the classical electron radius, allows the Compton and phase terms to be coupled together. This coupling allows us to solve for the attenuation and phase, or alternatively the projected electron density, using PBI images acquired using at least two energies as a matrix equation:

$$\begin{aligned} \begin{bmatrix} E_A^{-3} & \sigma _{KN}(E_A) + \chi _A k_{\perp {\text{PM}}}^2 \\ E_B^{-3} & \sigma _{KN}(E_B) + \chi _B k_{\perp {\text{PM}}}^2 \end{bmatrix} \begin{bmatrix} {\mathscr{F}}P \\ {\mathscr{F}}\int \rho _e \end{bmatrix}&= \begin{bmatrix} -{\mathscr{F}}\ln I(E_A) \\ -{\mathscr{F}}\ln I(E_B) \end{bmatrix}, \end{aligned}$$
(21)

where \(\chi _A = (h^2 c^2 r_e \Delta ) / (2\pi E_A^2)\) and \(\chi _B = (h^2 c^2 r_e \Delta )(2\pi E_B^2)\). Equation (21) can then be inverted to solve for P, the photoelectric component, and \(\int \rho _e\). Details of this inversion can be found in Schaff et al.22. Schaff et al. showed that the reconstruction of \(\int \rho _e\) is highly robust as the combined phase and Compton signals provide a solution that includes a low-pass Fourier filter term that smooths noise. Such smoothing is akin to the single image phase retrieval algorithms in Paganin et al.4 and Beltran et al.7. Here, we generalise Eq. (21) using a discretised version of the Fourier derivative theorem by replacing \(k_{\perp {\text{PM}}}^2\) with \(k_{\perp {\text{GPM}}}^2\). We focus our study on the stable electron density results.

To explore the benefit of the modified derivation in improving spatial resolution, we performed 180\(^{\circ }\) CTs of a dual-material phantom, again using the LAMBDA detector at IMBL. We used the same cylindrical PMMA and aluminium phantom described in the previous section. Given that the highly-attenuating aluminium rod was approximately \({5}\,{\text{mm}}\) in diameter, we opted for higher energies, \({30}\,{\text{keV}}\) and \({40}\,{\text{keV}}\), than in our previous studies, to reduce beam hardening artefacts. Figure 7 shows slices of the electron density maps calculated from Eq. (21) using dual-energy recordings at \({30}\,{\text{keV}}\) and \({40}\,{\text{keV}}\). Figure 7a uses a propagation distance of \({1}\,{\text{m}}\) and Fig. 7b uses \({2}\,{\text{m}}\) propagation distance. The 1 m propagation distance provides a region entirely within the validity of the TIE, but possesses significant ring artefacts due to limited SNR enhancement from phase retrieval. The \({2}\,{\text{m}}\) distance provides a stronger SNR enhancement, while sitting at the edge of the TIE validity; signified by the slight fringes on either side of the Pearson VII fits in Fig. 7b. At a \({1}\,{\text{m}}\) propagation distance, Eq. (14) gives a \({16(6)}{\%}\) improvement in the PSF, while the \({2}\,{\text{m}}\) propagation gives a \({29(2)}{\%}\) improvement. Each PSF fit uses a \({30}^{\circ }\) arc of the outer PMMA edge, (a) using the angular range \({202.5}^{\circ }\) to \({247.5}^{\circ }\) and (b) using the range \({45}^{\circ }\) to \({90}^{\circ }\). Overall, we find using the discretized Fourier derivative theorem provides increased spatial resolution when applied to the Alvarez–Macovski-based dual-energy reconstruction of electron density and would likely improve resolution for material decomposition too.

Figure 7
figure 7

PSF comparison of the dual-energy phase retrieval images showing reconstructed slices of electron density maps shown in the graph inset. PSFs have been fit with Pearson VII functions for quantifying the PSF FWHM. (a) Shows a spatial resolution improvement of \({16(6)}\%\) at \({1}\,{\text{m}}\) propagation distance, while (b) shows the same phantom at \({2}\,{\text{m}}\) propagation distance showing a \({29(2)}{\%}\) improvement, likely due to the higher SNR afforded by low-pass filtering in the phase retrieval algorithm.

Full size image

Conclusion

Paganin et al.18 presented a rederivation of the PBI phase retrieval algorithm of Paganin et al.4, implementing a discretized form of the Laplacian, and provided a theoretical grounding for its ability to increase spatial resolution in processed images. We expand on this result by first considering how the point spread function affects the algorithm. Using theory and simulation, we demonstrated that for detectors with a PSF larger than the pixel size, the spatial frequency filter used in discretized phase retrieval becomes comparable to the standard approach, hence the total gain in spatial resolution may be minimal. We verified experimentally that the magnitude of any improvement is highly dependent upon the width of the detector PSF using PBI-CT scans with both direct and indirect detectors. For detector systems with PSF widths of several pixels or more, typical of unbinned indirect detectors, we saw negligible improvement to spatial resolution due to the broad PSF suppressing high spatial frequency content before it could be boosted; however, by rebinning the same data to reduce the effective PSF width (expressed in pixels), the GPM achieved up to a \({16(2)}\%\) improvement in spatial resolution compared to PM, only after binning by a factor of 8. Conversely, detector systems effectively possessing single pixel width PSFs, commonly found in indirect detectors, showed up to \({17(2)}{\%}\) improvement in spatial resolution. We next extended the concept of the more accurate discrete Fourier formalism to other stable phase retrieval algorithms. For the two material algorithm of Beltran et al.5, we found the spatial resolution was improved by \({11(3)}\%\) using a direct detector. Furthermore, 3D maps electron density recovered using the dual-energy algorithm of Schaff et al.22 with \({29(2)}\%\) resolution improvement. We speculate that similar spatial resolution improvements could be gained for other algorithms that use the Fourier derivative theorem in a discrete setting.

Although the improvements in resolution afforded by the discrete phase retrieval algorithms may only be modest, the most important benefit is the knowledge that the object has been reconstructed more faithfully than was previously achieved. Given the increasing popularity of direct, photon-counting detectors for dose reduction and spectral imaging, and the superior performance of these discrete Fourier transform-based algorithms when using photon-counting detectors, we anticipate that these rederived phase retrieval algorithms will be of increasing importance in future phase contrast work.