Low-exposure, high-quality multimodal speckle X-ray imaging via an intrinsic gradient-flow approach

TL;DR

This paper introduces a gradient-flow variation of the MIST algorithm for speckle-based X-ray imaging, reducing data needs and improving dark-field image quality, especially for low-density samples, with potential for low-exposure applications.

Contribution

The paper presents a novel gradient-flow MIST method that enhances speckle X-ray imaging by requiring less data and producing higher quality dark-field images for weakly attenuating samples.

Findings

Reduces data requirements for image retrieval.

Improves dark-field image quality for low-density samples.

Demonstrated on experimental data of carbon fibers at synchrotron.

Abstract

We present a new approach for retrieving dark-field, phase shift, and attenuation images from speckle-based X-ray imaging data. Speckle-based X-ray imaging (SBXI) exploits sample-induced alterations to a reference near-field speckle pattern produced by a randomly structured mask. Attenuation images allow materials of different densities to be visualised. Phase-shift images are useful because they reveal how materials in a sample refract the X-ray beam, providing contrast between similar low-density structures that are difficult to reconstruct in attenuation images. Dark-field images convey information about structures that are smaller than the spatial resolution and thus invisible in both attenuation and phase-shift images. In previous works, we presented the Multimodal Intrinsic Speckle-Tracking (MIST) algorithm, which recovers the three complementary imaging modes from SBXI data by solving the associated Fokker--Planck equation. In this work, we present a variation of MIST, called ``gradient-flow MIST", which (1) reduces the amount of SBXI data required for image retrieval, (2) maintains the full generality of the X-ray Fokker--Planck equation, and (3) recovers dark-field images with higher quality than the previously proposed variants for weakly attenuating (i.e., low density) samples. We demonstrate the new gradient-flow MIST approach on experimental SBXI data of a knotted bundle of carbon fibres acquired at the Australian synchrotron. This approach is anticipated to be useful in phase-contrast and dark-field applications that require simplicity in experimentation and low sample X-ray exposure.

Equations

$$I_R({\text{r}}_{\perp },Z)I_{\text{ob}}({\text{r}}_{\perp },0)-I_S({\text{r}}_{\perp },Z)=$$

$${\nabla }_{\perp }\cdot \left[\frac{Z}{k}{I}_R({\text{r}}_{\perp },Z)I_{\text{ob}}({\text{r}}_{\perp },0){\nabla }_{\perp }\varphi ({\text{r}}_{\perp },0)\right]-$$

$$Z^2{\nabla }_{\perp }^2\left[D({\text{r}}_{\perp },0)I_R({\text{r}}_{\perp },Z)I_{\text{ob}}({\text{r}}_{\perp },0)\right].$$

$$\frac{I_R{I}_{\text{ob}}-I_S}{Z}-{\nabla }_{\perp }\cdot \left[\frac{1}{k}{I}_R{I}_{\text{ob}}{\nabla }_{\perp }\varphi -Z{\nabla }_{\perp }\left(D{I}_R{I}_{\text{ob}}\right)\right]=0,$$

$${\text{J}}_{\perp }=\frac{1}{k}{I}_R{I}_{\text{ob}}{\nabla }_{\perp }\varphi -Z{\nabla }_{\perp }\left(D{I}_R{I}_{\text{ob}}\right).$$

$$\frac{1}{k}{I}_R{I}_{\text{ob}}{\nabla }_{\perp }\varphi -Z{\nabla }_{\perp }\left(D{I}_{\text{ob}}{I}_R\right)={\nabla }_{\perp }V+\text{rot}(A).$$

$$\frac{I_R{I}_{\text{ob}}-I_S}{Z}={\nabla }_{\perp }^{2}V.$$

$$I_{\text{ob}}^{(1)}\approx {\mathcal{F}}^{-1}\frac{1}{1+\frac{\gamma Z}{2k}{\text{k}}_{\perp }^2}\mathcal{F}\left(\frac{I_S}{I_R}\right),$$

$$I_R{\alpha }_1+Z({\partial }_x{I}_R){\alpha }_3\approx {\partial }_{x}V,$$

$$I_R{\alpha }_2+Z({\partial }_y{I}_R){\alpha }_3\approx {\partial }_{y}V,$$

