Single-shot x-ray speckle-based imaging of a single-material object

TL;DR

This paper introduces a single-shot x-ray speckle-based imaging method for measuring the projected thickness of single-material objects using a novel approach that requires only one image and leverages speckle deformation and the transport-of-intensity equation.

Contribution

It presents a new single-shot speckle-based phase imaging technique that simplifies the process by avoiding explicit speckle tracking and uses a closed-form solution for inverse problem solving.

Findings

Method is robust and efficient with noise

Applicable to x-ray phase--amplitude radiography and tomography

Suitable for dynamic and radiation-sensitive samples

Abstract

We develop a means for speckle-based phase imaging of the projected thickness of a single-material object, under the assumption of illumination by spatially random time-independent x-ray speckles. These speckles are generated by passing x rays through a suitable spatially random mask. The method makes use of a single image obtained in the presence of the object, which serves to deform the illuminating speckle field relative to a reference speckle field (which only needs to be measured once) obtained in the presence of the mask and the absence of the object. The method implicitly rather than explicitly tracks speckles, and utilizes the transport-of-intensity equation to give a closed-form solution to the inverse problem of determining the complex transmission function of the object. Implementation using x-ray synchrotron data shows the method to be robust and efficient with respect to noise. Applications include x-ray phase--amplitude radiography and tomography, as well as time-dependent imaging of dynamic and radiation-sensitive samples using low-flux sources.

Equations

$$\frac{\partial I_{z_1}(x,y)}{\partial z}=-\frac{1}{k}{\nabla }_{\perp }\cdot [I_{z_1}(x,y){\nabla }_{\perp }\phi (x,y)],$$

$$\phi (x,y)={\phi }_{\text{ob}}(x,y)+{\phi }_{\text{m}}(x,y).$$

$$I_{z_1+z_2}(x,y)=I_{z_1}(x,y)-\frac{z_2}{k}{\nabla }_{\perp }\cdot [I_{z_1}(x,y){\nabla }_{\perp }\phi (x,y)].$$

$$I_{z_1+z_2}({\mathbf{r}}_{\perp })=\left[I_{z_1}({\mathbf{r}}_{\perp })-\frac{z_2}{k}{\nabla }_{\perp }{I}_{z_1}({\mathbf{r}}_{\perp })\cdot {\nabla }_{\perp }\phi ({\mathbf{r}}_{\perp })\right]$$

$$-\frac{z_2}{k}{I}_{z_1}({\mathbf{r}}_{\perp }){\nabla }_{\perp }^2\phi ({\mathbf{r}}_{\perp })$$

$$\approx I_{z_1}\left({\mathbf{r}}_{\perp }-\frac{z_2}{k}{\nabla }_{\perp }\phi ({\mathbf{r}}_{\perp })\right)-\frac{z_2}{k}{I}_{z_1}({\mathbf{r}}_{\perp }){\nabla }_{\perp }^2\phi ({\mathbf{r}}_{\perp }).$$

$$I_{z_1+z_2}({\mathbf{r}}_{\perp })$$

$$I_{z_1}({\mathbf{r}}_{\perp })$$

$${\tilde{I}}_{z_1+z_2}({\mathbf{r}}_{\perp })$$

$${\tilde{I}}_{\text{ob}}({\mathbf{r}}_{\perp })$$

$$\frac{I_{z_1+z_2}\left({\mathbf{r}}_{\perp }\right)}{{\tilde{I}}_{z_1+z_2}({\mathbf{r}}_{\perp })}\approx I_{\text{ob}}\left({\mathbf{r}}_{\perp }\right)-\frac{z_2}{k}{\nabla }_{\perp }\cdot [I_{\text{ob}}\left({\mathbf{r}}_{\perp }\right){\nabla }_{\perp }{\phi }_{\text{ob}}\left({\mathbf{r}}_{\perp }\right)].$$

$$n=1-\Delta +i\beta =1-\gamma \beta +i\beta =1+\beta (i-\gamma ),$$

$$\gamma =\frac{\Delta }{\beta }.$$

$${\phi }_{\text{ob}}\left({\mathbf{r}}_{\perp }\right)=-k\Delta t\left({\mathbf{r}}_{\perp }\right)=-k\gamma \beta t\left({\mathbf{r}}_{\perp }\right)$$

