Near-field speckle-scanning-based x-ray imaging

The x-ray near-field speckle-scanning concept is an approach recently introduced to obtain absorption, phase, and dark-field images of a sample. In this paper, we present ways of recovering from a sample its ultrasmall-angle x-ray scattering distribution using numerical deconvolution. We also show how to access the 2D phase gradient signal from random step scans, the latter having the potential to elude the flat-field correction error. Each feature is explained theoretically and demonstrated experimentally at a synchrotron x-ray facility.


I. INTRODUCTION
When x-ray imaging was discovered, it took only a few months for the scientific community to express a great interest, as acknowledged by the award of the first Nobel prize to Röntgen. Since then fast and tremendous progress has been observed in diverse applications over domains including medicine, materials science, paleontology, and many other fields. For more than a century great minds have realized the importance of x-ray imaging and contributed to its development with, for instance, the introduction of the dose concept [1,2] and the advent of computed tomography [3,4]. In parallel, constant progress in instrumentation has allowed the development of x-ray imaging systems with ever better contrast, resolution, and efficiency [5][6][7].
While questing for means of collecting more information and reducing the dose necessary to image a specimen, the appearance of the first coherent x-ray sources favored the emergence of phase-contrast imaging methods [6,8,9]. Indeed, these were motivated by the fact that, in the x-ray domain, the material optical index n = 1 − δ − iβ presents usually a refractive index δ with a magnitude several orders larger than its absorption counterpart β [10]. Thus, since the inception of coherent third-generation synchrotron sources, different near-field techniques of x-ray phase sensing emerged, either inspired by their analogs in other spectral domains or specially developed [11]. One may categorize those imaging techniques depending on whether they are sensitive to the Laplacian or the gradient of the phase. The propagation-based methods, which belong to the first category, make use of the property of contrast enhancement of edges when increasing the distance to the detector from the specimen under study [12][13][14][15]. The second category brings together instruments such as deflection-based instruments, including x-ray grating interferometers [16,17], Hartmann-like sensors [18,19], coded apertures [20,21], analyzer-based systems [22], and specklebased methods [23,24]. Other full-field imaging approaches available in the x-ray domain and using different principles * berujon@esrf.eu include for instance ptychography [25] and Zernike phase microscopy [26,27].
Propagation-based methods offer the advantage of a high sensitivity to sharp changes in electronic density, the sample high-frequency features being rendered with high contrast thanks to the phenomenon of optical interference. Conversely, they often suffer from reconstruction artifacts when imaging materials presenting a slowly varying density and low spatial frequencies due to the very nature of the recorded signal, i.e., a Fresnel diffraction image [28]. In practice, a phasegradient-sensitive method would be preferable for imaging a homogenous material, although this usually results in a spatial resolution lower than the one obtained with a propagationbased technique, all other conditions remaining equal. For their implementation, phase-gradient-sensitive methods require the use of a specific optics (grid, grating, crystal, etc.) to modulate the wave front, which can, sometimes, be expensive or difficult to get.
Dark-field x-ray contrast is an additional contrast mode inspired from its counterpart in the visible light and made possible thanks to the availability of crystal analyzers [22,29] and of grating interferometers [30]. Dark-field imaging permits revealing the structure and the orientation of features having a characteristic size much smaller than the detector resolution. For instance, the mapping of the dark-field signal with a grating interferometer allows one to determine the local reduction of the interference amplitude caused by x-ray beam propagation through the sample. The interference amplitude actually reflects the progressive degradation of the beam transverse coherence properties caused by scattering induced upon photon propagation through the sample's submicroscopic features.
The x-ray near-field speckle-scanning technique [31,32] is a recent and particularly attractive deflection-based technique as it provides both phase and dark-field contrast images in addition to the traditional absorption one. Moreover, it uses only a random phase object as wave front modulator. This phase object is used to generate a speckle intensity pattern in the near-field domain caused by mutual interference of the scattered waves. The near-field domain is particularly interesting because the speckle pattern distortion is only ruled by the wave front transformation upon propagation in the Fresnel zone [33,34], a region easily accessible with x rays [35].
The x-ray near-field speckle-tracking technique previously developed and easily implemented [23,24] could provide phase information from one single exposure of the sample.
However, its resolution is limited by the speckle grain size and artifacts may appear in the case of complex objects. The speckle-scanning method appeared to overcome these limitations, nonetheless at the cost of additional image acquisitions. The technique advantages include a high sensitivity to the phase gradient, a moderate requirement for the beam monochromaticity, a high photon efficiency (no absorption in the system), and the use of a minimalist optical element in the form of a scattering membrane. Albeit the speckle-based methods were first developed at a synchrotron source, they were recently adapted with success to a laboratory source, thereby broadening the range of applications [36]. In this paper, several additional aspects of the speckle-scanning concept are presented. After having recalled the basic concept and the general experimental setup, we will expound a method capable of accessing the ultrasmall-angle x-ray scattering distribution of a sample from 2D raster scans. Next, a processing scheme to recover the phase gradient information in two dimensions from scans with random steps is demonstrated. Each aspect is illustrated with experimental applications. We also describe a way of eluding the "flat-field correction error" using the speckle-scanning technique. A short discussion dealing with the sensitivity and operational constraints concludes this paper.

