A quantitative correction for phase wrapping artifacts in hard X-ray grating interferometry

X-ray grating interferometry-based computed tomography is a phase contrast imaging technique that provides non-destructive, quantitative, and three-dimensional visualization with contrast superior to traditional absorption-based techniques, especially for materials primarily composed of low Z elements, such as biological tissues. However, it relies on measurements of the lateral shift of an interference pattern and is thus susceptible to so-called phase wrapping artifacts, which mainly occur at the sample-air interface. In this work, we present an algorithm for removal of such artifacts in the case of cylindrical samples and an experiment to verify its accuracy. The proposed algorithm is applied to the sinogram after phase retrieval and prior to reconstruction by finding sample edges with the absorption sinogram and replacing regions of the phase wrapped sinogram with modeled data. Our measurements show that the algorithm removes artifacts and produces more accurate δ values, as validated by measurements without phase wrapping. Our correction algorithm allows for measurements without submerging the sample in a water bath, simplifying the experimental setup and avoiding motion artifacts from gas bubbles.X-ray grating interferometry-based computed tomography is a phase contrast imaging technique that provides non-destructive, quantitative, and three-dimensional visualization with contrast superior to traditional absorption-based techniques, especially for materials primarily composed of low Z elements, such as biological tissues. However, it relies on measurements of the lateral shift of an interference pattern and is thus susceptible to so-called phase wrapping artifacts, which mainly occur at the sample-air interface. In this work, we present an algorithm for removal of such artifacts in the case of cylindrical samples and an experiment to verify its accuracy. The proposed algorithm is applied to the sinogram after phase retrieval and prior to reconstruction by finding sample edges with the absorption sinogram and replacing regions of the phase wrapped sinogram with modeled data. Our measurements show that the algorithm removes artifacts and produces more accurate δ values, as validated by measurements ...

