Precise phase retrieval for propagation-based images using discrete mathematics

The ill-posed problem of phase retrieval in optics, using one or more intensity measurements, has a multitude of applications using electromagnetic or matter waves. Many phase retrieval algorithms are computed on pixel arrays using discrete Fourier transforms due to their high computational efficiency. However, the mathematics underpinning these algorithms is typically formulated using continuous mathematics, which can result in a loss of spatial resolution in the reconstructed images. Herein we investigate how phase retrieval algorithms for propagation-based phase-contrast X-ray imaging can be rederived using discrete mathematics and result in more precise retrieval for single- and multi-material objects and for spectral image decomposition. We validate this theory through experimental measurements of spatial resolution using computed tomography (CT) reconstructions of plastic phantoms and biological tissues, using detectors with a range of imaging system point spread functions (PSFs). We demonstrate that if the PSF substantially suppresses high spatial frequencies, the potential improvement from utilising the discrete derivation is limited. However, with detectors characterised by a single pixel PSF (e.g. direct, photon-counting X-ray detectors), a significant improvement in spatial resolution can be obtained, demonstrated here at up to 17%.

where I 0 is the incident wavefield intensity, δ is the decrement component of the complex refractive index n = 1 − δ + iβ , and µ is the linear attenuation coefficient, which is related to the imaginary component of the complex refractive index β and X-ray wavelength by µ = 4πβ/ . The vector k 2 ⊥PM represents the perpendicular components of the wavevector, defined as 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 thermoplastics 9 , studies of sandstone 10 and coal micro-structure 11 as well as in animal imaging studies such as detecting iron oxide nanoparticles in mouse brains 12,13 , and dose optimization of lung microtomography 3,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 ANKAphase 15 version 2.1 software package implementation of the PM incorporating a deconvolution filter, and the addition of an unsharp mask to the pyHST2 implementation 16 . Alternative methods have also included adding high spatial frequency information from the phase contrast image back into the retrieved images 17 , 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 approaches 18 . 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) in 18 , 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 Laplacian 20,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 where DFT represents the discrete Fourier transform mapping from the image coordinates (x, y) to the Fourier space coordinates (k x , k y ) , and 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 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 algorithm 18 . Paganin et al. 18 provided a comparison between the PM and GPM spatial frequency filters by varying δ�/µ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 given as Eqs. (20) and (19) in Paganin et al. 18 , where H PM (k x , k y ) represents the amplification applied to each spatial frequency amplitude by the PM, and H 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) in 2 ), R(k x , k y ) = H GPM (k x , k y )/H 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 δ�/µW 2 , displaying the fractional difference between the GPM and PM as a shaded region bounded by the |k y | ∪ |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 k 2 x + k 2 y , leading to values above π 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 www.nature.com/scientificreports/ 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 | ∪ |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 G (k x , k y , Ŵ) . We describe the realspace detector PSF as an azimuthally symmetric Gaussian, so that G (k x , k y , Ŵ) is given by the Fourier transform where Ŵ is the a full width at half maximum (FWHM) in real space, measured in pixels, and the factors of 4π 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 H as allowing us to similarly define a new fractional difference as 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 δ�/µ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 Ŵ 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 δ�/µW 2 = 10 , while varying the PSF size in pixels ( Ŵ ). The 'no convolution' trend displays the fractional difference without incorporating a PSF, 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.
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 propagationbased 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, where x 0 is the peak position, Ŵ 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 ( Ŵ PM ) and GPM ( Ŵ GPM ) to calculate a fractional improvement in spatial resolution between the two algorithms, www.nature.com/scientificreports/ Our simulations used end-on cylindrical phantoms composed of water ( 0.998 g cm −3 ) and polymethyl methacrylate (PMMA, 1.19 g cm −3 ) with a radius of 900.5 pixels, created on a 2048 × 2048 pixel array with pixel size W = 25 µm . The wavefield directly after transmission through the phantom was constructed using the projection approximation 23 on a ×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 ( Ŵ = 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 G (k x , k y , Ŵ) component, simulated through a Gaussian blurring filter, applied after the TIE propagation 24 . From Fig. 2b we see, for these low-Z materials, a ∼ 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 ∼ 2% improvement is seen, and only a ∼ 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 operation 3 . 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 µm thick Gadox (P43) phosphor directly coupled to the sCMOS sensor, with 6.5 µm pixels in a 2048 × 2048 geometry. Combined with the imaging system at beamline 20B2 of the SPring-8 Synchrotron in Japan, with a propagation distance of = 2 m and energy of 24 keV , this produced a system PSF with FHWM Ŵ = 15.6 µ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, G (k x , k y , Ŵ) , 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 • 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 reso- www.nature.com/scientificreports/ lution 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). 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. 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 keV synchrotron radiation (SPring-8, BL20B2). The detector comprised a Timepix chip with a 55 µm pixel size and silicon sensor, with the threshold set to count X-rays above 12 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 ± 14 to 81 ± 16 ( 12 ± 35% ). In a second experiment, we used an X-Spectrum LAMBDA (Large Area Medipix Based Detector Array) 350K photon-counting detector with a 1000 µm thick cadmium telluride (CdTe) scintillator and a 55 µ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  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 m propagation distance and a beam energy of 26 keV at the IMBL. The energy threshold of the detector was set to 6 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 indirect [25][26][27] and direct detectors 28,29 , would likely benefit further from the GPM rederivation.

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 but in a discrete domain we need to incorporate the discrete representation of the Laplacian, described in 18 through the 5 point approximation, giving 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 components 30 . 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.  www.nature.com/scientificreports/ Two-material phase retrieval algorithm. We focus on the highly-stable two-material propagationbased 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 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 lowpass 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 Modifying Eq. (18) to use a discrete representation of the Fourier derivative theorem, Eq. (16), gives 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 m and an X-ray energy of 40 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.
Dual-energy material decomposition. Phase retrieval algorithms for PBI have recently been developed to use the phase and attenuation components to enable material decomposition 31,32 or to separate attenuation from phase effects 30 . 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 Das 30 , which uses the Alvarez-Macovski (AM) model of X-ray attenuation 33  , 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:  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 m , and spatial resolution analysis was performed using the aluminium-PMMA material interface. www.nature.com/scientificreports/ where χ A = (h 2 c 2 r e �)/(2πE 2 A ) and χ B = (h 2 c 2 r e �)(2πE 2 B ) . Equation (21) can then be inverted to solve for P, the photoelectric component, and ρ e . Details of this inversion can be found in Schaff et al. 22 . Schaff et al. showed that the reconstruction of ρ 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 2 ⊥PM with k 2 ⊥GPM . 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 • 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 mm in diameter, we opted for higher energies, 30 keV and 40 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 keV and 40 keV . Figure 7a uses a propagation distance of 1 m and Fig. 7b uses 2 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 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 m propagation distance, Eq. (14) gives a 16(6)% improvement in the PSF, while the 2 m propagation gives a 29(2)% improvement. Each PSF fit uses a 30 • arc of the outer PMMA edge, (a) using the angular range 202.5 • to 247.5 • and (b) using the range 45 • to 90 • . 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.

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 www.nature.com/scientificreports/ 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.

Data availability
The datasets used and/or analysed for this manuscript are available from the corresponding author on reasonable request.