II. EXPERIMENTAL SETUP
Let us introduce the experimental setup that was used to acquire the data presented further. The experiments were conducted at the beamline BM05 of the ESRF [37]. There, the photons are produced by synchrotron radiation from a bending magnet of 0.85 T installed on a storage ring operating with 6.04 GeV electrons. In the experiments presented, a beam with an energy E = 17 keV was selected thanks to a double-crystal Si(111) monochromator with a selectivity of E/E ∼ 10 −4 . The experimental setup sketched in Fig. 1 was installed in a lead-shielded hutch, at a distance from the scattering object to the source of R ≈ 55 m. Considering the source size of 29 μm rms vertically and ∼77 μm rms horizontally, the transverse coherence at the level of the membrane at E = 17 keV is expected to be at most ∼29 μm vertically by 11 μm horizontally.
The scattering object, a piece of sandpaper with a grit designation P800, was mounted on a motorized 2D piezo stage, itself mounted above stepping micromotors. Such a setup permitted moving the random phase object transversally to the beam with either a nanometer precision and a 100 μm range, or a micrometer accuracy over a millimeter range. The scatterer was placed at a distance l = 450 mm upstream of the sample and the detector at a variable distance d downstream from it. The imaging system was a FReLoN (fast readout low noise) CCD camera receiving a visible-light signal resulting from the luminescence produced by x-ray illumination of a scintillator and imaged by a microscope objective. The system effective pixel size was of s p = 5.8 μm. Near-field speckle of a size of a few pixels was observable in the recorded images with a contrast defined by a coefficient of variation of ∼14% when d = 850 mm [visible in Figs. 1(a) and 1(b)].