(Received 20 June 2018; accepted 3 August 2018; published online 31 August 2018) X-ray grating interferometry-based computed tomography is a phase contrast imaging technique that provides non-destructive, quantitative, and three-dimensional visualization with contrast superior to traditional absorption-based techniques, especially for materials primarily composed of low Z elements, such as biological tissues.However, it relies on measurements of the lateral shift of an interference pattern and is thus susceptible to so-called phase wrapping artifacts, which mainly occur at the sample-air interface.In this work, we present an algorithm for removal of such artifacts in the case of cylindrical samples and an experiment to verify its accuracy.The proposed algorithm is applied to the sinogram after phase retrieval and prior to reconstruction by finding sample edges with the absorption sinogram and replacing regions of the phase wrapped sinogram with modeled data.Our measurements show that the algorithm removes artifacts and produces more accurate d values, as validated by measurements without phase wrapping.Our correction algorithm allows for measurements without submerging the sample in a water bath, simplifying the experimental setup and avoiding motion artifacts from gas bubbles.V C 2018 Author(s).All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http:// creativecommons.org/licenses/by/4.0/).https://doi.org/10.1063/1.5045398Currently, X-ray grating interferometry (XGI) is especially promising among the phase tomography techniques because it reaches high density resolution (contrast), produces quantitative results in a model-independent manner without a priori knowledge of sample composition, and can be transitioned from synchrotron radiation facilities to laboratory and clinical environments.There are about a dozen teams worldwide that apply this phase tomography method to non-destructively visualize specimens with a resolution down to a few micrometers.
In most cases where high phase sensitivity is necessary, the specimen is fixed top-down by a manipulator so that it can be rotated in a water bath.This reduces the change in the index of refraction at the sample boundaries and eliminates phase wrapping artifacts.However, this setup is suboptimal because standard manipulators hold the sample from below.Also, in many cases, gas bubbles are formed at the waterspecimen interface, creating motion artifacts. 1,2Therefore, a measurement technique without the water bath would simplify and improve grating-based phase contrast imaging both at synchrotron radiation facilities and in laboratory and clinical settings as sensitivity increases.
We use the term "phase wrapping," as it is commonly found in the XGI community.Namely, it describes both the standard definition [3][4][5] and the more complex case arising from strong wavefront curvature wherein the interference pattern phase shift is not linearly related to the wavefront gradient.In the second case, even a perfect unwrapping of the interference pattern phase shift by standard methods will not produce a useable signal, and thus a dedicated correction algorithm is needed.One such approach has been to identify and correct phase wrapped pixels from regions where the derivative of the simultaneously recorded attenuation signal is large. 6,7This is limited for sensitive setups or certain materials where the phase signal is wrapped in regions where the absorption signal is uniform.Additionally, there has been no study to verify the quantitative accuracy of these corrections.Another approach uses maximum likelihood principles and multiple energy data to estimate the interferometer phase shift. 8However, in practice, this is time consuming and requires a flexible setup.
Our proposed algorithm works for cylindrical specimens and uses a priori knowledge of the sample geometry to model and replace the phase-contrast measurement around the entire edge of the sample.This approach prevents artifacts in the center of the sample, which are caused by the phase wrapped edges through the integration step during reconstruction.An earlier version of this technique was used in a dedicated fashion to retrieve the phase of a tumor in an Eppendorf tube taking advantage of the cylindrical shape. 9-11However, a prerequisite for the widespread adoption of this method is a general form of the correction and a detailed understanding of the quantitative accuracy of data treated by it.
It should be noted that in the context of computed tomography, cylindrical samples are very common.For example, samples that are immersed in a liquid such as formalin, phosphate-buffered saline, etc., are typically held by a cylindrical polymer container, for example, an Eppendorf container.Another common technique involves creating cylindrical punches from paraffin-embedded tissues, a standard preparation technique in histopathology for soft tissues including biopsies.
An XGI setup contains a beam splitter grating creating an interference pattern maximum at an odd fractional Talbot distance d downstream, where, in the case of a double grating setup, an analyzer grating and the detector are placed.A sample in the beam path causes a local shift, /, of this interference pattern, which is measured relative to a flat field by each detector pixel using e.g., a phase stepping method. 12his / is related to the differential phase shift of the X-ray wavefront and thus the sample composition via the decrement d of the real part of the refractive index where p 2 is the period of the analyzer grating, which is matched to the period of the interference pattern.After acquiring projections /(x, h) around 180 (or 360 ), integration along the x-axis and tomographic reconstruction can be performed to retrieve a two-dimensional distribution of dðx 0 ; y 0 Þ.This can be accomplished by a variety of methods, e.g., a filtered backprojection algorithm with a filter kernel modified to integrate the differential phase. 13ue to the periodic nature of the interference pattern cast on the detector, the measured shift of the interference pattern w is a restricted version of the actual shift: w ¼ / modulo 2p.A projection with pixels where /ðx; hÞ 6 2 ðÀp; p suffers from phase wrapping.The extent of phase wrapping can be found by solving Eq. ( 1) for j/ðx pw Þj !p under the assumption of a uniform cylindrical sample of radius R. 14 Thus, phase wrapping occurs at positions Reconstructed slices from phase wrapped projections contain cupping artifacts, wherein d values at the edges are lowered with respect to the center, as well as overall reductions in d.This is due to the wrapped edges with jwj < j/j and the integration step during reconstruction.For phase wrapped data, a quantitative measurement of d values is impossible and visualization is difficult because cupping is often larger than d variations between features.
Our proposed correction algorithm employs a priori knowledge of the sample geometry and is applied to the sinogram of the interference pattern phase shift (phase sinogram) before reconstruction.The absorption sinogram is used for identification of the sample edges.The correction is accomplished by the following procedure: Step 1: Locate the edges of the sample and the extent of phase wrapping.In the case of noisy data or slightly noncylindrical specimens, the edges can be fit with a sinusoidal function to avoid jagged edges.The edges determine the center position and the radius of the modeled cylinder, and together with the extent of phase wrapping they define a window where the phase sinogram will later be replaced.This window is extended by a buffer of pixels to account for jagged, noisy edges.This procedure is shown in Fig. 1(a), where the red dashed line indicates the edges and the yellow lines define the replacement window.
Step 2: Create a modeled phase sinogram for a uniform, cylindrical sample with radius R and center x 0 (h) determined from Step 1.This is accomplished using the projected thickness tðx; hÞ ¼ 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi R 2 À ðx À x 0 ðhÞÞ 2 q .A discrete derivative D x replaces the derivative in Eq. ( 1) and a preliminary guess for d gives Here, we have used the assumption that for a uniform sample, Ð dðx 0 ; y 0 Þdy ¼ d Â tðx; hÞ.A final d is refined in Step 5.
Step 3: Replace the phase sinogram with the modeled phase sinogram (from Step 2) within the window (defined in Step 1), as demonstrated in Fig. 1(c).
Step 4: Reconstruct the corrected slice from the sinogram in Step 3.
Step 5: Repeat Steps 2-4 and monitor the flatness of a uniform region of interest in the slice.The optimal model d value is the one which results in the flattest region of interest for a given replacement window, meaning that cupping has been completely removed.Figure 1 shows Steps 1 to 3 applied to sinograms from the paraffin-embedded cerebellum measurement.The correction has been implemented in MATLAB (The MathWorks, Inc., Natick, Massachusetts, United States) and made available in the supplementary material as well as at github.com/grodg-ers1/phase_wrapping_correctionalong with sample data.
To confirm the quantitative accuracy of the correction, two specimens were measured with an double-grating XGI setup both with and without a water bath, allowing for a comparison between a not phase wrapped (NW) dataset, a phase wrapped (W) dataset, and the wrapped dataset corrected with our algorithm (WC).The first specimen consisted of a stack of cylindrical, plastic phantoms including two varieties of polyvinyl chloride (PVC1, PVC2), polyoxymethylene (POM), and polypropylene (PP) with 1 mm in height and 2.5 mm in diameter.To demonstrate the correction in the context of a biological specimen, a cylindrical punch of a paraffin-embedded human cerebellum with a diameter of 6 mm was used.The brain was extracted postmortem and embedded in paraffin following standard histological preparation.All associated procedures were conducted in accordance with the Declaration of Helsinki and were approved by the ethics committee of the Medical School of the National and Kapodistrian University of Athens.
The measurements were performed at the Diamond Manchester Imaging Branchline (I13-2, Diamond Light Source, UK). 15 A double crystal monochromator was used to select a photon energy of 20 keV and images were recorded on a scintillator-coupled detector with an effective pixel size of 4.6 lm and a field of view of 9.0 Â 6.0 mm 2 .The grating interferometer consisted of an absorption grating as the beam splitter grating G1 with period p 1 ¼ 7 lm.The inter-grating distance was set to the first Talbot order of 0.8 m and the analyzer grating G2 had the same period as G1.Our setup achieved a mean visibility of 35% over a 7.6 Â 6.0 mm 2 field of view.
The phantom and cerebellum datasets were taken with 501 and 1201 projections, respectively, with 5 phase steps per projection around 180 .The exposure time was selected to maintain nearly equal counts within the specimen for projections of the same sample with and without water bath.
Figure 2 shows the measured mean d values and associated standard deviations of the phantoms in the WC (blue) and W (green) datasets against the values from the NW dataset.All d values are relative to PVC2 because the surrounding medium was different (water for NW and air for W and WC).In Fig. 2, the red line with a slope of one is a visual aid indicating perfect agreement with the NW dataset.
The paraffin-embedded cerebellum is visualized in Figs.3(a  artifacts in the W dataset (green), as well as the strong agreement between the NW and WC datasets (red and blue, respectively).The histograms for 25 masked slices are shown in (e), with the W dataset smeared out due to cupping and a strong resemblance between the NW and WC datasets.These experimental results confirm the quantitative and qualitative accuracy of the proposed correction algorithm.
The measured d values of the phantoms in the WC dataset match those of the NW dataset, while the W dataset produces inaccurate values with high variations due to cupping.This cupping is clearly visible in the paraffin-embedded cerebellum dataset, but has been removed by the correction algorithm, as demonstrated by slices and line profiles.The histograms of the three datasets demonstrate the improved density resolution after the application of our correction algorithm.
In conclusion, we have presented an algorithm to remove phase wrapping artifacts from phase tomography.This algorithm uses the absorption sinogram to determine the specimen geometry, then models and replaces the phase sinogram in the wrapped regions.The selection of d for the modeled sinogram can be determined automatically by minimizing cupping in the resulting reconstruction.Finally, we present experimental results from phase tomography datasets with phase wrapping, with phase wrapping and corrected with the proposed algorithm, and without phase wrapping.These results confirm that datasets corrected with our proposed algorithm are quantitatively and qualitatively consistent with datasets that had no phase wrapping.Our correction algorithm allows for measurements without submerging the sample in a water bath, improving the experimental setup and avoiding motion artifacts from gas bubbles.
While we used this algorithm with phase-stepping data from a double-grating setup, its application is not restricted to these modalities; the method could also be used, e.g., for single-grating setups or for data acquired in moir e mode.
See supplementary material for the Matlab source code of the proposed correction algorithm.It is also available along with sample data and examples at github.com/grodg-ers1/phase_wrapping_correction.

FIG. 1 .
FIG. 1. Proposed correction algorithm applied to cerebellum sinograms.Edges (red) are found from the absorption sinogram (a) and a replacement window (yellow) is defined, then transferred to the phase sinogram (b) where phase wrapping is present.(c) A model sinogram is created using this geometry and replaces the phase wrapped data within the window.The grayscale for (a) corresponds to the transmission values of [0.6, 1.2] and for (b) and (c) the interference pattern shift of [Àp, p].The vertical scale bar represents 250 lm.

FIG. 3 .
FIG.2.Decrement of the index of refraction Dd of each phantom (relative to PVC2) as measured in the wrapped (W) (green) and wrapped and corrected (WC) (blue) datasets, plotted against the dataset without phase wrapping (NW) (x-axis).The red dotted line of slope ¼ 1 indicates agreement with the NW dataset.Error bars correspond to the standard deviation of measured d values.