$$I_{ob}\left({\mathbf{r}}_{\perp }\right)=\exp [-2k\beta t\left({\mathbf{r}}_{\perp }\right)].$$

$$\frac{I_{z_1+z_2}\left({\mathbf{r}}_{\perp }\right)}{{\tilde{I}}_{z_1+z_2}({\mathbf{r}}_{\perp })}$$

$$\frac{I_{z_1+z_2}\left({\mathbf{r}}_{\perp }\right)}{{\tilde{I}}_{z_1+z_2}({\mathbf{r}}_{\perp })}-1=\left[1-\frac{\gamma z_2}{2k}{\nabla }_{\perp }^2\right]\left[-2k\beta t({\mathbf{r}}_{\perp })\right].$$

$$t(x,y)=$$

$$\frac{1}{\mu }{\mathcal{F}}^{-1}\left(\text{LFF}\left\{\frac{\mathcal{F}\left[1-I_{z_1+z_2}(x,y)/{\tilde{I}}_{z_1+z_2}(x,y)\right]}{1+\pi \gamma z_2\lambda (u^2+v^2)}\right\}\right).$$

$$\mu =2k\beta ,$$

$$t_{\text{HAM}}(x,y)=$$

$$-\frac{1}{\mu }{\log }_e{\mathcal{F}}^{-1}\left(\text{LFF}\left\{\frac{\mathcal{F}\left[I_{z_1+z_2}(x,y)/{\tilde{I}}_{z_1+z_2}(x,y)\right]}{1+\pi \gamma z_2\lambda (u^2+v^2)}\right\}\right).$$

$${\Delta }_{\text{eff}}=\Delta -{\Delta }_{\text{emb}}\text{ and }{\beta }_{\text{eff}}=\beta -{\beta }_{\text{emb}}$$

$${\mathcal{L}}_{\sigma }=1+\frac{1}{2}{\sigma }^2{\nabla }_{\perp }^2$$

$$\frac{I_{z_1+z_2}(x,y)}{{\tilde{I}}_{z_1+z_2}(x,y)}-1$$

$$=\left[1-\frac{1}{2}\left(\frac{\gamma z_2}{k}-{\sigma }^2\right){\nabla }_{\perp }^2\right][-2k\beta t(x,y)].$$

$$\gamma z_2⟶\gamma z_2-k{\sigma }^2,\gamma z_2\ge k{\sigma }^2.$$

$$\gamma z_2=k{\sigma }^2,$$

Full text

††preprint: APS/123-QED

I Introduction

X-ray Phase Contrast Imaging (PCI) (see e.g. Paganin Paganin (2006) and references therein) appeared a few decades ago and quickly demonstrated higher diagnostic potential compared to conventional attenuation-based tomography. With the emergence of partially coherent x-ray sources over twenty years ago, expectations regarding PCI became achievable with several methods developed mostly at synchrotrons, some of which were later adapted to laboratory sources Bravin et al. (2013). Among all PCI techniques that utilize additional optical elements, Speckle Based Imaging Berujon et al. (2012); Morgan et al. (2012); Zdora (2018) makes use of a very simple experimental set-up, and has already been demonstrated to be a good candidate for transfer outside synchrotron facilities. The set-up, in addition to its simplicity of implementation, presents the main advantages of having no field of view limitation other than that imposed by the detector (speckle-generating diffusers of any size are easy to manufacture), and relatively low requirements on the beam coherence. As the use of speckle as a wave-front modulator for phase contrast x-ray imaging was proven valuable less than a decade ago Berujon et al. (2012); Morgan et al. (2012), only a few publications have dealt with the corresponding phase-retrieval problem.