III. SCANNING TECHNIQUE WITH CONSTANT STEP
This speckle-scanning mode is mainly presented in Ref. [31] where the concept is described as the crossroads of the speckle tracking [23,24] and grating interferometry techniques [38,39]. Since we are going to present an extension of it, let us recall quickly the method principle. Figure 2 shows the geometry considered with the vectors (e x ,e y ) perpendicular to the beam propagation direction e z . The scattering object with random phasors is placed in a partially coherent x-ray beam, either right upstream or downstream from the investigated object. While small differences exist between these two configurations [24,40] they are negligible in the following experiments as the beam was nearly collimated. Nonetheless, we will rederive here the correct formulas for the magnifying geometries currently encountered, e.g., when using a focusing optics or a highly divergent laboratory x-ray source.
The method consists of coupling together data collected when raster-scanning the membrane with the sample present in the x-ray beam [e.g., array of Fig. 1(c)] or removed from it [ Fig. 1(b)]. For the technique to be effective, the membrane must be scanned with a positional repeatability p = (X,Y ) of a fraction of a speckle grain size and with a strictly constant step size s s = constant. A pair of patterns recorded in a pixel located at r = (x,y) is considered to be independent from the data collected in the neighboring pixel. The pixel marked in yellow in Fig. 2(a) with the pattern pair displayed in panels (b) and (c) is an example of data sets recorded at a given position r. Then the 2D intensity value array recorded when scanning the membrane with the sample in the beam path and denoted f [ Fig. 2(c)] is cross-correlated with the one obtained during the reference scan and denoted g [ Fig. 2(b)], i.e., with no sample in the beam. The outcome of the operation is a correlation map whose maximum peak location indicates the displacement vector v between the two arrays: Well-known recipes based on correlation coefficient curvefitting permit us to calculate v with a fraction of a step accuracy [41,42]. Then, from the angles α, γ , and ζ defined in Fig. 4 and for the two geometries (a) or (b), we have α = γ + ζ and so by calculation which is valid under the small-angle approximation encountered with x rays. For a collimated beam, = (R + d + l)/R reduces to 1 and α = ζ . In the following we use the del operator denoted ∇. As the angle α relates to the differential wave front gradient α = ∇W itself linked to the phase gradient ∇φ [43], we can write where k is the wave number equal to 2π/λ, with λ the photon wavelength. Most often, data recorded through 2D scans are followed by projection of v in order to extract the two orthogonal transverse phase gradient components.
One-dimensional scans with many fewer points have also proved to be efficient for recovering the 1D phase gradient [31]. The absorption, or here the transmission image T , can be correctly calculated free from the flat-field correction errors that occur in the presence of a strongly dephasing object and/or of a high magnification projection geometry [25,44]. Its formulas can be obtained by consideration of Fig. 4: where μ denotes the mean operator over the variable p. Finally, the dark-field image D f can be reconstructed by considering the decrease of speckle contrast that occurs when introducing the sample in the beam. More precisely one has to calculate the ratio of the coefficient of variations of f and g: with σ the standard deviation operator over p.
Using a blueberry as a sample, Figs. 3(a)-3(d) illustrate the various modes the technique offers. The set of two scans was recorded with d = 850 mm, s s = 3 μm, each scan being composed of an array of 25 × 25 points. The collection of such large data scans was performed in an attempt to optimize the calculation of the higher moments of the scattering distribution that require more statistics, as we shall see further. Otherwise many fewer points would be needed. Figure 3 presents the calculated transmission image in panel (a), in panels (b) and (c) the respectively horizontal and vertical wave front gradients ∇W · e x/y , in panel (d) the integrated phase image, and in panel (e) the dark-field image.

IV. ULTRASMALL-ANGLE X-RAY SCATTERING DISTRIBUTIONS
Recent works demonstrated a way of recovering the scattering distribution of a sample in each detector pixel from raster scans. The concept was first introduced by Yashiro et al. in [45] and [46] using a grating interferometer, wherein the scattering effect of the sample is treated as a local point spread function whose estimated width mirrors the local scattering strength. Later on, Modregger et al. [47,48] also used a grating interferometer but, this time, combined with numerical deconvolution to be able to extract the higher moments of the scattering distribution. In this section, we introduce by equivalence the recovery of the ultrasmallangle x-ray scattering distribution using near-field speckle scans similar to the ones recorded and used in the previous section.
Let us consider the light distribution f in the case of a sample present in the x-ray beam path. Upon propagation, and for a given pixel position r, we can consider f as the effect of the sample optical transfer function o(r,p) convoluted with the input reference signal g [31,47]: Henceforth, for each pair of patterns (with and without sample) (f,g)(r), a deconvolution procedure can be applied to retrieve the distribution o. For this purpose the iterative Richardson-Lucy deconvolution [49][50][51] algorithm is convenient because it is robust and its input does not require any assumption on the system noise [52]. Such processing was applied to the blueberry sample data of Sec. III with 100 iterations used in the Richardson-Lucy deconvolution for each (μrad) ∇W · e x (μrad) ∇W · e y (μrad) ∇W · e y (μrad) (rad) ∇W · e x (μrad) pixel. An illustration of the scattering distribution o resulting from the deconvolution of the data of Figs. 2(b) and 2(c) is displayed in Fig. 2(d).
Geometrically, the position p corresponds to the angular coordinate α = s s d p (see Fig. 4). By projection in this angular variable coordinate system, a scattering image o(α) can be constructed for each angular vector. Figure 5 shows an image array of the blueberry sample organized as a function of the scattering angular coordinates. Given a scan step of s s = 3 μm, the angular step size from image to image was of 3.5 μrad.
Moments of the scattering distribution contain physical information about the sample features and correspond, as we shall see for the first of them, to the traditional contrast modalities, i.e., the absorption, phase, and scattering signals. The central moments M mn of order (m,n),m + n > 1 of o(r) are calculated within the meaning of angular probability density functions:  where M 00 = o(r,p)dp (8) and (X r ,Y r ) are the centroid positions of o(r,p): We have then v(r) = (X r ,Y r ) by definition of v, from which we can calculate ∇φ using Eq. (3). Furthermore, the two second normalized moments √ M 02 and √ M 20 mirror the angular characteristic scattering width in the two transverse directions. Higher moments provide complementary information on the scattering distribution, with, for instance, the third and fourth moments related to its skewness and kurtosis. The wave front differential gradients ∇W · e x/y for the blueberry sample obtained using the scattering distribution centroid are shown in Figs

