Fast, non-iterative algorithm for quantitative integration of X-ray differential phase-contrast images

: X-ray phase contrast imaging is gaining importance as an imaging tool. However, it is common for X-ray phase detection techniques to be sensitive to the derivatives of the phase. Therefore, the integration of differential phase images is a fundamental step both to access quantitative pixel content and for further analysis such as segmentation. The integration of noisy data leads to artefacts with a severe impact on image quality and on its quantitative content. In this work, an integration method based on the Wiener filter is presented and tested using simulated and real data obtained with the edge illumination differential X-ray phase imaging method. The method is shown to provide high image quality while preserving the quantitative pixel content of the integrated image. In addition, it requires a short computational time making it suitable for large datasets.


Introduction
Differential interference contrast is an established microscopy technique for the visualization of unstained biological specimens with poor transmission or reflection contrast to visible light [1].However, the obtained images cannot be directly used to draw quantitative conclusions, as their intensity is not linearly dependent on the phase distribution [2].In recent years, phase imaging was extended to the X-rays regime [3].The enhanced contrast to soft tissues provided by phase imaging, combined with the high penetration depth of X-rays, is opening interesting perspectives for pre-clinical and clinical imaging [4].However, the detection of phase changes for X-ray is challenging, since phase gradients in a sample typically lead to refraction angles of the order of microradians or less.Several techniques have been developed to translate such a refraction angle into a change of intensity recorded on the detector [5].The simplest implementation is based on a pure interference effect, i.e. propagation-based, requiring a spatially coherent X-ray source such as readily available at third generation synchrotron radiation facilities.The combination of high coherence and high flux offered by such machines provides impressive images, e.g. of whole animal organs with resolution down to the micron scale, which can be of high value for pre-clinical investigations [6][7][8][9].In the recent years, a lot of effort has been put into translating X-ray phase imaging from synchrotrons to conventional X-ray sources, leading to new phase detection schemes such as edge illumination and gratings-based [3,5,10].The availability of phase contrast techniques based on conventional sources and compact enough to fit into a standard laboratory is opening the way to applications previously restricted to synchrotrons [11,12].However, similarly to differential phase microscopy with visible light, they do not provide immediate access to phase, but only to its first or second derivative.Therefore, additional processing steps are required to quantitatively retrieve the phase.While the retrieval of the phase from its second derivative is widely investigated for propagation-based phase imaging [13,14], gratings-based and edge illumination imaging provide access to the first derivative of the phase, in one or two directions depending on the experimental setup [15,16].Generally speaking, the phase from unidirectional differential phase contrast (DPC) images is retrieved by direct integration of each image row (or column depending on the direction of phase sensitivity; from now on, we will assume this to be the horizontal direction).However, noise in the DPC image and inevitable limitations in the sampling frequency of the phase signal, leads to the formation of streak artifacts, which severely affect the quality and the quantitative content of the integrated image.In the framework of X-ray DPC imaging, several algorithms have been proposed to solve this problem, mainly based on iterative approaches [17][18][19].While they work well and largely prevent streak artifacts, the required number of iterations can make these approaches very slow on large images, possibly hindering their application to large CT datasets which usually consist of thousands of projections.In this work, we present a non-iterative integration method for unidirectional X-ray DPC images, based on an image model initially developed for differential interference contrast microscopy with visible light [20,21].The proposed method is based on the use of the Wiener filter to estimate the phase signal while keeping the noise in the direction orthogonal to phase sensitivity under control, thus preventing streak artifacts.We used simulated and real edge illumination DPC images to compare the proposed integration method with direct and iterative integration, and to test its quantitative response.We found that the method outperforms direct integration while providing a reliable quantitative estimation of image values if used on DPC images with a sufficiently high signal-to-noise ratio (SNR).At lower SNR levels, it still provides good image quality, but larger variation in the quantitative response are observed.In addition, it provides comparable image quality to the mentioned iterative approaches, but it does so in a much shorter computational time, which makes it suitable to the integration of large datasets.Finally, we used the method to integrate an experimental DPC edge illumination image of a complex specimen.