$${\alpha }_1$$

$${\alpha }_2$$

$${\alpha }_3$$

$${\nabla }_{\perp }V={\mathcal{F}}^{-1}{f}_d{f}_L\mathcal{F}\left[\frac{I_R{I}_{\text{ob}}-I_S}{Z}\right]$$

$$f_d=\left(\frac{i\text{sin}(k_x\mathcal{W})}{\mathcal{W}},\frac{i\text{sin}(k_y\mathcal{W})}{\mathcal{W}}\right),$$

$$f_L=\frac{{\mathcal{W}}^2}{2\text{cos}(k_x\mathcal{W})+2\text{cos}(k_y\mathcal{W})-4}.$$

$$\left(1-\frac{\gamma Z}{2k}{\nabla }_{\perp }^2\right)I_{\text{ob}}=g_1+g_2$$

$$g_1({\text{r}}_{\perp })\equiv \frac{I_S}{I_R},$$

$$g_2({\text{r}}_{\perp })\equiv \frac{Z}{I_R}({\alpha }_1{\partial }_x{I}_R+{\alpha }_2{\partial }_y{I}_R$$

$$-Z[I_R{\nabla }_{\perp }^2{\alpha }_3+{\alpha }_3{\nabla }_{\perp }^2{I}_R+{\nabla }_{\perp }{\alpha }_3\cdot {\nabla }_{\perp }{I}_R]).$$

$$I_{\text{ob}}^{(2)}={\mathcal{F}}^{-1}\frac{1}{1+\frac{\gamma Z}{2k}{\text{k}}_{\perp }^2}\mathcal{F}\left(g_1+g_2\right).$$

$$D=-\frac{{\alpha }_3}{I_{\text{ob}}^{(2)}}.$$

$$\varphi =\frac{\gamma }{2}\text{ln}\left[I_{\text{ob}}^{(2)}\right].$$

$$D=-\frac{{\alpha }_3}{I_{\text{ob}}},$$

$$I_{\text{ob}}-\frac{Z}{k}\left[\nabla \cdot I_{\text{ob}}\nabla \varphi \right]=g_1+g_2,$$

$$\text{FRC}=\frac{{\sum }_{|{\text{k}}_{\perp }|}{F}_1({\text{k}}_{\perp })\cdot F_2^{\ast }({\text{k}}_{\perp })}{\sqrt{{\sum }_{|{\text{k}}_{\perp }|}|F_1({\text{k}}_{\perp }){|}^2\cdot {\sum }_{|{\text{k}}_{\perp }|}|F_2({\text{k}}_{\perp }){|}^2}}$$

$$\text{CNR}=\frac{|{\mu }_{\text{s}}-{\mu }_{\text{b}}|}{{\sigma }_{\text{s}}}$$

Full text

††preprint: APS/123-QED

I Introduction

Conventional attenuation-based X-ray imaging (like clinical radiography) provides low image contrast between similarly attenuating structures and yields little information on structures that are below the imaging system’s spatial resolution. Phase-contrast imaging (PCI) [10.1063/1.1146073] and dark-field imaging (DFI) [1992RScI...63..611M] are two techniques that address these respective limitations [Pagnin_2006, ENDRIZZI201888, Levine_Long_2004, Moosmann2013, Moosmann_Ershov_Weinhardt_Baumbach_Prasad_LaBonne_Xiao_Kashef_Hofmann_2014]. The basis of attenuation-based imaging has remained unchanged since its discovery in 1896 [Rontgen1895], relying solely on differences in X-ray attenuation along the X-ray’s trajectory to generate contrast. PCI, on the other hand, additionally exploits differences in phase shifts induced by different structures in the sample to generate contrast. The contrast in DFI arises from multiple contrast mechanisms, such as, small angle X-ray scattering from microstructures and Young–Maggi–Rubinowicz boundary-diffraction from edges [Glatter_Kratky_1982, Levine_Long_2004, Miyamoto:62]. These additional channels of information allow PCI and DFI to reveal features that are invisible in attenuation images [Pfeiffer_Weitkamp_Bunk_David_2006, Pfeiffer_Bech_Bunk_Kraft_Eikenberry_Brönnimann_Grünzweig_David_2008]. For example, in medical imaging, phase-retrieved phase-contrast images of the breast show better image quality than their respective attenuation images [Gunaseelan2023, f0dbd53b8e70479faffed8ff590edda7, TavakoliTaba2020, Gureyev:tv5076], and dark-field images allow pathophysiological changes in the lung to be more easily visualised [velroyen2015grating].