To detect the refraction of x rays due to their passage through a sample, the speckle technique has been primarily based on tracking transverse speckle displacement Berujon et al. (2012); Morgan et al. (2012); Zdora (2018) between (i) reference speckle images (i.e. with a membrane generating a random spatial modulation), and (ii) sample speckle images (i.e. with the membrane and the sample). Recently a technique was developed for speckle-tracking phase retrieval based on a geometric-flow formalism Paganin et al. (2018); this approach was inspired by both the transport-of-intensity equation of paraxial wave optics Teague (1983) and the optical-flow formulation of image processing Horn and Schunck (1981); Atcheson et al. (2009). The geometric-flow speckle-tracking method Paganin et al. (2018) has demonstrated good efficiency for retrieving phase shifts of a non-absorbing sample using only a single sample exposure. The method implicitly rather than explicitly tracks speckles, and therefore has potentially an added advantage of finer spatial resolution compared to methods that explicitly track speckles using e.g. correlation windows of at least a few pixels in width and height. The reason for this possibly finer spatial resolution is that, with implicit speckle tracking, the image does not need to be divided into correlation windows. The geometric-flow formalism works by solving a differential equation (continuity equation) that enforces local energy conservation at every pixel Paganin et al. (2018), and so does not need to define multi-pixel regions over which correlation coefficients are computed. This geometric-flow method has also very recently been translated to visible-light microscopy Lu et al. (2019). Here we extend the geometric-flow results of Paganin et al. (2018) to the case of an absorbing monomorphous Paganin et al. (2004); Gureyev et al. (2015) sample, which by assumption obeys at least one of the following approximations: (i) a single-material sample of possibly variable density Paganin et al. (2002); (ii) a sample for which the logarithm of attenuation is proportional to the corresponding phase shift Paganin et al. (2004); Wu et al. (2005), or equivalently a sample of possibly variable density for which the ratio of the real refractive index decrement, and the imaginary part of the complex refractive index, is constant within the volume occupied by the sample Paganin et al. (2004); Gureyev et al. (2015); (iii) a sample composed of light elements ($Z<10$) when irradiated with high-energy x rays (60$-$500 keV) Wu et al. (2005), this being a special case of the single-material assumption Paganin et al. (2002, 2004) in which the single material may be well approximated as electron density. As shall be shown below, an x-ray speckle-tracking experiment, that satisfies any one or more of the above assumptions, will be amenable to analysis using the same reconstruction algorithm.

We close this introduction with a brief overview of the remainder of the paper. Section  II is dedicated to the theoretical development of our method for single-shot x-ray speckle-based imaging of a single-material object. Section III describes the x-ray synchrotron experiment, used to obtain data to implement our method for x-ray phase imaging in both two and three spatial dimensions. The results of applying our theory to these data are given in Sec. IV. We discuss some broader implications of our work, together with possible future extensions and applications, in Sec. V. We conclude with a summary in Sec. VI.

II Theory

Consider unit-intensity monochromatic x-ray scalar plane waves, incident on a spatially random mask, which produce a speckle image downstream—see Fig. 1. Paraxial radiation is assumed throughout, and we use the term speckle to refer to any spatially-random intensity distribution, not necessarily corresponding to fully developed speckle 111As shall become clear later in the paper, the spatially-random intensity distributions that we consider are not fully-developed speckle. Much literature equates speckle with fully developed speckle, but we adopt a more broadly applicable usage of the term speckle to refer to any spatially-random intensity distribution.. If we further assume that the incident plane wave propagates along the positive $z$-axis, which we identify with the optical axis, then such a wave can be represented as $\exp(ikz)$, where $k=2\pi/\lambda$ and $\lambda$ is the wavelength in vacuum. The distance between the mask and the object is $z_{1}$, and the distance between the object and the downstream position-sensitive detector is $z_{2}$. We assume the object to have negligible extension in the $z$-direction.

We assume that the transport-of-intensity equation (TIE) Teague (1983) approximation is valid for this experiment (see Fig. 1). This is because the following validity conditions are, by assumption, met: (i) paraxial x-ray radiation; (ii) thin paraxially-scattering sample; (iii) object-to-detector distance that is sufficiently small for the longitudinal evolution of intensity to be linear with respect to the propagation distance; (iv) negligible inelastic scattering by the sample. Then the longitudinal derivative of the 2D intensity distribution just after the object, $I_{z_{1}}(x,y)$, with respect to the $z$-coordinate can be estimated via Teague (1983):

[TABLE]