V. SCANS WITH RANDOM STEPS
In this section, we present a processing method to recover the sample-induced 2D phase gradients from images acquired using random-step scans, i.e., scans where the images are recorded with the scattering membrane moved at random transverse positions (X,Y ) in the plan (e x ,e y ), with no particular translation distance from one image to the next. However, it is crucial that scans performed with and without the sample have their images taken with the scattering membrane located at the very same positions for each corresponding step.
Here, we will change slightly the notation, in order not to confuse the reader with the previous processing method. We note i(r,n) the intensity collected at the position r of the nth image of a scan containing N images when the sample is inserted into the beam. Then, writing i r,n = i(r,n), we have the vectors i r = (i r,1 , . . . ,i r,n , . . . ,i r,N ) built from N realizations of the random variables I r , induced by the variable position of the membrane. With equivalent notations, we also build the reference vectors ω r , from N realizations of W r obtained in absence of a sample.
The idea is here again to use correlation calculations between the speckle distribution in the various pixels to retrieve the beam wave front derivative through the calculation of the local light deflection angles α(r) as sketched in Fig. 6. So, we use the Pearson correlation coefficient ρ r generally defined by where E denotes the expected value. The factor ρ is the normalized covariance of the pixel signals distributions. It is used as figure of merit to evaluate the level of similarity between the data acquired in the different pixels of the two scans. For the vectors i and ω, this operation can also be seen as a scalar product. Thus, we calculate an estimate of each projection vector of i r onto vectors ω r+h located at a small distance h from r in the reference scans: where i r and ω r denote the mean values of the measurement vectors. For all i r , and with the correlation map ρ(i r ,ω r+h ) [see Fig. 6(c)], we can locate the vector position of ω r+h max rendering maximum correlation: Here also the vector h max is determined with a subpixel accuracy by numerical interpolation over the neighboring pixels [42,53]. Following the law of light propagation in homogeneous media, the phase gradient ∇φ is recovered by where s p is the detector pixel size. The corrected transmission signal is equal to and is valid for any geometry. With this processing method the dark-field signal is written which is fully equivalent to Eq. (5). This method can be seen as an oversampling version of the x-ray speckle-tracking technique. Indeed, the cross-correlation operation is operated between data collected in different pixels as performed with speckle tracking [23]. Moreover, the refraction angle α is here calculated at the sample position (see Fig. 6) and not at the detector plane level as in the case of the method of Sec. III with the angle ζ . An experimental demonstration of the processing scheme is provided in Fig. 7, using the same blueberry sample as before, for comparison purposes. A pair of 25-image scans was acquired by applying displacements to the membrane with steps larger than 20 μm using a micrometer precision linear  Figure 7(f) displays a plot of cuts along the marked lines of the three previous images normalized to an arbitrary scale. One can observe the way the different contrasts vary depending on the image modality.
This processing scheme is particularly interesting as it allows us to recover the 2D phase and scattering information for every pixel, as does the scheme with constant step, although requiring a lower number of images. With the constant-step method, O(N 2 ) images are necessary to build the 4D matrix f (r,p) while only O(N ) are required in this method for 2D phase sensing.