Currently, a variety of PCI and DFI techniques have been demonstrated, such as analyser-based [Förster_Goetz_Zaumseil_1980, Chapman_Thomlinson_Johnston_Washburn_Pisano_Gmür_Zhong_Menk_Arfelli_Sayers_1997], propagation-based [Paganin_Mayo_Gureyev_Miller_Wilkins_2002, Gureyev_2020], grating-based [Momose_Kawamoto_Koyama_Hamaishi_Takai_Suzuki_2003, Pfeiffer_Bech_Bunk_Kraft_Eikenberry_Brönnimann_Grünzweig_David_2008], edge illumination [Olivo_Arfelli_Cantatore_Longo_Menk_Pani_Prest_Poropat_Rigon_Tromba_2001, PhysRevApplied.19.054042], and speckle-based X-ray imaging (SBXI) [Morgan_Paganin_Siu_2012, Peverini_2012]. All these methods can combine attenuation imaging, PCI, and DFI; hence, they are multimodal. These multimodal techniques have also been shown to work with non-synchrotron X-ray sources, making them strong candidates for mainstream use, for example, in clinics [Russo_2017].

Since its realisation, SBXI has been rising in popularity because it is easy to implement experimentally. SBXI uses a speckle pattern in the near-field regime as an X-ray wavefront marker. An X-ray speckle pattern is generated by inserting a structured, usually random, mask with spatially absorbing features and/or phase-shifting features, such as sandpaper, in a partially coherent X-ray beam [Dainty_1975]. When a sample is placed in the beam downstream of the mask, the reference speckle pattern will be modulated according to the structural features in the sample. By mathematically tracking how speckles in the reference speckle pattern $I_{R}$ have transformed into the sample speckle pattern $I_{S}$, the structural and compositional features in the sample can be extracted [PhysRevLett.112.253903]. Different modulations to the reference speckle pattern are attributed to different sample-induced contrast mechanisms: a decrease in intensity reflects that the X-ray beam has been attenuated in that region of the sample, transverse speckle-shifts are due to a phase change of the beam in the sample, and blurring of the speckle pattern–which is also called a local decrease in visibility–is due to dark-field effects [Berujon_Wang_Sawhney_2012, Zdora_Thibault_Zhou_Koch_Romell_Sala_Last_Rau_Zanette_2017, cerbino2008x].

Initially, speckle modulations were tracked by mathematically comparing $I_{R}$ and $I_{S}$ in small regions [Zdora_2018, Celestre2025]. Examples of such methods include X-ray speckle tracking (XST) [Morgan_Paganin_Siu_2012, Peverini_2012], X-ray speckle scanning (XSS) [Peverini_2012] or unified modulated pattern analysis (UMPA) [Zdora_Thibault_Zhou_Koch_Romell_Sala_Last_Rau_Zanette_2017]. In 2018, another paradigm of speckle tracking was realised and is called the optical–flow method [Paganin_Labriet_Brun_Berujon_2018], which employed an equation similar in spirit to the widely adopted transport-of-intensity equation (TIE) [Teague_1983] to model the evolution of the reference speckles in $I_{R}$ into the sample-modified speckles in $I_{S}$. In comparison to the aforementioned XST, XSS, or UMPA, the speckles are tracked globally, instead of locally, through the entire image, by solving this TIE–like equation. The initial optical–flow method presented in Paganin_Labriet_Brun_Berujon_2018 was limited to the case of phase-objects, i.e., attenuation and dark-field effects were neglected. An extension, named ‘Multimodal Intrinsic Speckle Tracking’ (MIST), was realised in 2020 in Pavlov_Paganin_Li_Berujon_Rougé-Labriet_Brun_2020 which utilises solving the more general X-ray-Fokker–Planck equation [Paganin_Morgan_2019]. The X-ray Fokker–Planck equation extends the TIE by also including dark-field effects. The inclusion of dark-field effects allows the reconstruction of an image that reveals information about unresolved microstructures, while also improving the quality of the phase-shift image [Alloo_2025]. MIST therefore reconstructs three complementary images: attenuation, phase-contrast, and dark-field [Alloo_Paganin_Morgan_Kitchen_Stevenson_Mayo_Li_Kennedy_Maksimenko_Bowden2022].