where $\varphi(x,y)$ is the phase profile of the wave-field just after the object, $(x,y)$ are Cartesian coordinates in planes perpendicular to $z$ and $\nabla_{\perp}=(\partial/\partial x,\partial/\partial y)$. Here, the phase shift $\varphi(x,y)$ is caused by both the object (as indicated by the “ob” subscript) and the mask (subscript “m”):

[TABLE]

For large Fresnel numbers Paganin (2006), namely for sufficiently small sample-to-detector propagation distances $z_{2}$, Eq. (1) can be simplified as follows:

[TABLE]

Here $I_{z_{1}+z_{2}}(x,y)$ is the 2D intensity distribution after the wave-field propagates the distance $z_{2}$ from the object towards the detector.

Let us rewrite Eq. (II) in a different form:

[TABLE]

In writing the final line of the above equation, we have combined two terms on the right-hand side and neglected higher order terms; note also that we have introduced the notation $(x,y)\equiv{\bf{r}}_{\perp}$. Similar reasoning leads to:

[TABLE]

Here $\tilde{I}_{z_{1}}(x,y)$ is the 2D intensity distribution just before the object (see Fig. 1). As we consider a general case of the object, we may write:

[TABLE]

where ${I}_{\textrm{ob}}({\bf{r}}_{\perp})$ describes the attenuation caused by the object. We can also write an expression similar to Eq. (5) for the case when there is no object in the beam:

[TABLE]

Equations (5) and (7) provide a general formalism describing speckle-based imaging. They also describe the intensity distributions in certain other mask (grating) based phase-contrast imaging techniques. This formalism is somewhat similar to that for combined analyzer-based and propagation-based phase-contrast imaging Pavlov et al. (2004); Coan et al. (2005); Pavlov et al. (2005), where the terms responsible for the shift depend on the gradient of the phase, and the terms in brackets describe free-space propagation effects. Because cross-correlation approaches to phase retrieval (see e.g., Berujon et al. (2012), Morgan et al. (2012), Zdora (2018) and references therein) look simply at sample-induced distortions of the speckle pattern, the propagation-based edge effects described by the terms in brackets are not incorporated. However, ignoring these terms may produce artifacts in the reconstructed images. For instance, the usual procedure to obtain the object attenuation term (uncompensated for the effects of free-space propagation), $\tilde{I}_{\textrm{ob}}({\bf{r}}_{\perp})$, is based on dividing (i) the registered intensity when both the mask and object are in the beam (corrected to mitigate the refraction shift caused by the object), by (ii) the intensity when only the mask is present (see e.g., Morgan et al. (2013) and Zhou et al. (2015)). Thus we may write:

[TABLE]

The uncompensated free-space propagation term in brackets in Eq. (8) can cause fringes in the reconstructed attenuation images to be evident (see e.g. Fig. 2f in Morgan et al. (2013) and Fig. 3b in Zhou et al. (2015)), which require removal via a post-retrieval correction if the attenuation image is to be quantitative Groenendijk et al. (2020). However, if we consider a ratio of unshifted intensities, we obtain an expression similar to the standard TIE (cf. Eq. (II)) having only functions relating to the object in its right-hand side:

[TABLE]

The result shown in Eq. (II) is obtained in the case of a general object. In Eqs. (8) and (II) we have neglected terms containing the averaged scalar product of the gradient of a random field and the gradient of a slowly varying function. However, the solution of this differential equation significantly simplifies in the case of a monomorphous (e.g., single-material) object.

As we have assumed that the object is monomorphous, its complex index of refraction $n$ obeys

[TABLE]

where

[TABLE]

Here, the real numbers $\Delta,\beta$ respectively denote the refractive index decrement and the imaginary (absorptive) part of the complex refractive index. The value of $\gamma$ is considered known from tables or can be experimentally adjusted by trial and error to match the sample composition and density (see e.g., Sanchez del Rio and Dejus (2011)). With reference to Fig. 1, we seek to reconstruct the projected thickness $t(x,y)$ of the thin single-material object, which defines both the object’s phase shift

[TABLE]

and the object’s absorption term

[TABLE]

Then Eq. (II) can be rewritten in the following form:

[TABLE]

In the case of a weakly-attenuating-object approximation, which is relevant to the experimental data used in this paper, we may write $\exp(-2k\beta\,t({\bf{r}}_{\perp}))\approx 1-2k\beta\,t({\bf{r}}_{\perp})$ and hence obtain:

[TABLE]

This is surprising insofar as the data function defined by the left-hand side of Eq. (15) is proportional to the Laplacian $\nabla_{\perp}^{2}$ of the projected thickness, rather than to the components of its transverse first derivative. Finally, the projected thickness in the case of a weakly attenuating object is

[TABLE]

Here, the linear attenuation coefficient $\mu$ is given by

[TABLE]

$\mathcal{F}$ denotes Fourier transformation with respect to $x$ and $y$, $(u,v)$ are Fourier coordinates dual to $(x,y)$, $\mathcal{F}^{-1}$ denotes inverse Fourier transformation with respect to $u$ and $v$, and LFF is a optional high-pass filter. The projected thickness in the case of a relatively highly attenuating monomorphous (HAM) object is:

[TABLE]

If the single-material object is embedded in another material with complex refractive index $n=1-\Delta_{\textrm{emb}}+i\beta_{\textrm{emb}}$, the effective values of

[TABLE]

should be used (cf. Gureyev et al. Gureyev et al. (2002)).

We close this section by considering how the blurring effects of finite source size and detector-based smearing may be taken into account in the reconstruction. Following Subbarao et al. Subbarao et al. (1995), Gureyev et al. Gureyev et al. (2004) and Beltran et al. Beltran et al. (2018), we begin by noting that the linear differential operator

[TABLE]

may be used to approximate the operation of convolving an image with a normalized rotationally symmetric Gaussian point spread function having a standard deviation of $\sigma$. Thus $\mathcal{L}_{\sigma}$ is a diffusion-type linear operator that blurs over a transverse length scale of $\sigma$. It may therefore be used to approximately model the effects of image smearing due to both finite source size and detector-induced blur. Applying this blurring operator $\mathcal{L}_{\sigma}$ to the right side of Eq. (15), expanding the resulting operator product, and then discarding the term containing the bi-Laplacian $\nabla_{\perp}^{2}\nabla_{\perp}^{2}$, gives

[TABLE]

Since this differs from Eq. (15) only through a different coefficient for the Laplacian, the corresponding solution in Eq. (16) still holds, upon making a substitution similar to that given by Beltran et al. in a different context Beltran et al. (2018):

[TABLE]

Interestingly, if the blurring is such that

[TABLE]

the effects of blurring exactly compensate the phase contrast, hence in this particular case the denominator in Eq. (16) may be replaced with unity. This phenomenon is directly analogous to that discussed by Gureyev et al. Gureyev et al. (2004) in the context of propagation-based x-ray phase contrast, whereby the sharpening phase-contrast effects of coherent propagation may be exactly balanced against the blurring effects due to source-size and detector-induced smearing Mittone et al. (2017).

We emphasize that the preceding calculations, which are applied in the present context of our implicit approach to x-ray speckle imaging which does not rely on the use of correlation windows, are more broadly applicable. In particular, the above calculations and their evident Fokker–Planck-type generalizations Risken (1989); Morgan and Paganin (2019); Paganin and Morgan (2019) may also be employed in contexts such as the “X-ray Speckle-Vector Tracking” (XSVT) formalism of Berujon and Ziegler Berujon and Ziegler (2016), together with the “Unified Modulated Pattern Analysis” (UMPA) formalism of Zdora et al. Zdora et al. (2017). Thus, while only one particular aspect of our formalism is implemented in the present paper, we believe it to be applicable to and indeed consistent with the two major formalisms currently employed for x-ray speckle tracking, namely XSVT and UMPA.

III Experiments

Two distinct x-ray experiments, one of which employed monochromated radiation and the other of which employed broad-band polychromatic radiation, were conducted at the ID17 Biomedical Beamline of the European Synchrotron (ESRF).