Algorithms
The algorithm proposed in this work is based on a model originally introduced to describe images acquired with visible light differential phase microscopy [20,22].However, the model can be easily adapted to data acquired using X-ray systems with differential phase sensitivity, such as edge illumination, taking into account that phase sensitivity is achieved in one direction only [23,24].In a unidirectional X-ray DPC image, the intensity value at the pixel location (x, y) is equal to the refraction angle ∆θ, which is related to the unit decrement of the real part of the refractive index δ at sample location (x,y,z) by the equation: where ϕ = ∫ δ(x, y, z)dz.Its differential nature means the DPC image can be described as the difference between the phase images calculated at sample position (x + ∆x, y) and (x, y) where ∆x is equal to the sampling step of the image in the x direction.An equivalent formulation is given by replacing ϕ(x + ∆x, y) − ϕ(x, y) with ϕ(x + ∆x/2, y) − ϕ(x − ∆x/2, y), since a shift within the sampling step does not change the value of the function.Eq. (1) can also be written as: where * denotes convolution and d(x, y) = δ(x + ∆x/2, y) − δ(x − ∆x/2, y) is the difference of Dirac's delta functions at two image positions separated by one sampling step.To obtain ϕ(x, y), the Fourier transform of 2 is considered: where ∆Θ(u, v), D(u, v) and Φ(u, v) are the Fourier transforms of ∆θ(x, y), d(x, y) and ϕ(x, y), respectively.In particular: The integrated phase ϕ(x, y) can thus be retrieved by taking the inverse Fourier transform of ∆Θ/D.However, division by D may enhance regions of the frequency spectrum where noise is dominant.Furthermore, since D is known analytically, a problem occurs at the frequencies for which 1/D has a discontinuities.A solution to this problem is offered by the Wiener filter.Given a signal in the frequency domain in the form Y(f)=H(f)X(f)+N(f), the Wiener filter W(f) provides a least square estimate of the function X(f) from the measured signal Y(f) while taking into account the noise N(f) in the input image and a convolved function H (in real space) [25].
Using the Wiener filter we can write: where W(u, v) is the Wiener filter in the frequency domain defined as: and SNR(u, v) is the signal-to-noise ratio as function of the spatial Combining Eqs. ( 6) and ( 5), the integrated phase image can be obtained as: where F and F −1 denote the Fourier transform and its inverse.We note that for practical uses SNR(u, v) is a-priory unknown, and is therefore typically approximated with an analytic function, with a cutoff parameter that allows controlling the frequency response of the filter.While 2D Gaussian functions have previously been used for this purpose [20], we propose the use of a Butterworth filter because of its better response at high spatial frequencies [26].We restricted the filter to the vertical frequency only, considering the unidirectional property of the streak artifact we aim to reduce, which also prevents any loss of details in the horizontal (phase sensitivity) direction where edge-illumination allows achieving higher resolutions by "dithering" (a form of oversampling, see [27]).The Butterworth filter used in this work is defined as: where v 0 is the vertical cutoff frequency, and n is an integer.The Wiener method discussed above is compared with direct integration and with an iterative approach [17].Direct integration is defined as the cumulative sum running along the image rows.Given a M × N DPC image, the pixel at location (x m , y n ) of the directly integrated image is obtained as: with n and m running from 1 to N and 1 to M, respectively.The iterative approach used as a comparison aims to find the integrated differential profile ϕ(x, y) as the minimum of the objective function: where ||.|| identifies the ℓ 2 norm operator, ∆θ is the experimental differential image, D x and D y are the operators taking the derivative along the horizontal and vertical directions, respectively, and λ is a regularization parameters controlling the noise in the vertical direction of the optimized phase image [17].