Since its introduction in 2020 by Pavlov_Paganin_Li_Berujon_Rougé-Labriet_Brun_2020, several variants of the MIST algorithm have been proposed [Alloo_Paganin_Morgan_Kitchen_Stevenson_Mayo_Li_Kennedy_Maksimenko_Bowden2022, Pavlov_Paganin_Morgan_Li_Berujon_Quénot_Brun_2021, Alloo_Morgan_Paganin_Pavlov_2023] to improve the stability of the solutions and increase the generality of the retrieval approach. Each variant differs in its approach to solving the X-ray Fokker–Planck equation and requires at least two or more speckle image pairs. A speckle image pair consists of the images $I_{R}$ and $I_{S}$; different pairs are obtained by shifting the speckle-generating mask and collecting $I_{R}$ and $I_{S}$ images at each mask position. The first variant of MIST [Pavlov_Paganin_Li_Berujon_Rougé-Labriet_Brun_2020] assumes the sample is a) pure-phase and b) the dark-field is spatially slow-varying, to simplify solving the X-ray Fokker–Planck equation. Assumption a) was relaxed in Alloo_Paganin_Morgan_Kitchen_Stevenson_Mayo_Li_Kennedy_Maksimenko_Bowden2022, as well as extending MIST to tomography. Assumption b) may lead to artefacts at sharp interfaces in the reconstructed dark-field images and was later relaxed in Alloo_Morgan_Paganin_Pavlov_2023. The most recent, superior in terms of image quality, and most general MIST variant for isotropic dark-field [Alloo_Morgan_Paganin_Pavlov_2023] requires a minimum of four speckle image pairs, which is an increase from the original variant by a factor of two [Pavlov_Paganin_Li_Berujon_Rougé-Labriet_Brun_2020]. This may limit its applicability in scenarios where image acquisition time should be minimised, such as medical imaging.

In this paper, we present a new MIST variant called gradient-flow MIST (GF-MIST). GF-MIST requires only two speckle image pairs, is robust to noise, and maintains the same generality as [Alloo_Morgan_Paganin_Pavlov_2023], making it applicable to a broad range of samples. We derive GF-MIST by solving the Fokker–Planck equation for SBXI, and then demonstrate the method’s effectiveness on SBXI data of a sample of a knotted carbon fibre collected at the Australian synchrotron.

II Theory

Consider the SBXI setup for monochromatic paraxial illumination as shown in Figure 1. The evolution of the reference-speckle pattern $I_{R}$ into the sample-speckle pattern $I_{S}$ can be modelled by the finite-difference form of the X-ray Fokker–Planck equation [Alloo_Paganin_Morgan_Kitchen_Stevenson_Mayo_Li_Kennedy_Maksimenko_Bowden2022]:

[TABLE]

Above, $\textbf{r}_{\perp}\equiv(x,y)$ denotes the Cartesian coordinates perpendicular to the optical $z$ axis, $\nabla_{\perp}\equiv(\partial_{x},\partial_{y})$ is the 2-D transverse gradient operator, $k$ is the wavenumber given by $2\pi/\lambda$, where $\lambda$ is the wavelength, $I_{\text{ob}}(\textbf{r}_{\perp},0)$ is the sample’s transmission at the sample’s exit-plane $z=0$, $I_{R}(\textbf{r}_{\perp},Z)$ and $I_{S}(\textbf{r}_{\perp},Z)$ are the reference and sample speckle images, respectively, collected at a distance $z=Z$ from the sample, $\phi(\textbf{r}_{\perp},0)$ is the sample-induced phase-shift, and $D(\textbf{r}_{\perp},0)$ is the dimensionless effective diffusion coefficient which characterises local isotropic dark-field effects. Equation II is the SBXI Fokker–Planck equation, accounting for contributions from attenuation, phase-shift, and dark-field effects. This equation can be rearranged into a form that is more suggestive of a continuity equation,