The first experiment aimed at testing the method under ideal conditions, i.e. with a quasi monochromatic collimated beam. A double Si(111) crystal system in a Laue–Laue configuration was used to select x-ray photons from the synchrotron light source with an energy $E=52$ keV and a monochromaticity defined by the relative energy spread $\Delta E/E\simeq 10^{-4}$. Collimation was defined by the natural divergence of the 21-pole wiggler source which was less than $\sim$1 mrad horizontally and $\sim$0.1 mrad vertically. The photons had first to pass through the speckle-generator membrane before traversing the sample located $z_{1}=900$ mm further away as sketched in Fig. 1. Eventually, the x rays impinged onto a detector placed $z_{2}=12$ m further downstream of the sample. Our detector consisted of an indirect detection system, based on a scientific CMOS camera coupled to magnifying optics imaging a scintillator screen made of a 60 $\mu$m thick gadolinium oxysulfide sheet as active layer. The resulting pixel size was $\simeq$ 6.1 $\mu$m. For this first experiment, a home made phantom composed of nylon wires of 150 $\mu$m in diameter was imaged to assess the accuracy of the method. For these nylon wires at 52 keV we used the following parameters for the complex refractive index: $\beta=0.2809\times 10^{-11};\Delta=0.96842\times 10^{-7}$ Stepanov .