Simulation
The proposed integration algorithm has been initially tested on a simulated DPC image of Shepp-Logan numerical phantom obtained with an edge illumination system [23,24].The simulation is based on geometric optics, assuming δ = 7 • 10 −7 and β = 1 • 10 −9 for the brightest region of the phantom.The values in the other regions are scaled according to their relative intensity.The sample thickness is assumed to be 1 mm.The simulated edge illumination setup is based on a detector mask with 20 µm aperture size and 50 µm pitch, and a matching sample mask demagnificated by M = 1.25× [11].The source is assumed to be monochromatic with energy 20 keV and to have a Gaussian shaped focal spot size with 70µm FWHM.The detector has a pixel size of 50 µm and is characterized by a Gaussian shaped point spread function with a FWHM of 50 µm.Noise is assumed to be Poisson distributed.A full scan of the transmitted intensity as a function of the relative mask displacement, referred to as illumination curve, is preliminary obtained as: where A 1 is the the pre-sample mask and S is the source shape, both re-scaled at the detector plane, and A 2 is the detector mask [28].From the computational point of view, the masks have been defined as square waves and all the function involved have been calculated with a sufficiently fine sampling step [28].Different images at different relative position on the IC (i.e.different masks displacements) are then simulated using the Shepp-Logan phantom according to the algorithm described in Algorithm 1 in Supplement 1. Specifically, five images have been simulated with the mask shifted by ±12, ±9 and 0 microns with respect the aligned (i.e.0) position (see Supplement 1), and used to retrieve the DPC image by pixel-wise Gaussian fit [29].

Experimental setups
Data obtained with different X-ray systems have been used in order to test the reliability of the proposed method under different conditions, encompassing both synchrotron and conventional sources.A 2.9 mm diameter acrylic rod was scanned at the SYRMEP beamline of the Italian ELETTRA synchrotron (Trieste) [30], using 20 keV monochromatic X-rays.The detector pixel size was 50µm and 60µm in vertical and horizontal directions, respectively.A 1 mm diameter nylon wire was scanned at the ID17 beamline of the European Synchrotron Radiation Facility (ESRF) in Grenoble, France, using 50 keV monochromatic X-rays [31].Phase sensitivity was achieved by an implementation of edge illumination based on two free-standing tungsten masks.A full description of the setup, image acquisition procedure and phase-retrieval method has been provided in a previous work [32].SNR values equal to 22.5 and 7.9 have been found for the DPC images of the acrylic rod and nylon wire, respectively; these were calculated as SNR = I fringe /σ bck , where I fringe is the value of the brightest phase peak and σ bck is the standard deviation of the air background, both expressed in radians.Finally, a beetle was scanned with an edge illumination setup based on a laboratory source.This is a scanning system employing the same masks, detector and system geometry previously described [33] but with the tungsten source replaced by a Rigaku M007 molybdenum source operated at 40 kVp and 30 mA.

Simulated data
The proposed integration method has been tested on simulated DPC images of the Shepp-Logan phantom obtained from an edge illumination system using different photon statistics.The comparison with the ground truth, with direct integration and with an iterative integration algorithm is shown in Fig. 1. Figure 1(a) illustrates the result of the three integration methods applied to a simulated DPC image with high photon statistic.10 8 photons have been used in this case, leading to a standard deviation in the background of DPC image of σ = 1 • 10 −9 (in a 50 × 50 ROI).All methods show good agreement with the ground truth as shown in the line profiles.The horizontal line profiles perfectly match the ground truth for both approaches, indicating no loss of detail in the horizontal direction.In the vertical direction both the proposed and the iterative algorithms fail to return a perfectly flat signal and a maximum difference of approximately 24% with respect to the ground truth is found (see stars in the vertical line profile).This is due to the suppression of some of the high frequencies in the vertical direction that are required to produce a completely flat profile.The result of direct integration gets significantly worse as the photon statistic decreases.).The direct integration now shows prominent streak artifacts that are absent in the image integrated with the proposed and with the iterative algorithms.This is again easily visible in the line profiles.In particular, the horizontal line profile shows that all the integration methods leads to a good agreement with the ground truth, although a discrepancy can now be observed in the flat region where the signal is expected to be zero (see black star).The vertical profile reveals that the directly integrated image is now dominated by noise introduced by streak artifacts.While the overall signal trend can be still distinguished, some fine details are lost (see features indicated by the yellow and red triangles).Conversely, the line profiles extracted from images integrated with the proposed algorithm and with the iterative one still show a very good match with the ground truth, with fine details still clearly detectable (see yellow and red triangles).Remarkably, no evident difference is found between the images obtained with the proposed algorithm and with the iterative one, neither in terms of quantitative values nor of preservation of fine details.This situation is exacerbated in the images simulated with 10 4 photons (background standard deviation σ = 1 • 10 −7 ) in Fig. 1(c).In this case, the directly integrated image is severely affected by streaks artifacts including in terms of quantitative content.This is well visible in the line profiles.The horizontal profile for direct integrated image is fully dominated by a linear drift due to the artifacts overlapping with the sample signal.In the vertical profile neither details nor a clear trend can be identified.Conversely, both profiles obtained from the proposed and the iterative methods still show a good agreement with the ground truth, preserving image quantitative values and fine details (see yellow and red triangles).Also in this case the match between the image obtained from the proposed algorithm and the iterative one is very good.The vertical line profile reveals how the iterative algorithm returns a slightly underestimated value compared to the proposed algorithm that is closer to the ground truth value.