[TABLE]

where the functional dependencies are omitted for clear presentation. The conserved quantity in this continuity equation  (2) is the integral of the quantity $I_{R}I_{\text{ob}}-I_{S}$ over a 2-D region (energy), and the conserved current $\textbf{J}_{\perp}(\textbf{r}_{\perp},0)$ is given by the term inside the divergence operator:

[TABLE]

The current $\textbf{J}_{\perp}$ describes the flow of energy in the transverse direction as $I_{R}I_{\text{ob}}$ is transformed into $I_{S}$. This current term comprises two parts that govern the energy flow: a coherent and a diffusive transport part [Paganin_Morgan_2019]. These two flow mechanisms together give the forward model describing the evolution of $I_{R}$ into $I_{S}$, and by solving the associated inverse problem, the attenuation, phase-shift, and dark-field can be recovered.

To begin solving the inverse problem, that is retrieving $\phi$ and $D$, we take the 2-D Helmholtz decomposition of $\textbf{J}_{\perp}$ into an irrotational and rotational part, characterised by $V(\textbf{r}_{\perp},0)$ and $A(\textbf{r}_{\perp},0)$, respectively (cf. Teague_1983), namely [Schmalz_Gureyev_Paganin_Pavlov_2011]

[TABLE]

Above, the rotation operator rot maps $A$ into the vector $(\partial_{y}A,-\partial_{x}A)$ [Schmalz_Gureyev_Paganin_Pavlov_2011]. In Eq. (4), we neglect the rotational component $\textbf{rot}(A)$, which is a valid assumption if the sample’s attenuation does not vary rapidly [Schmalz_Gureyev_Paganin_Pavlov_2011]. By applying the 2-D divergence to both sides of Eq. (4) and comparing with Eq. (2), the function $V$, which is related to the phase-shift and dark-field images we aim to recover, can be obtained by numerically solving the following Poisson equation:

[TABLE]

Equation (5) must be solved subject to appropriate boundary conditions. If the sample lies within the entire field of view of the detector, one can impose fixed Dirichlet boundary conditions [Paganin_Pelliccia_Morgan_2023] as there is no phase-shift or dark-field outside the sample, and hence, at the boundary of the images. The transmission image $I_{\text{ob}}$ required to solve Eq. (5) can be obtained by acquiring a contact image, i.e., an attenuation-only image, during imaging (i.e., $Z=0$). If the sample is approximately single-material, a first order approximation for $I_{\text{ob}}$ from a phase-retrieval solution derived from the TIE for SBXI can be used [Pavlov_Li_Paganin_Berujon_Rougé-Labriet_Brun_2020]:

[TABLE]

where $\mathcal{F}$ is the 2-D Fourier transform, $\textbf{k}_{\perp}$ is the 2-D spatial frequency vector, and $\gamma$ (known a priori) is the ratio of the real part of the sample’s refractive index to its imaginary, $\gamma=\delta/\beta$. Note that the refractive index, $n$, is given by $n=1-\delta+i\beta$, where $\delta$ and $\beta$ characterise the sample’s refraction and attenuation properties, respectively [Pagnin_2006]. We label $I_{\text{ob}}^{(1)}$ with a superscript ‘1’ to indicate that it is a first-order approximation derived using the TIE and does not account for dark-field effects. An iterative procedure that reintroduces these image-blur-inducing dark-field effects into the solution for $I_{\text{ob}}$ is described later, thereby improving the transmission retrieval.

Once $V$ is retrieved, it can be used to establish the following set of linear equations,

[TABLE]

in which the unknown variables are linked to the phase-shift and dark-field images to be retrieved via $\alpha_{1}$–$\alpha_{3}$, which are defined as

[TABLE]

The variables $\alpha_{1}$ and $\alpha_{2}$ are interpreted as being associated with the refraction angles in the $x$ and $y$ directions, respectively, which consider dark-field effects. The third variable, $\alpha_{3}$, is a thickness-dependent dark-field term. Numerically, $\nabla_{\perp}V$ can be computed on a square array of pixels with pixel size $\mathcal{W}\times\mathcal{W}$ using Fourier filtering,