To investigate further the potential applicability of the methodology to more polychromatic beams, another experiment was carried out using a pink beam tomographic configuration. After filtration of the white synchrotron beam by 0.5 mm of Al and 0.35 mm of Cu, the resulting utilized beam exhibited a peak with a mean energy of $E=37.3$ keV and a spectral bandwidth of approximately 20 keV ($\Delta E/E\simeq 0.5$). This second experiment aimed at testing the methodology by imaging in 3D a healthy mouse knee. Such a sample is composed of different tissue types with strongly varying refraction and absorption indices, and may therefore test under real-life conditions the efficacy of the algorithm, even for samples that are not truly weakly absorbing single-material objects. One can recall for instance the method of Paganin et al. Paganin et al. (2002) which eventually proved efficient as well for a much wider range of samples than the class initially set by the strict conditions articulated in that paper. In accordance with Directive 2010/63/EU, the experiments were performed in a mutually agreed animal facility (C3851610006) which was evaluated and authorized by an Ethical Committee for Animal Welfare (APAFIS #13792-201802261434542 v3). Mice were monitored for 10 weeks for any obvious locomotor disability by visual observation and Rotarod testing (Bioseb, Vitrolles, France) before being sacrificed. The imaged mouse knee was fixed with a 4% formaldehyde solution for 48h after its removal and later embedded in a phosphate buffer solution with 2% agarose. For these mouse-knee data the value for $\gamma$ used for the reconstruction was $\gamma=738$, which was obtained from tabulated values of $\Delta=2.98\times 10^{-7}$ and $\beta=4.04\times 10^{-10}$ known for bones Sanchez del Rio and Dejus (2011).

The detector system used to image the mouse knee was equivalent to that which was employed for the first experiment, but this time configured with more strongly magnifying optics which provided a resulting pixel size of $\simeq$ 3 $\mu$m. Overall, the setup layout was also equivalent to the one employed during the first experiment and sketched in Fig. 1, this time with the sample-to-detector and the membrane-to-sample distances respectively set to $z_{2}=4$ m and $z_{1}=1$ m. The tomography data consisted of 3000 projections collected periodically during a 360 degree scan of the sample. The center of rotation was transversely off-centered by 200 pixels to operate a so-called half acquisition tomography scan which is used to extend the 3D field of view. The phase images were calculated for each projection using Eq. (16). The 3D computed-tomography (CT) reconstruction was performed using a filtered back-projection algorithm whose implementation is described in Mirone et al. (2014). For the sake of comparison the same raw data were processed using our earlier algorithm Paganin et al. (2018) and the recovered phase subsequently used for volume reconstruction.

IV Results

Figure 2 presents the results of the calculated thickness on the nylon wire phantom, obtained using the monochromated beam with mean energy 52 keV and energy spread $\Delta E/E\simeq 10^{-4}$. The retrieved thickness is in good agreement with the nominal value of 150 $\mu$m wire diameter, showing the quantitativeness of our approach. In this case the LFF used in Eq. (16) was a Gaussian filter, applied in Fourier space to the zero-padded image, with standard deviation $\sigma^{\prime}=7.5\delta_{f}$. Here $\delta_{f}$ is the step size of the variables u and v in the Fourier space. The $\sigma$ used in the modified form of Eq. (16) as described in Eq. (22) to obtain Fig. 2 is $150~{}\mu m$. The blurring effect, implied by Eq. (22), causes the broadening of the wire’s reconstructed line profile as evident from Fig. 2.

Results of the tomography experiments, obtained using the broad-band polychromatic beam with mean energy $E=37.3$ keV and energy spread $\Delta E/E\simeq 0.5$, are presented in Figs. 3 and 4. Figure 3(a) is a 3D volume rendering of the relative (to the reference material) density of the mouse knee and the subsets are 3$\mu$m thick axial (b), coronal (c) and sagittal (d) plane sections with anatomical annotations provided. The sensitivity of the technique permitted a clear differentiation of the various tissues composing the articulation. Although the hypothesis of a single-material object is not fulfilled here, one can nevertheless easily differentiate the hard tissues such as the bones from the softer tissue, e.g. the ligaments. Moreover, tiny microcalcifications were also made visible with the technique. These microfeatures are usually investigated using conventional preclinical imaging modalities, but they nevertheless remain very difficult to detect. Indeed, the more traditional imaging modalities still encounter strong limitations in the detection of early cartilage and bony changes despite the crucial roles they play in the development of new therapeutic options Bravin et al. (2013). Such a feature of our current method can therefore occupy a place of value within the suite of x-ray phase contrast imaging techniques.

Figure 4 presents for comparison a sagittal thin slice from two 3D reconstructions obtained from the same raw data but using (a) the method of Paganin et al. Paganin et al. (2018) to obtain the function proportional to $\Delta\times k$ of the mouse knee and (b) the herein presented method (without using the LFF) to obtain the relative density of the mouse knee. Insets of certain regions are presented in the right column of the figures.

The yellow square in Fig. 4 marks a joint space region where the boundaries of the different tissues (bone, calcified cartilage, articular cartilage and the intra-articular space) are made visible thanks to the present method while most features remain barely distinguishable with the previous method of Paganin et al. Paganin et al. (2018). The method developed in the present paper allowed us in that case to measure articular cartilage thicknesses, which is of tremendous importance when evaluating for instance osteoarthritis. Such results are indicators of the present method being better able to handle the large difference between the bone ($\Delta=2.98\times 10^{-7}$ and $\beta=4.04\times 10^{-10}$) and cartilage ($\Delta=1.8\times 10^{-7}$ and $\beta=9.03\times 10^{-11}$) refractive index decrements Sanchez del Rio and Dejus (2011). The red insets are located on one meniscus where the number of detectable chondrocytes (cartilage cells depicted as black round structures) appears higher than in the top image. This higher resolution in the hard tissue is confirmed also with the green inset covering a fat-matrix region. Therein, the honeycomb-like matrix is better resolved with the current method than with the former one.

V Discussion

The x-ray speckle-tracking method may be viewed as a form of x-ray Hartmann–Shack sensor Mayo and Sexton (2004); Berto et al. (2017), which employs distortions in random rather than regular arrays of reference intensity maxima, to measure the phase of an x-ray wave front. Moreover, in contrast to typical approaches to both Hartmann–Shack sensing and speckle tracking—in an x-ray setting and more broadly—our method implicitly rather than explicitly tracks speckles Paganin et al. (2018). As noted earlier, this potentially increases the intrinsic resolution of the reconstruction since correlation windows, which must be at least a few pixels by a few pixels in size, are not needed. Further, our method is computationally simple and rapid to implement when compared to approaches based on optimizing error metrics in high-dimensional parameter spaces: regarding the latter, see references contained in the recent review by Zdora Zdora (2018). We also note some parallels with earlier work in the visible-light regime, by Massig Massig (1999, 2001) and Perciante et al. Perciante and Ferrari (2000).

Interestingly, the Fourier filtration in Eq. (18) is mathematically identical in form to the method of Paganin et al. Paganin et al. (2002) (PM) for propagation-based x-ray phase contrast imaging of a single-material object. Two conclusions immediately follow. (i) Observe an interesting mathematical similarity, in the sense remarked upon earlier in this paragraph, between propagation-based phase contrast methods (which use no mask and do not rely on speckle) and speckle-tracking phase contrast methods (which use a mask to generate illuminating speckles). These two classes of method have typically been treated as distinct hitherto. (ii) The previously noted similarity opens the possibility for working in very low-dose environments, on account of the PM’s very high robustness with respect to noise. For example, for the PM, 2D imaging with only one imaging quantum per pixel is possible Clark et al. (2019). Several papers Beltran et al. (2010, 2011); Nesterets and Gureyev (2014); Gureyev et al. (2014); Kitchen et al. (2017); Gureyev et al. (2017) show that signal-to-noise ratio (SNR) boosts well in excess of 100 are possible with the PM, which may be traded against dose (acquisition time) reductions of four orders of magnitude or more. This is of particular importance for contexts such as (a) dynamic samples and real-time x-ray imaging; (b) radiation-sensitive samples; (c) low-flux sources such as laboratory x-ray sources; (d) high-throughput scenarios such as industrial-product quality-control inspection.

While the reconstructions in the present paper were of sufficient quality to not necessitate further refinement, it is also important to note in very-low-flux contexts that our algorithm could be used as an initial estimate to seed iterative refinement engines, thus enabling improved reconstructions that explicitly leverage relevant a priori knowledge such as real-space sparsity, Fourier-space sparsity etc. This is aligned with the usual trade-off that can be made where appropriate, in which an increase in suitable a priori knowledge may be traded off against the amount of data that one obtains (and therefore the dose incurred): more a priori knowledge allows one to make use of less data at the price of reduced generality, whereas less a priori knowledge requires more data to achieve reconstructions that are valid under more-general conditions. Recent achievements in phase-retrieval utilizing both compressive sensing and neural networks Kemp (2018); Rivenson et al. (2018); Jiang et al. (2018) are particularly relevant in this regard, albeit beyond the scope of the present paper.

It is also interesting to note that our method overcomes a key issue for TIE imaging for x rays without the assumption of a single-material object, namely the need to acquire images in parallel planes for different $z$ values along the optical axis in order to retrieve quantitative sample information. While methods exist for single-plane x-ray TIE imaging, based e.g. on measurements taken at different energies Li et al. (2018); Gürsoy and Das (2013), the method developed in the present paper is considerably more stable with respect to noise, again on account of its similarity in mathematical form to the PM.

We close this discussion by noting that source code used for all analyses in this paper is freely available at:

https://github.com/labrieth/spytlab.

VI Conclusion

We have developed and tested a means for x-ray speckle-tracking phase imaging, which extends our previous work approaching the problem using the concept of geometric flow Paganin et al. (2018). The method decodes the deformations in a single speckle field, created by a mask and subsequently deformed by an object, by applying transport-of-intensity concepts to this speckle-field scenario and then solving the inverse problem associated with the resulting equations. The ability to work with a single exposure, once the initial reference-speckle image has been taken in the absence of a sample, allows the method to be applied to time-dependent data. The stability with respect to noise of the method, which may be traded off against reduced dose and therefore reduced data acquisition time, means that our single-image method may be applicable to time-dependent imaging of dynamic samples. This ability to work with dynamic samples might be improved still further by the fact that coherence is not used explicitly in our method, which essentially relies on geometric optics alone, with associated very relaxed temporal coherence requirements. This enables it to be used with broad-band polychromatic beams, whose higher flux relative to monochromated beams may enable frame rates to be increased still further, in the context of imaging time-dependent samples. More generally, we believe that our method for implicitly tracking x-ray speckles has the capacity to make x-ray speckle tracking into a significantly more practical, dose-efficient, stable, rapid technique. Lastly, given the very recent translation of our previously published method for implicit speckle tracking Paganin et al. (2018) to visible-light microscopy Lu et al. (2019), we anticipate that the method developed in the present paper may also be of utility in the visible-light domain.

Acknowledgements.

We acknowledge useful discussions with Claudio Ferrero, Timur Gureyev, Andrew Kingston, Marcus Kitchen, Sheridan Mayo, Kaye Morgan, and Gary Ruben. Financial support from the Experiment Division of the European Synchrotron Radiation Facility (ESRF) for DMP to visit in early 2018 is gratefully acknowledged. EB acknowledges support from LabEx PRIMES (ANR-11-LABX-0063/ANR-11-IDEX-0007).