Experimental data
The quantitative reliability of the proposed integration algorithm has been tested using DPC images of an acrylic rod and a nylon wire acquired at different synchrotron radiation facilities using monochromatic X-rays of different energy.The results are shown in Fig. 2. Panels (b) and (d) show the results for the integration of DPC image of the acrylic rod acquired with 20 keV monochromatic radiation, with panel (a) showing the corresponding DPC image.The relatively low X-ray energy, results in a large refraction angle, and, therefore, in a high SNR value.Panel (d) show line profiles extracted from averaging the rows in the areas highlighted in panels (b) and (c) for the integrated phase obtained by direct integration (dashed red line) and with the proposed algorithm (solid blue line), respectively, compared to the ground truth (solid black line).The proposed algorithm has been applied with the following parameters: v 0 = 0.1mm −1 , s = 3 • 10 4 , n = 1.Both integration methods provide a good quantitative response, with a maximum discrepancy of only about 8% from the theoretical values.This difference may be due to material inhomogeneity.The results for the nylon wire, acquired with 50 keV monochromatic radiation, are shown in Fig. 2(e) to (f).It should be noted that the overall image quality in this The standard deviation of the background (σ) is reported on all DPC images.Line profiles in the horizontal and vertical directions are compared against the ground truth for all the methods and for all photon statistic.Black stars highlight disagreement in flat regions while red and yellow triangles indicate a small feature which disappear in direct integration as soon as the noise is increased, but is preserved by the proposed integration method and by the iterative approach.The same region is also shown in the zoom-in insets.
case is lower, with several noisy regions visible in the DPC image in panel(e).Therefore, analysis has been performed in a homogeneous region averaging over more rows than the previous case (see red and blue rectangles in panels (f) and (g)).As already observed in the simulated case, the lower SNR has a significant impact on the directly integrated image.While the maximum value is close to the theoretical one, the profile of the wire is completely distorted because of the streak artifacts (see panel (h)).On the other hand, the proposed integration method reproduces the expected profile shape much better, while overestimating the integrated phase value of about 9% compared to the theoretical value, as shown in panel (h).The proposed algorithm has been used with v 0 = 0.08mm −1 , s = 3 • 10 4 , n = 1.It is worth noting that the background value doesn't go to zero on the right hand side of the wire and this is mainly due to a non perfect sampling of the phase peak.The proposed algorithm requires selecting the cut-off frequency, a user defined parameter typically chosen by visual inspection of the integrated image.The change in the quantitative value of the image as a function of v 0 therefore needs to be addressed.Panel (i) and (j) report the percentage variation of the maximum of the integrated image as function of the cut-off frequency v 0 for both samples.This plot is usually referred to as the L-curve [34].In the plot for the acrylic rod (panel(i)), a variation of the cut-off frequency from 1 mm −1 to 0.001 mm −1 leads to a variation of the maximum of the integrated phase value from 9% to about 10.5% compared to the theoretical value.For cut-off frequencies around 1 mm −1 an insufficient amount of high frequencies is filtered out leaving streak artifacts in the image, while around 0.001 mm −1 the image appears to be over-smoothed (see Supplement 1).Optimal results seem to be obtained for v 0 around 0.1 mm −1 , which correspond to the curvature change of the L-curve, usually assumed as the best value for a regularization parameter [34].A similar analysis is shown in Fig. 2(j) for the nylon wire.In this case, a larger variation is found, probably because of the significant lower SNR.Decreasing the cut-off frequency from 0.5 mm −1 to 0.05 mm −1 leads to variations from −80% to 20% over the theoretical values.It is worth noting, again, that the outermost points represent the extreme cases, where the image is either still affected by streak artifacts, or over-smoothed (see Supplement 1).Remarkably, a good image quality and a reasonable quantitative estimation can be still found around the same 0.1 mm −1 value, e.g. the underestimation of about 9% observed in panel (h) was obtained with v 0 = 0.08mm −1 .Even though a certain degree of freedom in the choice of the filter parameters exists, the previous analysis reveals how the quantitative results obtained with the proposed integration method depend significantly on the quality of the initial DPC image, which highlights the importance of using a robust algorithm for phase retrieval [35].Finally, the proposed algorithm has been tested on a more complex biological specimens.The results are shown in Fig. 3, which shows a ground beetle acquired using a scanning-based edge illumination system implemented with a laboratory source [33].Direct integration provides a very poor results, as shown in Fig. 3(b): the image is dominated by streak artifacts, to the extent that the insect is barely recognisable despite the strong DPC signal visible in panel (a).Conversely, the result of integration using the proposed method (panel (c)) leads to a significant advantage in terms of image quality.Details that were completely invisible in the directly integrated image are now clearly visualized, as made even clearer by the zoomed-in regions on the right-hand side of panels (b) and (c).Remarkably, the image integrated with the proposed method is very similar to the one obtained with the iterative method shown in panel (d), with minimal differences being immediately observable mainly in the centre of the specimen where a lower phase sensitivity is achieved.This is also confirmed by the line profiles extracted from both images and reported in panel (e).Finally, the application of the proposed method to a biological specimen producing a weaker phase signal is reported in Supplement 1 for completeness.Fig. 3. Panels (a) to (d) shows the DPC image of a beetle, the result of direct integration, of integration using the proposed algorithm and the iterative approach, respectively [17].Zoomed-in regions corresponding to the red rectangles are reported on the right-hand side of all images.Panel (e) shows line profiles extracted from images integrated with the proposed and the iterative algorithms.The proposed algorithm is used with the following parameters: v 0 = 0.05mm −1 , n = 1 and s = 3 • 10 4 .The iterative algorithm is used with λ = 0.01 and 200 iterations.

