The effect of the spatial sampling rate on quantitative phase information extracted from planar and tomographic edge illumination x-ray phase contrast images

The effect of the spatial sampling rate on the quantitative phase information that can be retrieved from planar and tomographic edge illumination (EI) x-ray phase contrast imaging (XPCi) with a laboratory-based prototype scanner is analysed. The study is conducted on simulated and experimental data from a custom-built phantom. Optimal sampling rates, i.e. the minimum ones allowing the unambiguous extraction of quantitative phase measurements from the acquired data, are identified for planar and tomographic imaging. One of the key outcomes of this study is the demonstration that the optimal sampling rate in tomographic EI XPCi is low, allowing the implementation of low-dose volumetric imaging without having to compromise on quantitative accuracy.


Introduction
X-ray phase contrast imaging (XPCi) has been the focus of intensive research over the past two decades [1][2][3][4][5]. The most recent developments in XPCi are the transition from specialised synchrotron facilities into standard research laboratories [6,7], and the combination with the principles of computed tomography (CT), enabling 3D imaging [8][9][10][11][12]. The interest in XPCi originates from the fact that improved image contrast can be obtained for weakly attenuating specimen such as biological soft tissue [13]. The contrast in XPCi arises due to phase (refraction) effects, which, for photon energies used in biomedical imaging, can be much larger than the attenuation effects exploited by conventional radiography [13]. Both effects can be described by the complex refractive index: where E is the photon energy and δ and β drive phase (refraction) and attenuation effects, respectively. The edge illumination (EI) XPCi method [14][15][16][17][18][19] can be implemented both with synchrotron and conventional, nonmicrofocal x-ray sources [7], and is therefore well-suited to laboratory based applications. The working principle of EI XPCi is schematically shown in figure 1: an x-ray mask positioned upstream of the sample separates the incoming (cone) beam into individual beamlets, with gaps between them sufficiently wide as to prevent them from interfering. A second x-ray mask creates insensitive areas between the pixels of a detector. When the beamlets are aligned with the edges of these insensitive areas, an edge illumination condition is achieved, making the system sensitive to refraction in addition Journal of Physics D: Applied Physics The effect of the spatial sampling rate on quantitative phase information extracted from planar and tomographic edge illumination x-ray phase contrast images to attenuation effects [7]. In the absence of small angle scattering, the angle of refraction along x, given by: can be extracted when two projections are acquired while illuminating opposite edges on the detector aperture, and processed together by means of a dedicated algorithm [17]. In equation (2), which in the following will be referred to as a 'differential phase projection', k is the wave number and Φ is the phase shift given by: d.
x y s ℓ( , ; ) The line ℓ(x, y;s) describes the trajectory of an x-ray hitting the detector plane at position (x, y). The phase shift can be recovered from equation (2) via a 1D integration along the transverse (x-) detector direction. The constant of integration can be fixed if a region outside the object exists where δ is constant and known.
When EI XPCi is operated in CT mode, equation (2) becomes a differential phase sinogram that shows the first derivative of the Radon transform (R) of the refractive index decrement δ: where θ is the angle by which the sample is rotated around an axis aligned with the y-direction. From equation (4) a tomographic map of δ can be reconstructed, e.g. via the filtered back-projection (FBP) formula after integration of S along the transverse direction. Equivalently, a reconstruction of δ can be obtained without prior integration of S when the ramp filter in the FBP is replaced by the Hilbert filter, which performs the integration implicitly in Fourier space [20].
In both the planar and CT cases, the retrieval of quantitative phase information, i.e. of the phase shift Φ (planar) and of the refractive index decrement δ (CT), requires a one-dimensional integration. This integration is generally regarded as problematic, as the propagation of local measurement errors can cause artefacts [21]. In addition to measurement errors, the integrated differential phase projections and sinograms can also be affected by the spatial sampling rate. Differential data typically contain high frequency components, which, according to the Nyquist-Shannon theorem, need to be sampled at a high rate to be fully captured [22,23]. If sampled at too low a rate, an artefact-free integration is impossible. Not only do the resulting artefacts corrupt the integrated differential phase projections and sinograms qualitatively, but they also may affect their quantitative meaning.
The spatial sampling rate of a differential phase projection obtained from EI XPCi data is determined by the demagnified detector pixel size (i.e. the centre-to-centre distance by which the beamlets are separated); however, in the transverse direction it can be increased by dithering, i.e. shifting the sample by sub-pixel amounts and taking multiple projections which are then combined into a single, upsampled projection. The dithering direction is indicated by the dotted arrow in figure 1. While effectively improving the spatial resolution, dithering comes at the cost of increased scan time and dose by a factor equal to the number of dithering steps.
In a recent publication [24], we demonstrated that, when operated in CT mode, EI XPCi can provide reconstructions of the refractive index decrement δ which are accurate within the polychromaticity constraints imposed by the broad energy spectrum of a conventional x-ray source [25]. One of the main observations was that accurate results were obtained without dithering, i.e. with a relatively low sampling rate and, hence, with a low dose. In this article, the effect of the sampling rate on the quantitative phase information that can be recovered from planar and tomographic EI XPCi measurements acquired with a specific laboratory based setup is analysed for a custom-built wire phantom. We show that, in the tomographic case, the sampling rate and the extracted quantitative information are virtually independent, which confirms the observation reported in our previous publication [24]. This is an important result with respect to in vivo applications (e.g. small animal imaging), as it implies the possibility to perform CT imaging at low doses while remaining quantitative.