[TABLE]

where $f_{d}$ is the noise-robust derivative filter [Fornberg_1987],

[TABLE]

and $f_{L}$ is the inverse $\nabla_{\perp}^{2}$ filter [Paganin_Favre-Nicolin_Mirone_Rack_Vilanova_Olbinado_Fernandez_daSilva_Pelliccia_2020, Press2007],

[TABLE]

The filter $f_{L}$ is numerically unstable at the origin due to division by zero. However, as we only require Laplacian inversion to recover the gradient of $V$, we can simply set the solution of the inversion to zero at the origin. This corresponds to shifting $V$ by a constant value, which is possible since the gradient of $V$ gives the same solution as an un-shifted variant. Equations (7)–(8) are equivalent to Eq. (4) but with the orthogonal $x$ and $y$ components forming one equation each, and the rotational part (i.e, $\textbf{rot}(A)$) neglected. This linear system is underdetermined, i.e., it contains more unknowns than equations, and therefore no unique solution exists. An approach applied in all MIST algorithms is to construct additional sets of equations by using multiple sets of speckle image pairs. A solution that minimises the residuals of this overdetermined system can be computed using a linear least squares approach [Alloo_Morgan_Paganin_Pavlov_2023], such as Tikhonov regularised-QR decomposition [tikhonov1977solutions]. However, if only one set of equations is available, then the system is underdetermined.

When Eq. (6) is used to approximate $I_{\text{ob}}$, an extra step may be required to correct for dark-field effects if the sample is strongly scattering, as the TIE-retrieved transmission image will be significantly blurred in these areas. To consider dark-field effects in the attenuation retrieval step, Eq. (6) should be extended to include the diffusion term in the X-ray Fokker–Planck equation. This is achieved by expressing a closed-form solution to $I_{\text{ob}}$ from Eq. (2) instead of from the TIE. Making use of the projection approximation, $\phi=-k\delta T$, where $T$ is the projected thickness of the sample, and Beer’s law of absorption, $I_{\text{ob}}=\text{exp}(-2k\beta T)$ [Pagnin_2006, Paganin_Nugent_1998, Pavlov_Li_Paganin_Berujon_Rougé-Labriet_Brun_2020], we rewrite Eq. (2) as

[TABLE]

where we define the source terms

[TABLE]

[TABLE]

Equation (15) has a closed-form solution for the transmission which considers dark-field effects:

[TABLE]

Here, we have separated the dominant term, $g_{1}$, from $g_{2}$, which represents small second-order correction terms to $g_{1}$. The unknown variables $\alpha_{1}$, $\alpha_{2}$, and $\alpha_{3}$ used in the right-hand side of Eq. (18) are computed with $I_{\text{ob}}^{(1)}$ first. The diffusion coefficient dark-field can then be retrieved directly using the solution $\alpha_{3}$ and the updated transmission $I_{\text{ob}}^{(2)}$, that is,

[TABLE]

The phase-shift can then be retrieved using

[TABLE]

In the case when a contact image is used to obtain an exact value of $I_{\text{ob}}$, then dark-field is retrieved from

[TABLE]

and the phase-shift is found by numerically solving the following differential equation:

[TABLE]

for example, using finite difference methods [Press2007]. Equation (22) was formed by expressing Eq (2) in terms of $g_{1}$ and $g_{2}$.

III Methods

To validate our approach, SBXI data of a knotted carbon fibre was collected at the Australian synchrotron’s Imaging and Medical beamline (IMBL) using a setup similar to that shown in Fig. 1. The knotted carbon fibre was placed $Z=3$ m from the detector and $Z_{MS}=1$ m from the speckle-generating mask. The mask consisted of a stack of six sandpapers with grit size P240. X-rays of 25 keV with an energy bandwidth of $\Delta E/E=10^{-3}$ were used to illuminate the mask and the sample. Twenty sets of speckle image pairs were collected by transversely shifting the sandpaper. The effective pixel size for the detector was 10.2 $\mu$m, while the speckle pattern had an effective speckle size of 105.9 $\mu$m, which was computed by fitting the autocorrelation of $I_{R}$ to a Gaussian function and taking the full-width-half-maximum [goodman2007speckle, 88d937a3476a42f6b9b5490418b4da06].