Discussion
The presented results clearly show that the proposed method outperforms direct integration, especially when the SNR of the differential image is limited.Moreover, it provides very similar image quality to previously proposed iterative phase integration approaches [17], but with the advantage of a short computation time, which makes it more suitable to integrate large datasets such as those obtained from CT scans.For example, a Matlab implementation of the proposed algorithm running on a workstation based on a CPU Intel Xeon Gold 6134 3.20GHz, was tested on 420 × 420 and 2320 × 2320 images taking 12 ms and 250 ms, respectively.It is worth noting that, while iterative integration methods and the one proposed here look different, they share the same theoretical basis.Iterative approaches are based on the minimization of an error function.This task may require many computational resources depending on the number of parameters slowing down the convergence of the algorithm.The image integrated with the proposed method is also obtained as a least square solution, but calculated through the application of the Wiener filter.Therefore, both approaches ultimately, minimize an error function [25].This also explains why the proposed method is much faster than an iterative one while providing similar results, since the minimisation process is undertaken "a priori" by the use of the Wiener filter.Control over the frequency content of the output image after application of the Wiener filter, is achieved by approximating the signal to noise ratio required by its definition with an analytical function.
In the present work, the Butterworth filter was used.Three parameters needs to be optimized, i.e. the order n, the filter amplitude s and the cut-off frequency v 0 .While these are quantitatively related to the signal and noise content of the input image, the easiest way to find the best values is by trial and error, by visually optimizing image quality for each integrated sample.We found that the most critical parameter is the cut-off frequency v 0 , which has an impact on both the image quality and the quantitative values of the image.An excessively low cut-off frequency results in residual streak artifacts, while an excessively high one leads to an over-smoothing of the image in the vertical direction (see Supplement 1).In addition, when the SNR of the input DPC image is high, the quantitative content of the integrated image is almost independent from the cut-off frequency if it is chosen within a reasonable range.Notably, the proposed method still delivers a good image quality when the SNR of the DPC image is poor; however the variation in the quantitative content of the integrated image as a function of the cut-off frequency becomes larger, suggesting that care must be taken where quantitative conclusion are drawn.This notwithstanding, a range of cut-off frequencies for which the algorithm provides a good quantitative estimate can be still found.The amplitude s is a less critical parameter, and it is related to the value of the SNR of the input image.A low amplitude is equivalent to an input image dominated by noise, for which application of the Wiener filter returns an almost empty image as the "best estimate" of the original signal.On the other hands, an excessively high value for the amplitude makes the filter ineffective, resulting in residual streak artifacts (see Supplement 1).Finally, the order n affects the sharpness of frequency cut-off.In all the cases we explored, we found that a low value around 1 or 2 works provides the best results, with higher values leading to an heavily oversmoothed image in the vertical direction (see Supplement 1).The presented integration algorithm has been developed and tested using unidirectional DPC images obtained from edge illumination, however, other X-ray phase contrast techniques with unidirectional differential sensitivity, such as gratings or crystal based, would also benefit from its application.