Materials and methods
The phantom consisted of three wires (polyether ether ketone (PEEK), diameter = 450 µm, sapphire, diameter = 250 µm, aluminium, diameter = 250 µm) oriented parallel to the apertures, which were held in place by a plastic cylinder (diameter = 21 mm, wall thickness = 1.5 mm). To assess the effect of the sampling rate on the extractable quantitative information, a combination of experimental and simulated EI XPCi data of the phantom was considered.
The experiments were performed using an EI XPCi scanner prototype with a source-to-detector distance of 2 m. The Anrad SMAM a-Se flat panel detector (Analogic Canada Corporation, Canada) was employed, which has a pixel size of 85 µm. The two x-ray masks, fabricated by electroplating gold strips onto a graphite substrate (Creatv Microtech Inc. Potomac, MD), were located at 1.6 m and 1.96 m from the source, and their periods were 66.8 µm and 83.5 µm respectively, matching the demagnified pixel size. The masks' slits were 12 µm and 20 µm wide, respectively. The edge illumination condition was achieved by shifting the first mask such that 50% of each beamlet fell onto the insensitive detector areas created by the second mask. The sample was positioned immediately downstream of the first mask, at approximately 1.65 m from the source (z so ). The object-to-detector distance (z od ) was therefore 35 cm. The x-ray source was the Rigaku MicroMax 007 HF (Rigaku Corporation, Japan) with a rotating molybdenum target, operated at 25 mA/40 kVp (planar acquisitions) and 25 mA/35 kVp (CT). The source focal spot (σ) was measured to be 70 µm. In each experiment, the exposure time was 5 s per projection and dithering step.
The simulated data were obtained using an EI XPCi simulation code based on a rigorous wave optics approach [26]. The input parameters for the simulation matched the experimental ones; however, while the experiments were performed with a polychromatic, unfiltered molybdenum spectrum, for the simulation a photon energy of 18 keV (the source's mean energy) was assumed. A numerical wire phantom was created that closely matches the real one.
The following data acquisition and analysis strategy was applied to both simulated and experimental data. Planar and tomographic EI XPCi data were initially obtained with a very high sampling rate (128 and 32 points/pixel). From these densely sampled data, quantitative phase information was extracted: the phase shift Φ and the refractive index decrement δ of the phantom were extracted in the planar and tomographic case, respectively. The values of Φ and δ extracted in this way were considered the 'reference', as the used sampling steps (approximately 0.5 µm and 2 µm) were much smaller than the projected source focal spot size demagnified to the sample plane (σ· z od /(z so + z od ) ≈ 14 µm), justifying the assumption that these values are unaffected by under-sampling artefacts. These densely sampled data were then consecutively sub-sampled by factors of two (i.e. if the sampling rate in the initial data was 128 points/pixel, the sampling rates in the sub-sampled data were 64, 32, 16 points/pixel and so on) and Φ and δ were extracted again. This sub-sampling continued until a sampling rate of 1 point/ pixel was reached. As illustrated in figure 2, the sub-sampling factor determines how many sub-sampled datasets exist (e.g. sub-sampling by a factor of two produces two choices for the sub-sampled dataset, and so on). In practice, these different datasets correspond to different locations of the sampling points, determined by the relative position of the sample and the beamlets. All these sub-sampled datasets were taken into account, i.e. Φ and δ were extracted for each one of them. The minimum, mean and maximum values of Φ and δ that were extracted for a given sub-sampling factor were plotted against the corresponding sampling rate. Finally, the sampling rate for which the extracted minimum, maximum and mean Φand δ-values overlapped was considered optimal, in the sense that no further improvement is obtained by increasing the sampling rate. This also means that, for the optimal sampling rate, unambiguous Φ-and δ-values are extracted.

Generation of planar EI XPCi data
Projections of the wire phantom were simulated/experimentally acquired while illuminating opposite edges of the apertures in the detector mask, with an initial sampling rate of 128 points/pixel. From these projections, differential phase data (equation (2)) were extracted [17], and integrated using a simple sum rule in order to obtain images of the phase shift Φ: The coordinates (i, j) are a discretisation of the coordinates (x, y) and Δs is the sampling interval that corresponds to a given sampling rate. The value of the phase shift image at the centre of each wire in the phantom was extracted. The finely sampled data were then consecutively sub-sampled as described above, and the same procedure was repeated. Due to the polychromaticity of the source, the phase shifts extracted from experimental data refer to effective energies. These depend on the imaging system as well as on the attenuation properties of the wires themselves, and therefore their estimation is not straightforward [25]. As a consequence, the multiplication by the wave number k was omitted when the phase shift images were calculated from the differential phase projections. For the sake of consistency, the same approach was used with the monochromatic simulated data. Omitting the factor k has no influence on the effect of the sampling rate, and therefore on the conclusions of our study.

Generation of tomographic EI XPCi data
Projections were simulated/experimentally acquired while illuminating opposite aperture edges on the detector. The tomographic scan involved rotating the sample over a total angular range of 180 degrees, with an angular increment of 1 degree. The axis of rotation was aligned with the y-direction. The initial sampling rate for each projection was 32 points/pixel. This lower initial sampling rate compared to the planar case was chosen due to constraints on computation and scan time. Sinograms in the form of equation (4) were computed from the set of projections [17], and integrated along the transverse direction using a simple sum rule (equation (5)). Tomographic maps of δ were reconstructed using FBP with the ramp filter, and average values of δ were extracted for each wire. As in the planar case, the data were then repeatedly sub-sampled, and the reconstruction process and extraction of δ was repeated.

Results for planar EI XPCi
The experimentally acquired/simulated planar EI XPCi data and the results of the analysis are shown in figure 3. The experimental data (top row), which were acquired with a polychromatic source and, hence, correspond to effective energies, are indicated by hat symbols α Φk ( , , ) in the vertical axis labels. The simulated data are shown in the bottom row. The plots in figures 3(a) and (d) show profiles extracted from an individual differential phase projection, sampled at a rate of 128 points/ pixel. The bright and dark fringes visible at interfaces within the sample, which can be clearly appreciated in these profiles, are a characteristic feature of differential phase projections. The strong fringes at the edges of the profile are due to the cylinder supporting the wires, and the three fringe pairs at the centre originate from the three wires (the materials of which are explicitly indicated in panel figure 3(b), which shows the wires in the same order). A comparison of the profiles in figures 3(a) and (d) reveals a good match between experiment and simulation, the only difference being a slightly weaker signal in the experimental profile. This is explained by the fact that the effective energies for the wire materials lie above the source's mean energy of 18 keV [25]. The plots in figures 3(b) and (e) show the phase shift obtained by integrating the differential phase profiles. Again a good match is observed when comparing experiment and simulation.
The plots in figures 3(c) and (f) show the effect of the sampling rate on the phase shift values extracted for the three wires. The error bars in the experimental analysis plot (figure 3(c)) represent the standard deviation of the background noise in the phase shift image. The three curves for each material describe the minimum, mean and maximum phase shift values that were extracted for each sub-sampling of the initial data. All curves show the same behaviour for experimental and simulated data: for a low sampling rate, the minimum and maximum values diverge significantly from the mean ones, but progressively converge towards them when the sampling rate is increased. The large discrepancy between maximum and minimum values implies a large potential error in the retrieved phase shift when a low sampling rate is used: for example, for a sampling rate of 1 point/pixel (no dithering) the error can be on the order of a hundred per cent, making the extraction of accurate quantitative phase information impossible. Moreover, negative values can also appear, which has no physical meaning. These occur for example when the positive peaks of the cylinder in the differential phase profiles are more severely affected by under-sampling artefacts than the negative ones. Although the wires have different refractive index decrements, the curves describing their mean extracted phase shift values are close to each other. This is not contradictory considering that the phase shifts of the wires are additionally determined by their diameters, inclination, rotation and location within the cylinder.
In the experimental case, the optimal sampling rate (i.e. the one at which minimum, maximum and mean values for all wires converge, and it is therefore possible to identify the phase shift univocally) appears to be 32 points/pixel (see arrow in figure 3(c)). In the simulated analysis plot (figure 3(f )), this can already be observed for a sampling rate of 16 points/pixel. We believe that the slightly higher sampling rate in the experimental case is due to noise and other possible experimental errors that could have affected the analysis.
In a recent publication [27], it was shown that, when an extended laboratory based source is used, the spatial resolution in a differential phase projection acquired with EI XPCi is determined by the smallest between projected source size rescaled to the object plane and width of the pre-sample mask slits. In the setup used here, the sample mask slits (12 µm) are smaller than the downscaled projected source size (approximately 14 µm), and, therefore, determine the spatial resolution. This can also be understood intuitively: areas of the object which are not illuminated by the beamlets cannot contribute to the measured signal. According to the Nyquist-Shannon theorem [22,23], a band-limited signal must be sampled at least at twice the limiting frequency for it to be fully captured (Nyquist rate). The Nyquist rate for the used setup would correspond to a sampling interval of 6 µm. Considering that only sampling rates which are powers of two were considered in the analysis leading to figure 3(f), an optimal sampling rate of 16 points/pixel (corresponding to a sampling interval of approximately 4 µm) is in agreement with Nyquist's prediction. The next lower sampling rate of 8 points/pixel (corresponding to a sampling interval of approximately 8 µm) would have been an under-estimation.  In summary, the plots in figures 3(c) and (f) show that, in planar EI XPCi, spatial resolution and the reliability of the retrieved phase shift values are coupled, in the sense that a high sampling rate is required to obtain unambiguous results. It should be noted that the optimal sampling rate as stated above applies to the used scanner prototype. For other experimental setups the optimal sampling rate could be different, as the spatial resolution depends on the source blurring and the mask design. For example, if the source blurring and the sample mask slits were larger, a lower sampling rate would be required. Likewise, a lower sampling rate could be sufficient if the sample structures are larger, as in this case under-sampling artefacts would not affect the quantitative phase shift to the same degree of severity. It is also noteworthy that the link between the sampling rate and the extractable quantitative information applies to the phase shift images, but not to differential phase projections, i.e. to images showing the refraction angle distribution. In this case, as no integration is needed, increasing the sampling rate is only a means to improve the spatial resolution. This is especially important with regards to radiation dosimetry: while a precise quantitative retrieval of the phase shift using planar EI XPCi may require an increased dose, non-integrated differential phase projections and phase-enhanced projections (i.e. those acquired in only one of the edge illumination configurations, hence containing a mix of phase and attenuation contrast) can be acquired at low doses, as demonstrated in [19].

Results for tomographic EI XPCi
The results of the tomographic EI XPCi data acquisition/simulation are shown in figure 4. The images in figures 4(a) and (c) show slices of the refractive index decrement δ, reconstructed from data acquired with a spatial sampling rate of 32 points/ pixel. Please note again that the grey scale in the experimental slice (figure 4(a)) refers to δ-values at effective energies, whereas the scale in the simulated slice ( figure 4(b)) refers to δ-values at 18 keV. The plots in figures 4(b) and (d) show the effect of the sampling rate on the retrieved δ-values for the three wires. The error bars in the experimental curves correspond to one standard deviation of all pixels that were averaged to obtain the results. The plots show that, although the lines describing minimum, maximum and mean retrieved δ-values diverge for a low spatial sampling rate, the discrepancy is significantly smaller than in the planar case (figures 3(c) and (f )). In fact, the worst possible error on δ that can occur for the lowest sampling rate of 1 point/pixel (no dithering) is less than 10% for all three wires. Moreover, the three lines converge very quickly, such that already for a sampling rate of 2 points/pixel a good agreement among minimum, maximum and mean retrieved δ-values is observed. Again, the error intervals for this rate are slightly larger in the experimental than in the simulated plot, which we believe is due to noise and other experimental errors affecting the analysis.
These results demonstrate that, in the CT case, δ can be extracted unambiguously already when a sampling rate of 2 points/pixel, corresponding to a sampling interval of approximately 33 µm, is used. In this sense, a sampling rate of 2 points/pixel can be defined optimal. If small errors on δ (on the order of some per cent) can be tolerated, a sampling rate of 1 point/pixel can be considered sufficient. The latter can be achieved without any dithering. It should be noted again that the optimal sampling rate refers to the used experimental setup and may vary for a different setup. The lower optimal sampling rate in tomographic compared to planar EI XPCi can be understood intuitively. In planar EI XPCi, a low sampling rate is not sufficient since not all parts of the sample are illuminated by the beamlets, and, therefore, do not contribute to the differential phase image. When the sampling rate is increased via dithering, the sample is shifted through the beamlets until it has been entirely illuminated. This means that all parts of the sample contribute to the differential phase image, the integration of which yields artefact free phase shift images. During a tomographic acquisition, the sample moves along a circular trajectory. As a result, sample features which are not illuminated by the beamlets at one rotation angle are illuminated at another. Therefore, although no dithering is performed, all parts of the sample contribute to the image formation, in other words, the sample rotation effectively replaces the effect of dithering, and the reconstructed CT slice is unaffected by under-sampling artefacts. Figure 5 illustrates this concept by highlighting the trajectory of a specific sample feature during planar (dithered) and tomographic (undithered) acquisitions for a single aperture pair.
At the same time, as in planar EI XPCi, increasing the sampling rate during tomographic acquisitions can be used to improve the spatial resolution in the reconstructed δ-maps [27]. However, our results show that, unlike in the planar case, in CT the spatial resolution and the accuracy of the extracted quantitative information are not coupled. This is confirmed by figure 6, which shows tomographic δ-maps reconstructed from data sampled with a rate of 32 points/pixel (a) and with the optimal rate of 2 points/pixel (b). The zoom of the PEEK wire reveals that the spatial resolution in (a) is higher than in (b); however, from the point of view of the quantitative retrieval of the δ-values, the two maps are equivalent. As an example for a more complicated object, figure 7 shows the same outcome (quantitatively identical δ-maps, reconstructed from data sampled at rates of 32 and 2 points/pixel) for a transverse crosssection of a domestic wasp. The insect was imaged under the same conditions as the wire phantom.

Conclusion
By analysing simulated and experimental images of a custom-built phantom, it was shown that, in planar EI XPCi, a high spatial sampling rate is needed to obtain unambiguous retrieval of the phase shift value (i.e. quantitatively correct integrated differential phase projections). Typically, this corresponds to a sampling step equal to half the maximum spatial resolution provided by the imaging system, equal to the smallest between slit width in pre-sample mask and projected source size downscaled to the sample plane [27]. On the contrary, when EI XPCi is performed in tomographic mode, a much lower spatial sampling rate of 1-2 points/ pixel was proven to be sufficient to reconstruct slices of the refractive index decrement δ, the quantitative value of which is practically unaffected by under-sampling artefacts. This implies that, in tomographic EI XPCi, an increased sampling rate can be used to enhance the spatial resolution, but is not essential to the quantitative reliability of the extracted results. As a consequence, object-specific quantitative imaging can be performed: if the sample does not require a very high resolution, a CT scan without dithering is sufficient, while for samples with smaller features the number of dithering steps can be adjusted depending on the size of the features to be resolved. This approach enables optimal dose control, which is crucial for many biomedical applications, especially in vivo. The optimal sampling rate depends on the imaging setup and, in this study, has been determined for a setup employing an extended laboratory source. The interplay between sampling rate and retrievable quantitative phase information when using coherent synchrotron sources, as well as the effect of continuous scanning schemes, will be part of future investigations.