III.1 Numerical Implementation

Our developed GF-MIST algorithm was implemented in MATLAB (version 2023b), with the script and test data openly available at [GF-MIST]. For this carbon fibre sample, we assumed it to be single-material under 25 keV X-ray illumination and therefore had a $\gamma$ value of $2600$, which was computed from the TS imaging calculator [AAA]. The transmission was then initially approximated using the TIE-derived solution (Eq. (6)). The overdetermined system of linear equations generated by twenty instances of Eqs. (7)-(8) (totalling 40) was solved with Tikhonov-regularised QR decomposition using a Tikhonov parameter of $4.6\times 10^{-5}$, determined using the L-curve method [Hansen_O’Leary_1993]. To then consider dark-field effects, thereby correcting for scatter-induced image blur, the transmission was updated using Eq. (18). The dark-field and phase-shift were then extracted using Eqs. (19) and (20) using the updated transmission.

III.2 Evaluation of Image Quality

We evaluated GF-MIST’s performance by comparing its retrieved images to those obtained using the MIST approach developed in [Alloo_Morgan_Paganin_Pavlov_2023], which we call ‘RV-MIST’ hereafter. We used the openly available Python script at [Alloo2023MIST] to obtain the RV-MIST reconstructions.

The Fourier ring correlation (FRC) was used to quantify the spatial resolution of the images, that is, the smallest length scale resolvable in the images. The FRC measures the correlation of two images in the frequency domain and is defined as

[TABLE]

where $F_{1}(\textbf{k}_{\perp})$ and $F_{2}(\textbf{k}_{\perp})$ are the Fourier transforms of image 1 and image 2, respectively, of the same sample. Images 1 and 2 were obtained by splitting the full twenty SBXI dataset in half and applying the GF-MIST/RV-MIST to each half. As the FRC curves were noisy, they were fitted to a complementary error function (erfc). The 0.143 FRC threshold criterion [VANHEEL2005250] was then used as a cutoff for the spatial resolution.

The contrast-to-noise ratio (CNR) was used to measure the visibility of structures within the images, which was computed using

[TABLE]

where $\mu_{\text{s}}$ and $\sigma_{\text{s}}$ are the mean and standard deviation of the intensity of an area of a flat signal, respectively, and $\mu_{\text{b}}$ is the mean background intensity.

IV Results and Discussion

Figure 2 shows the GF-MIST-retrieved dark-field, refraction angles in the transverse directions, phase-shift, and transmission images reconstructed using the entirety of the SBXI of the carbon knot, i.e., 20 speckle image pairs. Notably, GF-MIST is the first Fokker–Planck-equation-derived algorithm for SBXI to retain the speckle-shift term $\nabla I_{R}\cdot\nabla\phi$ [Pavlov_Li_Paganin_Berujon_Rougé-Labriet_Brun_2020], and hence can reconstruct the two orthogonal refraction angles (i.e., consider phase gradients within the sample) without requiring to calculate the phase-shift first. Due to their main orientation, the fibres are the most visible in the reconstruction of the $x$ refraction angle. However, they can also be seen on the $y$ refraction angle, phase-shift, transmission, and partially on the dark field images.

Figure 4 compares dark-field and phase-shift images retrieved using GF-MIST and RV-MIST, and shows how the retrieved images are affected by the number of input speckle image pairs. RV-MIST requires at least four speckle image pairs, whereas GF-MIST needs a minimum of two. The magnitude of the GF-MIST’s retrieved dark-field signal is $\sim 2$ times that of that obtained using RV-MIST, and hence, the maximum and minimum displayed greyscale values in Fig. 4 had to be set differently for each to allow for appropriate qualitative comparison. The difference in the retrieved magnitude of the dark-field signal may arise from the fact that GF-MIST computes first derivatives of $I_{R}$ when constructing its system of linear equations, whereas RV-MIST computes the second derivative (i.e., transverse Laplacian operator) of $I_{R}$. This reduction in one derivative order results in less high-frequency noise amplification, allowing more dark-field signal to be retrieved through the noise. In other words, from the continuous Fourier derivative theorem, the first derivative of $I_{R}$ is equivalent to linear scaling by $i(k_{x},k_{y})$ in the frequency domain, whereas the second derivative is quadratically scaled by $\textbf{k}_{\perp}^{2}$, which is less numerically robust to noise. For both methods, however, the phase-shifts quantitatively agree. There also appear to be vertical stripe artefacts that are reconstructed in the dark-field images of RV-MIST, but not in GF-MIST.