A. Sensitivity
The sensitivity of the speckle-scanning technique differs depending on the processing method. When applying the algorithm presented in Sec. III on data acquired with a scan of constant step, the accuracy σ (α) on the wave front gradient is This can be compared to the processing method introduced in Sec. V, for instance when using the same data: From these equations, with a similar experimental setup, it is the constant-step scanning method that delivers the highest sensitivity. As a matter of fact, considering an equivalent subpixel/step accuracy for the peak finder algorithm [41] in both processing methods, we have σ (v r ) ≈ σ (h r ) and we can set s s < s p to diminish σ 1 . When working with a collimated x-ray beam, the gain in sensitivity gets rapidly limited by the minimum step length necessary to get a sufficient intensity variation within a pixel from one step to the next [42,54]. On the contrary, with a magnification projection geometry setup, the gain can be tremendous since a tiny displacement of the membrane located near the x-ray beam focal point can generate a much larger displacement of the speckle pattern at the detector level. In this case, when generating a sufficient speckle statistic variation for the correlation algorithm to be efficient and accurate, the angular resolution may reach the nanoradian scale [55].
The convergence of the distribution estimators in Eq. (11) scales with 1/ √ N [56]. This implies that only the few first images are necessary to reach a suitable accuracy since √ Nσ (h) → 0 and √ Nσ (v) → 0 as N → ∞. From Eq. (10), considering two identical scans or areas of r where the x-ray light is neither refracted nor scattered by any sample, the distributions (I r ,W r ) should, by definition, render a theoretical value ρ(I r ,W r ) = 1. For our experimental data, when using the method of Sec. V, the deviation from unity of ρ over an area of 150×150 pixels amounts to σ (α) = 0.8 μrad for a scan of 25 images. When more than 150 images are used it tends to a limit of σ (α) = 0.35 μrad. This last figure is comparable to the one calculated with the two constant-step methods presented in Secs. III and IV, suggesting that the discrepancy may be due to monochromator instabilities, more felt in the vertical direction.

B. Tracking orientation
Section V explains and illustrates the way the refraction angle is recovered by examining the vector positions in the original reference scan data. This subsection aims at illustrating the importance of operating along this line rather than tracking the reference scan vectors in the data array collected when the sample is present in the beam. The problem raised here is due to the creation of optical vortices in the propagated beam when going through the sample discontinuities [57]. Figure 8 presents an experimental data set for a juniper berry sample imaged with x rays. That sample exhibits strong scattering features and turbid phase shift. Figures 8(a) and 8(b) successfully display the transmission and phase contrast images, respectively. Figure 8(c) presents the horizontal speckle displacement vectors calculated picking vectors within the sample scan and searching for their counterpart vectors in the reference scan. Conversely, Fig. 8(d) corresponds to the tracking of vectors from the reference scan across the sample scan data. Figures 8(e) and 8(f) zoom-in the square subsets of Figs. 8(c) and 8(d), while Fig. 8(g) shows the absolute intensity difference between these two subsets. The significant difference observed between the two calculated images results from the incapacity of the algorithm to seize reference vectors folded up in optical vortices upon propagation. As these vectors add up at singularity points they generate new distributions I r [57] that are no longer correlated to any reference vector. Obviously, as the detector pixel size and the speckle grains themselves get smaller, or as the propagation distance increases, the effect becomes more pronounced.

VII. CONCLUSION
Several advanced processing schemes were demonstrated within the context of our x-ray speckle-scanning imaging technique. Starting from the presently available scanning method with a constant step, we proposed ways of building error-free transmission images and of accessing the sample scattering distribution, pixel by pixel, through a numerical deconvolution of the scan data. The various processing schemes were demonstrated with experimental data, assessing thereby the excellence of the method. The recovery of the sample scattering distribution is of interest in various fields of materials science as well as in biological imaging where the orientations of subpixel structures are responsible for the functional behavior of the systems. In addition, a processing scheme with relaxed constraint on the step size was proposed. It can be seen as an oversampled version of the speckle-tracking technique with a resolution enhanced up to the one defined by the detector. The speckle-scanning imaging technique based on a random-step scheme reduces dramatically the number of acquisitions as compared to the previously available schemes. This point is often essential in the case of biological imaging, where the photon efficiency permits us to reduce the dose delivered to the sample. Future work will permit pushing and testing the limit of the method in this regard. The proposed approach will also promote the development of the speckle imaging technique with laboratory low-brilliance sources.