Funding
Royal Academy of Engineering; Engineering and Physical Sciences Research Council (EP/T005408/1).

Acknowledgment
This work is funded by EPSRC (grant EP/T005408/1).ME was supported by the Royal Academy of Engineering under the RAEng Research Fellowships scheme.AO was supported by the Royal Academy of Engineering under the Chairs in Emerging Technologies scheme.We thank Elettra Sincrotrone Trieste for access to SYRMEP beamline (Proposal No. 20140147) that contributed to the results presented here and the staff for support during data acquisition.The authors also thank the staff of the ID17 beamline at the European Synchrotron Radiation Facility in Grenoble, for their support during data acquisition.We also thank Mr.Savvas Savvidis from the Advanced X-ray Imaging group of the Department of Medical Physics and Biomedical Engineering, University College London, for providing the oesophagus data.

Fig. 1 .
Fig. 1.Comparison of the proposed integration algorithm with direct and iterative integration using simulated DPC images of the Shepp-Logan phantom.Panels (a),(b) and (c) shows DPC images simulated using 10 8 , 10 5 and 10 4 photons and the corresponding integration results.The standard deviation of the background (σ) is reported on all DPC images.Line profiles in the horizontal and vertical directions are compared against the ground truth for all the methods and for all photon statistic.Black stars highlight disagreement in flat regions while red and yellow triangles indicate a small feature which disappear in direct integration as soon as the noise is increased, but is preserved by the proposed integration method and by the iterative approach.The same region is also shown in the zoom-in insets.

Fig. 2 .
Fig.2.Panel (a) shows the DPC image with a line profile in the inset for the 2.9 mm acrylic rod imaged at ELETTRA using 20 keV monochromatic X-rays.Panels (b) and (c) show the result from direct integration and from the proposed algorithm using the following parameters: ∆x = 0.06mmm, n = 1, s = 3 • 10 4 and v 0 = 0.1mm −1 .Panel (d) reports a comparison of line profiles from both integrated images and the ground truth.Panels (e) to (h) show the same comparison for the 1 mm nylon wire, imaged at ESRF with 50 keV X-rays.The proposed algorithm has been used with the following parameters: ∆x = 0.01mm, n = 1, s = 3 • 10 4 and v 0 = 0.08mm −1 .For sake of comparison, for each sample, the line profiles have been normalized to the maximum value expected for the ground truth.Panel (i) and (j) show, for both samples, the percentage variation with respect to the ground truth, of the maximum value of the line profiles obtained with the proposed algorithm as a function of the cut-off frequency v 0 .