To compare the image quality of the dark-field and phase-shift images retrieved using RV-MIST and GF-MIST across different numbers of speckle image pairs, the spatial

resolution and CNR were computed. Figure 3 shows FRC curves for both approaches. The general trend for both approaches is that increasing the number of speckle image pairs shifts the FRC curve to the right, indicating improved spatial resolution.

Figure 5a plots the FRC-computed spatial resolution, as a function of speckle image pairs used to reconstruct the dark-field and phase-shift images using GF-MIST and RV-MIST. Figure 5b does the same but with the computed CNR. The dark field spatial resolution of RV-MIST approximately lags behind GF-MIST by two speckle image pairs. For example, the dark field image retrieved from GF-MIST with four pairs yields comparable resolution to that obtained from RV-MIST using six. This is particularly noteworthy, as it suggests that comparable image quality can be achieved with reduced sample X-ray exposure.

Regarding the phase shift, the resolution of GF-MIST is approximately one pixel lower than that of RV-MIST. This is because computing terms such as $\alpha_{3}\nabla_{\perp}^{2}I_{R}$ and $\nabla_{\perp}\alpha_{3}\cdot\nabla_{\perp}I_{R}$ in the correction factor $g_{2}$ (Eq. (18)), for reasons mentioned previously, amplifies noise. For the same reasons, the CNR for the GF-MIST-retrieved dark-field images is up to $\sim 10$ times greater than RV-MIST using the same number of speckle imaging pairs, and slightly below RV-MIST for the phase-shift images.

While we demonstrate GF-MIST’s superior dark-field performance over RV-MIST in this manuscript, this result corresponds to the ideal case for GF-MIST–namely, a low-attenuating sample. Evidently, there will be cases in which RV-MIST outperforms GF-MIST, and even cases where earlier versions of MIST, such as those proposed by Pavlov_Paganin_Li_Berujon_Rougé-Labriet_Brun_2020 and Alloo_Paganin_Morgan_Kitchen_Stevenson_Mayo_Li_Kennedy_Maksimenko_Bowden2022, may outperform both RV- and GF-MIST. Although earlier variants of MIST rely on stronger assumptions, they are mathematically simpler to solve, which can affect the stability and, consequently, the quality of the retrieved images. It would be valuable for future work to investigate the performance of RV-MIST, GF-MIST, and earlier variants across different types of samples, including highly attenuating specimens or those with hierarchical structures that vary rapidly in space.

V conclusion

We presented GF-MIST, a general and noise-robust MIST algorithm for SBXI, which uses a gradient-flow approach to solving the X-ray Fokker–Planck equation. This approach was applied to speckle data of a knotted carbon fibre to successfully recover the attenuation, phase-shift and dark-field images, in addition to the refraction angles from the sample. We compared the image quality of the dark-field and phase-shift images reconstructed from GF-MIST to the most recent MIST method, which we call ‘RV-MIST’ [Alloo_Morgan_Paganin_Pavlov_2023]. GF-MIST reconstructs dark field images at higher quality, at the cost of a slight decrease in the phase-shift image quality. GF-MIST also requires fewer speckle image pairs than RV-MIST to operate; therefore is suitable for applications where experimental simplicity and image acquisition speed are required.

Acknowledgements.

This research was funded by the French Agence Nationale de la Recherche and France 2030, grants MIAI Cluster (ANR-23-IACL-0006) and Labex CAMI (ANR-11-LABX-0004). Samantha J. Alloo acknowledges funding support from an AINSE Ltd. Early Career Researcher Grant (ECRG). We acknowledge the beamline scientists at the Imaging and Medical Beamline at the Australian Synchrotron for their support in acquiring the imaging data.