# Multimodal Intrinsic Speckle-Tracking (MIST) to extract rapidly-varying
diffuse X-ray scatter

Samantha J Alloo*,
School of Physical and Chemical Sciences, University of Canterbury,
Christchurch, New Zealand
<EMAIL_ADDRESS>
&Kaye S Morgan,
School of Physics and Astronomy, Monash University, Victoria, Australia
&David M Paganin,
School of Physics and Astronomy, Monash University, Victoria, Australia
&Konstantin M Pavlov,
School of Physical and Chemical Sciences, University of Canterbury,
Christchurch, New Zealand.
School of Physics and Astronomy, Monash University, Victoria, Australia.
School of Science and Technology, University of New England, Armidale,
Australia

###### Abstract

Speckle-based phase-contrast X-ray imaging (SB-PCXI) can reconstruct high-
resolution images of weakly-attenuating materials that would otherwise be
indistinguishable in conventional attenuation-based imaging. The experimental
setup of SB-PCXI requires only a sufficiently coherent source and spatially
random mask, positioned between the source and detector. The technique can
extract sample information at length scales smaller than the imaging system’s
spatial resolution; this enables multimodal signal reconstruction. “Multimodal
Intrinsic Speckle-Tracking” (MIST) is a rapid and deterministic formalism
derived from the paraxial-optics form of the Fokker-Planck equation. MIST
simultaneously extracts attenuation, refraction, and small-angle scattering
(diffusive-dark-field) signals from a sample and is more computationally
efficient compared to alternative speckle-tracking approaches. Hitherto,
variants of MIST have assumed the diffusive-dark-field signal to be spatially
slowly varying. Although successful, these approaches have been unable to
well-describe unresolved sample microstructure whose statistical form is not
spatially slowly varying. Here, we extend the MIST formalism such that there
is no such restriction, in terms of a sample’s rotationally-isotropic
diffusive-dark-field signal. We reconstruct multimodal signals of two samples,
each with distinct X-ray attenuation and scattering properties. The
reconstructed diffusive-dark-field signals have superior image quality
compared to our previous approaches which assume the diffusive-dark-field to
be a slowly varying function of transverse position. Our generalisation may
assist increased adoption of SB-PCXI in applications such as engineering and
biomedical disciplines, forestry, and palaeontology, and is anticipated to aid
the development of speckle-based diffusive-dark-field tensor tomography.

_K_ eywords phase-contrast X-ray imaging $\cdot$ speckle-based $\cdot$ dark-
field

##  1 Introduction

Phase-contrast X-ray imaging (PCXI)[1] techniques can reconstruct high-
resolution images of weakly-attenuating materials that would otherwise be
indistinguishable in conventional attenuation-based imaging, for example,
resolving soft-tissues in conventional radiography. PCXI techniques utilise
X-ray phase-shifts and attenuation, described by the material’s real and
imaginary components of its refractive index,
${n({\mathbf{r}}^{\prime})=1-\delta({\mathbf{r}}^{\prime})+i\beta({\mathbf{r}}^{\prime})}$,
respectively, to reconstruct images. Here ${\mathbf{r}}^{\prime}$ is the
three-dimensional position vector. This is a vital point of difference
compared to attenuation-based imaging, which only exploits X-ray attenuation.
The most common approach to converting sample-imposed phase effects into
measurable detector intensity differences is to introduce sophisticated
optical elements into the experimental setup, for example, via grating-
interferometry[2], analyser-based imaging[3], grid-based imaging[4, 5], and
edge-illumination imaging[6]. To eliminate the need for complex optics,
propagation-based[7, 8] PCXI has also been developed, however, it does require
sufficiently large sample-to-detector distances and has stringent coherence
requirements. The superiority of PCXI techniques over attenuation-based
imaging has already been proven in multiple disciplines. Such studies have
been performed using synchrotron and conventional laboratory X-ray sources for
applications in several fields, including biomedical[9, 10, 11, 12],
agricultural and food sciences[13], palaeontology[14, 15], and materials
science[16, 17].

PCXI techniques can also provide information regarding sample structures at
length scales smaller than the spatial resolution of the imaging system,
namely, through small-angle X-ray scattering (SAXS) (or diffusive-dark-field
(DF) imaging) from unresolved sample microstructure. PCXI techniques that
reconstruct the coherent, phase-contrast (PC), and diffuse, DF, flows of an
X-ray wavefield are often termed “multimodal”, in the sense that they can
extract multiple signals. Grating interferometry[18], analyser-based
imaging[19, 20], grid-based imaging[21], edge-illumination imaging[22], and
propagation-based imaging[23, 24] are examples of PCXI techniques that are
sensitive to SAXS, and can therefore reconstruct a sample’s DF signal.

Speckle-based (SB) PCXI, first developed in 2012[25, 26], is a particularly
appealing multimodal technique as it is experimentally simple, cost-effective,
and radiation dose-efficient. A speckle pattern is generated when a
sufficiently coherent wavefield propagates through a membrane with random
refractive index fluctuations[27]. A speckle pattern, in the context of SB-
PCXI, is used as an X-ray wavefront marker whose subsequent modification is
used to measure sample-induced speckle modulations, e.g., transverse spatial
shifts, attenuation, and blurring. The SB-PCXI experimental setup consists of
an X-ray source, speckle-generating mask, sample, and detector system,
positioned some finite distance downstream from the sample, as shown in Fig.
1. In this work, we consider a synchrotron X-ray source that produces a
monochromatic, paraxial wavefield having a high degree of both spatial and
temporal coherence. The X-ray wavefield is randomly modulated by propagation
through the speckle-generating mask (e.g., conventional sandpaper) and is then
registered by the position-sensitive detector positioned downstream. We
highlight two particularly attractive features of SB-PCXI. Firstly, the
speckle-generating mask can be any spatially-random medium. This removes the
restriction of needing precisely-manufactured optical elements, making the
experimental setup both easy to implement and flexible. Secondly, SB-PCXI
requires relatively low spatial and temporal coherence[28].

The inverse problem[29] of SB-PCXI involves reconstructing a sample’s
multimodal signals, given suitable reference-speckle and sample-reference-
speckle intensity images. The reference-speckle images resemble the
composition of the mask and the sample-reference-speckle images are captured
when a sample is placed into this reference-speckle field. This reference-
speckle pattern is modified, depending on the refractive properties of the
sample, and these speckle modifications are used to reconstruct sample
information. Transverse speckle shifts are associated with the PC signal,
whereas speckle blurring, or reduction in visibility, is correlated to the
sample’s DF signal. There are two distinct approaches in the research
literature to solve the multimodal inverse SB-PCXI problem, namely, extrinsic
and intrinsic approaches. X-ray speckle-vector tracking[25] (XSVT), mixed XSVT
approaches[30, 31], X-ray speckle-scanning[32, 33] (XSS), and unified
modulated pattern analysis[34] (UMPA) are examples of extrinsic approaches
that reconstruct multimodal signals using iterative pixel-wise methods on SB-
PCXI data acquired at multiple random-mask positions. In essence, XSVT and
mixed XSVT approaches perform a zero-normalised cross-correlation analysis,
between the reference-speckle and sample-reference-speckle images. UMPA
proposes a different methodology, based on least-squares minimisation between
a model and the measurement of the sample-reference-speckle pattern across all
mask positions. These techniques have the advantage of a relatively short
image acquisition time compared to XSS[32, 33] techniques. The SB-PCXI inverse
problem was reconceptualised in 2018 when Paganin et al.[35] proposed a
geometric-flow approach to reconstruct a sample’s PC signal using a single set
of SB-PCXI intensity data; this was the first realisation of so-called
intrinsic speckle-tracking, since it does not explicitly track individual
speckles, but rather solves a partial differential equation formulated at the
whole-of-image level. This speckle-tracking geometric-flow formalism[35] was
then combined with a Fokker-Planck-type[36, 37] generalisation of the
transport-of-intensity equation[38] of paraxial wave optics to allow for
multimodal intrinsic signal extraction[39], which was named “Multimodal
Intrinsic Speckle-Tracking” (MIST). The Fokker-Planck[40] expression is based
on local energy conservation, and it considers transverse radiation flows as a
combination of coherent and diffusive effects; these effects for the case of
multimodal signal extraction align with the PC and DF signals, respectively.
MIST is less computationally expensive than the alternative extrinsic
approaches[25, 26, 34, 30, 31, 32, 33]. This means that MIST can be used to
rapidly reconstruct tomographic data[41] of both the PC and DF signals,
therefore, to provide three-dimensional sample information that is
inaccessible in single-projection imaging.

MIST was first developed[39] under three key assumptions: (a) the sample is a
pure phase-object, such that X-ray attenuation can be neglected, (b) the
unresolved sample microstructure diffusely scatters the X-ray beam in a
rotationally-isotropic manner, and (c) the sample’s DF signal is a slowly
varying function of transverse position. The first assumption, (a), was
relaxed in our most recent work where we presented the case of attenuating
materials having rotationally-isotropic position-dependent diffuse
scatter[41]. Assumption (b) was relaxed when we considered anisotropically-
scattering attenuating materials[42], such that directional DF [43, 44]
signals could be obtained. Here, we should highlight that we have also
performed DF computed-tomography (CT) using MIST[41], as the reconstructed
two-dimensional DF signals had a sufficiently high signal-to-noise ratio (SNR)
and spatial resolution such that they could be used in standard CT
reconstruction algorithms, such as filtered back-projection[45]. Assumption
(c), which considers the sample’s DF signal to be spatially slowly varying, is
a condition that has remained in all of the MIST approaches to date. Although
both Pavlov et al.[39] and Alloo et al.[41] have demonstrated that this
approach to MIST is capable of reconstructing images with a high spatial
resolution, which is at least comparable to the alternative extrinsic
approaches, our recent investigations have found that this approximation may
break down at sharp interfaces. This assumption may hinder the full potential
of previous MIST approaches. In particular, samples that have rapidly varying
microstructure-autocorrelation functions may breach the domain of
applicability of all previous MIST approaches. In the present paper, we
generalise MIST to alleviate this restriction, thereby broadening its domain
of utility. Several studies have verified the broad applicability and
importance of DF imaging[46, 33, 47], with a significant focus on biomedical
clinical applications[48], e.g., using DF imaging for early-stage diagnosis of
lung diseases such as fibrosis[49, 50], pneumothorax[51], emphysema[52], and
breast cancer[23]. Our improved MIST approach could provide an alternative
experimentally versatile, low-dose imaging technique that can reconstruct
high-resolution multimodal signals in two- and three-dimensions.

This paper progresses as follows. First, we theoretically develop the new
generalised MIST approach, deriving an analytical solution for a sample’s
phase-shift (PC signal) and effective diffusion coefficient (DF signal).
Numerical stabilisation techniques are then discussed, before applying the
approach to synchrotron SB-PCXI data of two samples that have different X-ray
attenuation and SAXS characteristics. Moreover, we first consider a sample
that is weakly-attenuating (almost a pure phase-object), and then a more
attenuating object to investigate the breadth of applicability of our
approach. We then compare this new approach, qualitatively and quantitatively,
to both of the published rotationally-isotropic MIST approaches[39, 41]: (a)
Pavlov et al.’s approach[39], which neglects X-ray attenuation, and (b) Alloo
et al.’s approach[41], which considers it. Note that both previous approaches
approximate the DF to be spatially slowly varying. The paper is finished by
discussing potential future research avenues.

Figure 1: Experimental setup of speckle-based phase-contrast X-ray imaging
using a synchrotron X-ray source.

## Theoretical derivation of updated MIST approach

The SB-PCXI form of the Fokker-Planck equation models the coherent and
diffusive flows of X-rays in the case of a phase-object, as described by its
phase-shift, $\phi_{\textrm{ob}}(\textbf{r})$, and effective diffusion
coefficient, $D_{\textrm{eff, Phase}}(\textbf{r};\Delta)$, respectively. The
SB-PCXI Fokker-Planck equation is derived assuming the reference-speckle field
is spatially well-resolved, and that this field, alongside the speckle field
in the presence of the sample, obeys a Fokker-Planck[40] extension of the
geometric-flow formalism for speckle tracking[35]. The Fokker-Planck equation
in this instance is[36]

$I_{R}(\textbf{r})-I_{S}(\textbf{r})=\frac{\Delta}{k}\nabla_{\perp}\cdot\left[I_{R}(\textbf{r})\nabla_{\perp}\phi_{\textrm{ob}}(\textbf{r})\right]-\Delta\nabla_{\perp}^{2}\left[D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)I_{R}(\textbf{r})\right],$ (1)

where $I_{R}(\textbf{r})$ and $I_{S}(\textbf{r})$ denote the reference-speckle
and sample-reference-speckle intensities, respectively,
$\textbf{r}\equiv(x,y)$ denotes Cartesian coordinates in planes perpendicular
to the optical axis $z$, assuming plane-wave illumination, $\Delta$ is the
sample-to-detector distance, $k$ is the wavenumber,
$\nabla_{\perp}=(\partial/\partial x,\partial/\partial y)$ is the transverse
gradient operator in the $(x,y)$ plane, and $\nabla_{\perp}^{2}$ is the
transverse Laplacian operator. The reference-speckle intensity,
$I_{R}(\textbf{r})$, is obtained with the speckle mask in and the object out
of the X-ray beam, and the sample-reference-speckle intensity,
$I_{S}(\textbf{r})$, is obtained by placing a phase- and SAXS-inducing object
into the speckle-modulated X-ray beam. In previous rotationally-isotropic MIST
approaches[39, 41], the Fokker-Planck equation has been simplified (for phase-
objects) to a linear equation in terms of the effective diffusion coefficient
and the Laplacian of the phase-shift. This simplification resulted by assuming
the effective diffusion coefficient to be spatially slowly varying, such that
various second-order terms could be neglected. In the present work, we do not
make this assumption, and instead, the effective diffusion coefficient can
accurately describe the diffuse scattering signal from unresolved sample
microstructure for which the autocorrelation function can be a rapidly varying
function of transverse position.

The present approach begins by expanding the coherent flow term, namely the
first term on the right side of the preceding equation, into its two
components that describe the lensing and prism-like effects[36]:

$\frac{\Delta}{k}\nabla_{\perp}\cdot\left[I_{R}(\textbf{r})\nabla_{\perp}\phi_{\textrm{ob}}(\textbf{r})\right]=\frac{\Delta}{k}\left[I_{R}(\textbf{r})\nabla_{\perp}^{2}\phi_{\textrm{ob}}(\textbf{r})+\nabla_{\perp}I_{R}(\textbf{r})\cdot\nabla_{\perp}\phi_{\textrm{ob}}(\textbf{r})\right].$
(2)

Following the approximation described in Pavlov et al.[39], we neglect the
scalar product of the gradient of the random rapidly varying wavefield
intensity $I_{R}(\textbf{r})$ with the more slowly changing gradient of the
wavefield phase $\phi_{\textrm{ob}}(\textbf{r})$, namely, we assume that

$\frac{\Delta}{k}\nabla_{\perp}\cdot\left[I_{R}(\textbf{r})\nabla_{\perp}\phi_{\textrm{ob}}(\textbf{r})\right]\approx\frac{\Delta}{k}I_{R}(\textbf{r})\nabla_{\perp}^{2}\phi_{\textrm{ob}}(\textbf{r}).$
(3)

Hence equation (1) becomes (cf. Ref. 39)

$I_{R}(\textbf{r})-I_{S}(\textbf{r})=\frac{\Delta}{k}I_{R}(\textbf{r})\nabla_{\perp}^{2}\phi_{\textrm{ob}}(\textbf{r})-\Delta\nabla_{\perp}^{2}\left[D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)I_{R}(\textbf{r})\right].$ (4)

The associated inverse problem of multimodal SB-PCXI, in this work, involves
solving equation (4) for $D_{\textrm{eff, Phase}}(\textbf{r};\Delta)$ and
$\phi_{\textrm{ob}}(\textbf{r})$, given the SB-PCXI intensity data,
$I_{R}(\textbf{r})$ and $I_{S}(\textbf{r})$. We expand the final term on the
right side of equation (4), to give the following alternative, but still
equivalent, form of a phase-object’s SB-PCXI Fokker-Planck expression:

$\displaystyle\frac{1}{\Delta}\left[I_{R}(\textbf{r})-I_{S}(\textbf{r})\right]=I_{R}(\textbf{r})\nabla_{\perp}^{2}\left[\frac{1}{k}\phi_{\textrm{ob}}(\textbf{r})-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\right]-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\nabla_{\perp}^{2}I_{R}(\textbf{r})$ (5)
$\displaystyle-2D_{\textrm{eff,
Phase}}^{x}(\textbf{r};\Delta)I_{R}^{x}(\textbf{r})-2D_{\textrm{eff,
Phase}}^{y}(\textbf{r};\Delta)I_{R}^{y}(\textbf{r}).$

This is a linear equation in terms of four unknowns—namely,
$\nabla_{\perp}^{2}\left[\frac{1}{k}\phi_{\textrm{ob}}(\textbf{r})-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\right]$, $D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)$, $D_{\textrm{eff, Phase}}^{x}(\textbf{r};\Delta)$,
and $D_{\textrm{eff, Phase}}^{y}(\textbf{r};\Delta)$—where the partial
derivatives in the spatial coordinates x and y are denoted with superscripts x
and y, respectively. These four functions can then be employed to reconstruct
the sample’s true phase-shifts and effective diffusion coefficient. To obtain
unique solutions for four unknown variables we require four equations; to do
this, four unique forms of equation (5) can be generated by taking four
independent measurements of $I_{R}(\textbf{r})$ and $I_{S}(\textbf{r})$. This
can be achieved by, for example, transversely shifting the mask to generate a
new reference-speckle pattern. Four independent measurements, denoted with
subscripts $a,b$, $c$, and $d$, will induce coherent and diffuse X-ray flows
that satisfy the following system of equations:

$\displaystyle\frac{1}{\Delta}\left[I_{R_{a}}(\textbf{r})-I_{S_{a}}(\textbf{r})\right]=I_{R_{a}}(\textbf{r})\nabla_{\perp}^{2}\left[\frac{1}{k}\phi_{\textrm{ob}}(\textbf{r})-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\right]-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\nabla_{\perp}^{2}I_{R_{a}}(\textbf{r})$
$\displaystyle-2D_{\textrm{eff,
Phase}}^{x}(\textbf{r};\Delta)I_{R_{a}}^{x}(\textbf{r})-2D_{\textrm{eff,
Phase}}^{y}(\textbf{r};\Delta)I_{R_{a}}^{y}(\textbf{r}),$ (6)
$\displaystyle\frac{1}{\Delta}\left[I_{R_{b}}(\textbf{r})-I_{S_{b}}(\textbf{r})\right]=I_{R_{b}}(\textbf{r})\nabla_{\perp}^{2}\left[\frac{1}{k}\phi_{\textrm{ob}}(\textbf{r})-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\right]-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\nabla_{\perp}^{2}I_{R_{b}}(\textbf{r})$
$\displaystyle-2D_{\textrm{eff,
Phase}}^{x}(\textbf{r};\Delta)I_{R_{b}}^{x}(\textbf{r})-2D_{\textrm{eff,
Phase}}^{y}(\textbf{r};\Delta)I_{R_{b}}^{y}(\textbf{r}),$ (7)
$\displaystyle\frac{1}{\Delta}\left[I_{R_{c}}(\textbf{r})-I_{S_{c}}(\textbf{r})\right]=I_{R_{c}}(\textbf{r})\nabla_{\perp}^{2}\left[\frac{1}{k}\phi_{\textrm{ob}}(\textbf{r})-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\right]-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\nabla_{\perp}^{2}I_{R_{c}}(\textbf{r})$
$\displaystyle-2D_{\textrm{eff,
Phase}}^{x}(\textbf{r};\Delta)I_{R_{c}}^{x}(\textbf{r})-2D_{\textrm{eff,
Phase}}^{y}(\textbf{r};\Delta)I_{R_{c}}^{y}(\textbf{r}),$ (8)
$\displaystyle\frac{1}{\Delta}\left[I_{R_{d}}(\textbf{r})-I_{S_{d}}(\textbf{r})\right]=I_{R_{d}}(\textbf{r})\nabla_{\perp}^{2}\left[\frac{1}{k}\phi_{\textrm{ob}}(\textbf{r})-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\right]-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\nabla_{\perp}^{2}I_{R_{d}}(\textbf{r})$
$\displaystyle-2D_{\textrm{eff,
Phase}}^{x}(\textbf{r};\Delta)I_{R_{d}}^{x}(\textbf{r})-2D_{\textrm{eff,
Phase}}^{y}(\textbf{r};\Delta)I_{R_{d}}^{y}(\textbf{r}).$ (9)

The phase-shift $\phi_{\textrm{ob}}(\textbf{r})$, together with the effective
diffusion coefficient $D_{\textrm{eff, Phase}}(\textbf{r};\Delta)$ and its
transverse derivatives, are exclusively defined by the object and are hence
equal for equations (Theoretical derivation of updated MIST
approach)–(Theoretical derivation of updated MIST approach). Equations
(Theoretical derivation of updated MIST approach)–(Theoretical derivation of
updated MIST approach) can be solved using, for example, a standard Gaussian-
elimination. Here, it is important to highlight that previous MIST approaches
for rotationally-isotropic diffuse scatter [39, 41] require SB-PCXI data from
just two mask positions. Within the present method, we require four mask
positions since we extract two additional quantities, namely $D_{\textrm{eff,
Phase}}^{x}(\textbf{r};\Delta)$ and $D_{\textrm{eff,
Phase}}^{y}(\textbf{r};\Delta)$. These terms are small for samples whose
unresolved microstructure can be considered to have an autocorrelation that is
a slowly varying function of transverse position. However, these terms are
significant in the contrary case where the unresolved sample-microstructure-
autocorrelation is rapidly varying. These terms are also important to
correctly reconstruct sample edges.

A sample’s true effective diffusion coefficient can then be reconstructed by
aggregating the three extracted quantities, $D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)$, $D_{\textrm{eff, Phase}}^{x}(\textbf{r};\Delta)$,
and $D_{\textrm{eff, Phase}}^{y}(\textbf{r};\Delta)$. We employ a technique
adopted in various differential imaging techniques[53, 54, 55, 5] whereby a
two-dimensional function, $g(x,y)$, can be calculated from its two spatial
derivatives, $g^{x}(x,y)$ and $g^{y}(x,y)$, using the identity:

$g(x,y)=\mathcal{F}^{-1}\left[\frac{\mathcal{F}\left(g^{x}(x,y)+ig^{y}(x,y)\right)}{ik_{x}-k_{y}}\right].$
(10)

Here, $\mathcal{F}$ denotes the two-dimensional Fourier transformation with
respect to x and y, with $k_{x}$ and $k_{y}$ being the corresponding Fourier-
space variables, respectively. This expression is unstable near the Fourier-
space origin, $(k_{x},k_{y})=(0,0)$, hence the solution will diverge at low
spatial frequencies unless sufficiently regularised. An effective diffusion
coefficient can be calculated by equation (10), using the two calculated
spatial derivatives $D_{\textrm{eff, Phase}}^{x}(\textbf{r};\Delta)$ and
$D_{\textrm{eff, Phase}}^{y}(\textbf{r};\Delta)$. This solution will be stable
at high spatial frequencies, that is, far from the Fourier-space origin, and
the opposite is true for the effective diffusion coefficient $D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)$ that is extracted directly from the system of
linear equations (equations (Theoretical derivation of updated MIST
approach)–(Theoretical derivation of updated MIST approach)). These two
solutions can be combined to consider only the stable spatial frequencies in
each DF solution. In particular, a true effective diffusion coefficient,
$D_{\textrm{eff, Phase}}^{\textrm{True}}(\textbf{r};\Delta)$, can be
calculated by implementing the following Fourier-space weighted filtering (cf.
Ref. 56), where $\rho$ is a cut-off parameter:

$D_{\textrm{eff,
Phase}}^{\textrm{True}}(\textbf{r};\Delta)=\mathcal{F}^{-1}\left[e^{-\rho(k_{x}^{2}+k_{y}^{2})}\mathcal{F}\left(D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\right)+\frac{1-e^{-\rho(k_{x}^{2}+k_{y}^{2})}}{ik_{x}-k_{y}}\mathcal{F}\left(D_{\textrm{eff,
Phase}}^{x}(\textbf{r};\Delta)+iD_{\textrm{eff,
Phase}}^{y}(\textbf{r};\Delta)\right)\right].$ (11)

This completes our description of the method for reconstructing a phase-
object’s true effective diffusion coefficient. The sample-induced phase-shift
term, $\phi_{\textrm{ob}}(\textbf{r})$, can then be reconstructed by utilising
equation (4), via

$\phi_{\textrm{ob}}(\textbf{r})=\nabla_{\perp}^{-2}\left[\frac{k}{\Delta
I_{R}(\textbf{r})}\left(I_{R}(\textbf{r})-I_{S}(\textbf{r})+\Delta\nabla_{\perp}^{2}\left[D_{\textrm{eff,
Phase}}^{\textrm{True}}(\textbf{r};\Delta)I_{R}(\textbf{r})\right]\right)\right].$
(12)

In the above expression,

$\nabla_{\perp}^{-2}=-\mathcal{F}^{-1}\frac{1}{k_{x}^{2}+k_{y}^{2}}\mathcal{F}$
(13)

is the inverse Laplacian operator, derived using the two-dimensional Fourier
derivative theorem[1]. The above PC signal extraction is more numerically
stable than utilising the
$\nabla_{\perp}^{2}\left[\frac{1}{k}\phi_{\textrm{ob}}(\textbf{r})-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\right]_{\textrm{Recon}}$ term reconstructed from QR
decomposition[57], as that solution would inherently suffer from the numerical
instabilities associated with $D_{\textrm{eff, Phase}}(\textbf{r};\Delta)$,
$D_{\textrm{eff, Phase}}^{x}(\textbf{r};\Delta)$ and $D_{\textrm{eff,
Phase}}^{y}(\textbf{r};\Delta)$, whereas equation (12) considers the
stabilised $D_{\textrm{eff, Phase}}^{\textrm{True}}(\textbf{r};\Delta)$.

Up until this point, X-ray attenuation by the sample has been neglected. We
now extend our analysis to the case of a weakly-attenuating object and
calculate its effective diffusion coefficient. We do this based on a
relationship obtained in Alloo et al.[41] (see equations (24) and (18)
therein) in which an attenuating-object’s effective diffusion coefficient,
$D_{\textrm{eff, Atten}}^{\textrm{True}}(\textbf{r};\Delta)$, can be
calculated from the phase-object approximation, $D_{\textrm{eff,
Phase}}^{\textrm{True}}(\textbf{r};\Delta)$, using

$D_{\textrm{eff,
Atten}}^{\textrm{True}}(\textbf{r};\Delta)=\frac{D_{\textrm{eff,
Phase}}^{\textrm{True}}(\textbf{r};\Delta)}{I_{\textrm{ob}}(\textbf{r})}.$
(14)

Above, $I_{\textrm{ob}}(\textbf{r})$ is the object’s attenuation term
describing the intensity at the exit surface of the sample, $z=0$, after the
object has attenuated the incident X-ray beam of unit intensity. To calculate
$I_{\textrm{ob}}(\textbf{r})$, we consider a single-material object such that
the projection approximation[1] can be written as
$\phi_{\textrm{ob}}(\textbf{r})=-k\delta t(\textbf{r})$ and
$I_{\textrm{ob}}(\textbf{r})=\textrm{exp}\left[-2k\beta t(\textbf{r})\right]$,
where $t(\textbf{r})$ is the projected thickness of the object along the
direction $z$ of the X-rays[58]. Hence

$I_{\textrm{ob}}(\textbf{r})=\textrm{exp}\left[\frac{2\phi_{\textrm{ob}}(\textbf{r})}{\gamma}\right],$
(15)

where $\gamma=\delta/\beta$ for the single-material object. This attenuation
term, obtained using the phase-shift term from equation (12), can then be used
in equation (14) to reconstruct $D_{\textrm{eff,
Atten}}^{\textrm{True}}(\textbf{r};\Delta)$. Although in theory this
attenuation extraction is restricted to single-material objects, it can be
extended to multi-material objects by taking the difference in the refractive
index components[59, 60] for composite materials. Furthermore, in a
tomographic context, it has been proven that this approximation does not
affect the reconstructed attenuation coefficient,
$\beta({\mathbf{r}}^{\prime})$, far away from material interfaces[61], and
hence, this restriction would only be adverse in a sample that has several
composite materials with significantly differing attenuation and refraction
properties.

## Stabilising the SB-PCXI multimodal inverse-problem

The inverse-problem[29] of reconstructing the effective diffusion coefficient
requires appropriate numerical regularisation. Here, we apply the common
approach[62, 63] of a Tikhonov regularisation[64], which sufficiently
stabilises the signal reconstruction. In its simplest form, a Tikhonov
regularisation of the quotient of two functions $A$ and $B$ can be employed,
using

$\frac{A}{B}\rightarrow\frac{AB}{B^{2}+\alpha},$ (16)

where $\alpha\geq 0$ is a regularisation parameter whose magnitude is
sufficiently small compared to $B^{2}$.

Within this paper, the effective diffusion coefficient is reconstructed by
solving a full-rank system of four linear equations, equations (Theoretical
derivation of updated MIST approach)–(Theoretical derivation of updated MIST
approach). The numerical stability of this solution can be improved by
utilising SB-PCXI data from more mask positions. Namely, for N mask positions,
we can generate a system of N linear equations following the form of equation
(5). The system generated can then be solved in a least-squares sense using
pixel-wise QR decomposition[57]. Tikhonov’s regularisation method can also be
applied to an ill-posed least-squares problem (QR factorisation), as described
by Zhu[65]. Namely, rather than using QR decomposition to solve the linear
system $A\bar{x}=\bar{b}$ for the least-squares solution $\tilde{x}$, QR
decomposition can instead be performed on the system $\begin{pmatrix}A;\alpha
I\end{pmatrix}\bar{x}=\begin{pmatrix}\bar{b};0\end{pmatrix},$ where $M\times
N$ is the coefficient matrix $A$, “;” denotes a new row, $I$ is the $N\times
N$ identity matrix, $\alpha$ is the chosen regularisation parameter, and the
right-hand-side vector $\bar{b}$ is filled with zeroes to reach the size of
$(M+N)\times 1$. The described Tikhonov-regularised QR decomposition can be
used to solve the overdetermined system of linear equations for the four
unknown variables,
$\nabla_{\perp}^{2}\left[\frac{1}{k}\phi_{\textrm{ob}}(\textbf{r})-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\right]$, $D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)$, $D_{\textrm{eff, Phase}}^{x}(\textbf{r};\Delta)$,
and $D_{\textrm{eff, Phase}}^{y}(\textbf{r};\Delta)$. For reasons described in
the preceding text, equation (11) is numerically stable and therefore does not
need to be regularised. The subsequent phase extraction, that is equation
(12), is ill-posed close to the Fourier-space origin $(k_{x},k_{y})=(0,0)$ and
hence an appropriate Tikhonov regularisation should be applied following
equation (16). For severely ill-posed cases, the phase extraction can be
further stabilised by utilising instances of equation (4), substituting in the
true effective diffusion coefficient using the method described above, and
performing Tikhonov-regularised QR decomposition to solve for
$\nabla_{\perp}^{-2}\phi_{\textrm{ob}}(\textbf{r})$ before applying the
inverse Laplacian operator to reconstruct $\phi_{\textrm{ob}}(\textbf{r})$.

For the Tikhonov regularisation to operate successfully, the regularisation
parameter, $\alpha$, needs to be selected appropriately for the given data. If
$\alpha$ is chosen to be too large, the computed solution will be over-
smoothed and will therefore lack fine detail. In the case of image
reconstruction, this means the computed solution will have poor spatial
resolution, although a high SNR. Conversely, if $\alpha$ is chosen to be too
small, the computed solution will be severely contaminated with errors
resulting from numerical instabilities. Evidently, optimising the
regularisation parameter is critical to successfully extract multimodal
signals. There are algorithms in the current research literature that optimise
the Tikhonov regularisation parameter for a given ill-posed problem, see Park
et al.[66] and references therein. However, in image reconstruction, image
quality metrics can be used to determine the optimal regularisation parameter.
In the present work, we used four metrics:

  1. 1.

Naturalness Image Quality Evaluator (NIQE)[67]: The NIQE is a blind image
quality assessment based on an image’s measurable deviations from statistical
regularities observed in natural images. A lower NIQE reflects an image with a
higher perceived image quality.

  2. 2.

Azimuthally Averaged Power-Spectrum[68]: The two-dimensional power-spectrum of
an image can be calculated by taking the absolute square of the Fourier-
transformed image. This can then be azimuthally-averaged, with the centre at
the Fourier-space origin $(k_{x},k_{y})=(0,0)$, to calculate a one-dimensional
power-spectrum that shows the contribution of all spatial frequencies in an
image. The noise in the image is reflected by the so-called “noise-floor”,
which typically makes up the majority of the signal at high spatial
frequencies, and the spatial resolution can be gauged by the “knee” of the
power-spectrum, namely the frequency at which noise becomes a significant
contribution to the signal.

  3. 3.

SNR: The SNR measures the magnitude of a signal relative to background noise.
It is used to quantify signal quality in an image, and is defined as

$SNR=I_{\textrm{avg}}/\sigma,$ (17)

where $I_{\textrm{avg}}$ is the signal strength, which can be measured as the
average pixel value within a region of approximately uniform signal, and
$\sigma$ is the noise as measured by the standard deviation of pixel values.
Note that if the noise characteristics are the same inside and outside an
object, then it may be easiest to measure $\sigma$ in the region outside the
object to avoid variations in the signal that come with a complex object.

  4. 4.

Human visual perception: Although subjective, meaningful image quality
measurements can be made by a human observer’s visual assessment of an image.

## Applying the theoretical approach to synchrotron SB-PCXI data

To test the proposed MIST approach, we extracted multimodal signals from two
samples, a wattle flower (WF) and red currant (RC). These two samples had
different X-ray attenuation coefficients, $\beta(\textbf{r}^{\prime})$, and
thicknesses, and hence attenuate the X-ray beam differently. In particular,
the WF was weakly-attenuating (almost a pure phase-object) and the RC was non-
negligibly-attenuating.

### Experimental procedures

SB-PCXI data of the WF were collected in experimental hutch 3B of the Imaging
and Medical Beamline (IMBL) at the Australian Synchrotron, similar to the
setup shown in Fig. 1. The entrance window to hutch 3B was located 135 m from
the source. A virtually monochromatic 25 keV X-ray beam was used for imaging.
The Ruby detector[69], which has a single pco.edge sensor and lens-coupled
scintillator, was positioned $\Delta=2$ m downstream from the sample. The
pixel array was 2560 $\times$ 2160 and a 105 mm macro lens was used to achieve
an effective pixel size of 9.9 $\mu$m. The speckle-generating mask was located
around 60 cm upstream of the sample, and a combination of sandpapers of grit
P40, P80 and P120 were simultaneously used. A combination of sandpapers was
used in this case such that the generated speckle pattern had a range of high-
visibility feature sizes. The speckle-generating sandpaper was placed on a
translation motor stage such that it could be brought in and out of the X-ray
beam path and also translated transversely to allow the beam to pass through
different parts of the mask. The effective speckle-size was 136.4 $\mu$m, as
measured by the average full-width at half-maximum of the autocorrelation
function[27] in the horizontal, x, and vertical, y, directions of the
reference-speckle field.

The RC sample was imaged using SB-PCXI at the European Synchrotron Radiation
Facility (ESRF) beamline BM05; these data were obtained and originally
published by Berujon and Ziegler[30]. The setup was similar to that shown in
Fig. 1, with an X-ray energy of 17 keV and spectral bandwidth of $\Delta
E/E\approx 10^{-4}$, which was produced using a double crystal Si(111)
monochromator located 27 m from the X-ray source. The RC was placed on a stage
55 m from the source. The detector system consisted of a Fast Read-Out Low-
Noise (FReLoN) e2V camera coupled to an optical imaging thin scintillator[70,
71]. This detector was placed $\Delta=1$ m downstream from the sample and the
effective pixel size of the optical system was 5.8 $\mu$m. The speckle-
generating sandpaper with grit size P800 was placed 0.5 m upstream from the
sample and had an effective speckle-size of 20.4 $\mu$m.

Image acquisitions for all data were collected using a similar procedure for
both samples: dark-current (no X-ray beam) and flat-field (sample and mask not
in the beam) exposures were collected before and after the scan, and
reference-speckle images with only the speckle mask in the beam were collected
before and after the scan. The speckle mask was transversely shifted in the
$x$-direction perpendicular to the optical axis (see Fig. 1) to acquire
multiple unique sets of SB-PCXI data. Seven and fifteen sets of SB-PCXI data
were collected for the RC and WF, respectively. The collected SB-PCXI data
were then processed using a Python3 script to implement our multimodal signal
extraction algorithm. An open-access repository for this code is on
GitHub[72].

### Multimodal signal extraction

Multimodal signals were extracted for the WF and RC using our new generalised
MIST approach. The entirety of the available SB-PCXI data for each sample,
that is the maximum number of masks, were used to reconstruct the multimodal
signals. Although fewer could be used, this work focuses solely on the new
theoretical development; a quantitative analysis of the influence of the
number of SB-PCXI data sets is given in Pavlov et al.[39] (see, in particular,
Fig. 2 in Ref. 39). The established system of linear equations for each sample
was solved using the Tikhonov-regularised QR decomposition described above. It
was found that the standard deviation of the coefficient matrix divided by
$10^{4}$ provided the optimal Tikhonov regularisation parameter for the case
of each sample. The phase-object approximation of the sample’s effective
diffusion coefficient, $D_{\textrm{eff, Phase}}(\textbf{r};\Delta)$, and its
spatial derivatives were calculated using this method, from which the true
phase-object approximation of the DF signal, $D_{\textrm{eff,
Phase}}^{\textrm{True}}(\textbf{r};\Delta)$, was computed. Next, the sample’s
phase-shifts, $\phi_{\textrm{ob}}(\textbf{r})$, attenuation term,
$I_{\textrm{ob}}(\textbf{r})$, and true attenuating-object effective diffusion
coefficient, $D_{\textrm{eff, Atten}}^{\textrm{True}}(\textbf{r};\Delta)$,
were calculated. It should be noted that the calculation of
$\phi_{\textrm{ob}}(\textbf{r})$ had to be regularised as equation (13) is
unstable at the Fourier-space origin; it was regularised using equation (16)
with $\alpha=0.0001$, which was suitable in all cases. The variable
$\gamma=\delta/\beta$ was required to calculate the object’s attenuation term,
as in equation (15). For the WF, the generic elemental composition of a plant
stem was used within the TS imaging calculator[73] to determine its complex
refractive index at 25 keV, producing $\gamma_{\,\textrm{WF, 25~{}keV}}=1403$.
For the RC sample, $\gamma$ was taken to be that of water at 17 keV, that is,
$\gamma_{\,\textrm{RC, 17~{}keV}}=\delta_{\,\textrm{water,
17~{}keV}}/\beta_{\,\textrm{water, 17~{}keV}}=1146$.

Multimodal signals were also calculated using our previously-published MIST
approaches, that is, equation (6) in Pavlov et al.[39] and equations (16)–(18)
in Alloo et al.[41], to provide a point of comparison for the new approach. As
described earlier, the approach published by Pavlov et al.[39] neglects X-ray
attenuation, but Alloo et al.[41] considers it. In both of these approaches,
the effective diffusion coefficient is assumed to be spatially slowly varying.
The multimodal signal extraction described by Alloo et al.[41] involved
numerical stabilisation via a Tikhonov-regularised “Weighted Determinant”
approach. Within this work, the multimodal signals extracted using Pavlov et
al.’s[39] algorithm were stabilised in an identical way to that described in
Alloo et al.[41], that is a Tikhonov-regularised “Weighted Determinant”
approach. The optimal regularisation parameter, in both instances, was equal
to the mean of the denominator in equation (23) of Ref. 41, divided by 100,
for both samples. Here we do not perform a comparison with our directional DF
approach[42].

### Weakly-attenuating sample: Wattle flower (WF)

We begin by investigating the WF sample as this sample conforms most closely
to the underlying assumptions of the derived theoretical approach, as it is a
weakly-attenuating object. Therefore, the multimodal signals (DF and PC)
should be superior using the present approach compared to those extracted
using our previous MIST approaches[39, 41]. Using the methodology described
above, an over-determined system of fifteen linear equations was solved using
Tikhonov-regularised QR decomposition, from which the true phase-object DF
approximation, $D_{\textrm{eff, Phase}}^{\textrm{True}}(\textbf{r};\Delta)$,
attenuation term, $I_{\textrm{ob}}(\textbf{r})$, and true attenuating-object
DF signals, $D_{\textrm{eff, Atten}}^{\textrm{True}}(\textbf{r};\Delta)$, were
calculated. The cut-off parameter, $\rho$, was determined by investigating the
SNR and NIQE of the reconstructed $D_{\textrm{eff,
Phase}}^{\textrm{True}}(\textbf{r};\Delta)$ for various cut-off parameter
values. Figure 2 shows how the value of the cut-off parameter influences the
image quality of the reconstructed DF, as measured by the NIQE and SNR of the
entire DF reconstruction. As an image with a better-perceived image quality
has a lower NIQE, Fig. 2 is shown as the reciprocal of the NIQE such that it
follows the same “bigger means better” convention as the SNR shown on the
secondary axis. The optimal cut-off parameter is 34 $\mu\textrm{m}{{}^{2}}$
and 198 $\mu\textrm{m}{{}^{2}}$ as measured by the SNR and reciprocal NIQE,
respectively. As we have described, the NIQE is a blind image quality metric
that measures the so-called perceived image quality. Within this work, the
NIQE was calculated using an untrained NIQE function, called niqe(A), in
MATLAB. This untrained NIQE model may mistake artefacts for signal. For
example, enhanced edges due to residual PC effects, or Fresnel fringes, are
considered qualitatively better by the NIQE as they have an increased
sharpness and contrast. Such features are present in the reconstructed DF
signal when the cut-off parameter is too large, and hence, explains why the
NIQEs for these images indicate higher perceived image quality. It follows
that the optimal cut-off parameter should be selected by appropriately
considering the SNR and NIQE simultaneously. We selected the optimal cut-off
parameter to be $\rho=$27 $\mu\textrm{m}{{}^{2}}$ as this balances the local
reciprocal NIQE maximum (at approximately $\rho=$18 $\mu\textrm{m}{{}^{2}}$),
global SNR maximum, and also the human observer’s verdict.

Figure 2: Influence of cut-off parameter, $\rho$, on the image quality of the
wattle flower’s reconstructed phase-object approximation of the effective
diffusion coefficient, $D_{\textrm{eff,
Phase}}^{\textrm{True}}(\textbf{r};\Delta)$. Image quality is measured by the
(blue) reciprocal of the Naturalness Image Quality Evaluator (NIQE)[67] and
(orange) signal-to-noise ratio (SNR). The dashed red vertical line denotes the
optimal cut-off parameter that appropriately considers both metrics and the
verdict of a human observer.  Figure 3:  Solutions of the system of linear
equations, equations (Theoretical derivation of updated MIST
approach)–(Theoretical derivation of updated MIST approach), and the wattle
flower’s reconstructed multimodal signals; (a) is the effective diffusion
coefficient, $D_{\textrm{eff, Phase}}(\textbf{r};\Delta)$, (b) is
$\nabla_{\perp}^{2}\left[\frac{1}{k}\phi_{\textrm{ob}}(\textbf{r})-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\right]_{\textrm{Recon}}$, (c) and (d) are the two
spatial derivatives of the effective diffusion coefficient, $D_{\textrm{eff,
Phase}}^{y}(\textbf{r};\Delta)$ and $D_{\textrm{eff,
Phase}}^{x}(\textbf{r};\Delta)$, respectively, (e) is the reconstructed
attenuation term, $I_{\textrm{ob}}(\textbf{r})$, and (f) is the wattle
flower’s optimally filtered true effective diffusion coefficient,
$D_{\textrm{eff, Atten}}^{\textrm{True}}(\textbf{r};\Delta)$. The greyscale
bars in subfigures (a) and (f) are $\times 10^{-5}$ $\mu$m, (b) is $\times
10^{-7}$ $\mu$m-1, and (c) and (d) are $\times 10^{-6}$.  Figure 4: Comparison
of wattle-flower multimodal signals extracted using the MIST approaches for
rotationally-isotropic diffuse scatter: (a)–(c) are reconstructed effective
diffusion coefficients, and (d)–(f) are phase-shifts (as the wattle flower is
a pure phase-object). (a) and (d) use Pavlov et al.’s[39] approach, which
neglects X-ray attenuation, and (b) and (e) use Alloo et al.’s[41] approach,
which considers it; both approaches assume the effective diffusion coefficient
to be slowly varying. (c) and (f) are calculated using the present approach,
which considers weak X-ray attenuation and does not approximate the diffusion
coefficient as slowly varying. The greyscale bars in (a)–(c) are $\times
10^{-5}$ $\mu$m and (d)–(f) are $\times 10^{0}$ rads.

Figure 3 shows the computed solutions,
$\nabla_{\perp}^{2}\left[\frac{1}{k}\phi_{\textrm{ob}}(\textbf{r})-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\right]$, $D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)$, $D_{\textrm{eff, Phase}}^{x}(\textbf{r};\Delta)$,
and $D_{\textrm{eff, Phase}}^{y}(\textbf{r};\Delta)$, and the multimodal
signals, $I_{\textrm{ob}}(\textbf{r})$ and $D_{\textrm{eff,
Atten}}^{\textrm{True}}(\textbf{r};\Delta)$, for the WF sample, where the
window and level of each image were set to optimise the respective greyscale
range. Although this sample has a weak DF signal, Figs. 3b–3d reveal the
first- and second-order derivatives of the diffuse scattering signal that had
been neglected in prior MIST approaches. It is obvious that the gradient term
is stronger at material interfaces on a global and local scale, that is the WF
leaf edges and the filaments that make up each flower, respectively. This
follows our initial prediction that the assumption of the DF signal being
slowly varying at material interfaces was insufficient; this is furthermore
supported in the case of the RC sample which is discussed later. Here, we
emphasise that the WF’s reconstructed $I_{\textrm{ob}}(\textbf{r})$ is close
to unity, confirming that it is a weakly-attenuating object.

Figure 4 compares the present MIST approach to those in the current MIST
research literature[39, 41]; the reconstructed DF signals are shown in Figs.
4a–4c and the phase-shifts in Figs. 4d–4f. Figures 4a and 4b are qualitatively
equivalent with regard to the resolvability of the WF’s leaves and filaments.
Quantitatively the reconstructed DF signal in Fig. 4b is larger than that in
Fig. 4a, owing to the former being a phase-object approximation and the latter
considering X-ray attenuation, and hence there is a higher reconstructed DF
signal in regions that attenuated the X-ray beam more; this agrees with what
was found in Alloo et al. [41]. The DF signal computed using the present
approach, Fig. 4c, initially appears to have a weaker reconstructed signal
within the WF’s leaves and filaments. However, after inspecting the
reconstructed phase-shifts using all MIST approaches, Figs. 4d–4f, it becomes
apparent the DF signals reconstructed using the alternative approaches have
predominant blurring at the WF’s flowers. Moreover, the WF’s DF signal from
the filaments that make up each flower cannot be well-described by assuming
the DF is slowly varying, as assumed in Pavlov et al.[39] and Alloo et
al.[41]. For this sample, the previous MIST formalisms are unable to detect
the rapidly varying features inside the WF filaments and its edges. This is
apparent in the reconstructions of the phase-shifts and effective diffusion
signal. This weakly-attenuating sample demonstrates that our new approach
gives reconstructed multimodal signals that are qualitatively superior
compared to the previously-published MIST approaches. Note that a quantitative
comparison of the reconstructed signals’ image quality is provided in a
following section.

Figure 5: Solutions of the system of linear equations, equations (Theoretical
derivation of updated MIST approach)–(Theoretical derivation of updated MIST
approach), and the red currants’s reconstructed multimodal signals; (a) is the
effective diffusion coefficient, $D_{\textrm{eff, Phase}}(\textbf{r};\Delta)$,
(b) is
$\nabla_{\perp}^{2}\left[\frac{1}{k}\phi_{\textrm{ob}}(\textbf{r})-D_{\textrm{eff,
Phase}}(\textbf{r};\Delta)\right]_{\textrm{Recon}}$, (c) and (d) are the two
spatial derivatives of the effective diffusion coefficient, $D_{\textrm{eff,
Phase}}^{y}(\textbf{r};\Delta)$ and $D_{\textrm{eff,
Phase}}^{x}(\textbf{r};\Delta)$, respectively, (e) is the reconstructed
attenuation term, $I_{\textrm{ob}}(\textbf{r})$, and (f) is the red currants’s
optimally filtered true effective diffusion coefficient, $D_{\textrm{eff,
Atten}}^{\textrm{True}}(\textbf{r};\Delta)$. The greyscale bars in subfigures
(a) and (f) are $\times 10^{-5}$ $\mu$m, (b) is $\times 10^{-7}$ $\mu$m-1, and
(c) and (d) are $\times 10^{-6}$.

### Attenuating sample: Red currant (RC)

We now turn to the non-negligibly-attenuating RC sample, which tests the
breadth of applicability of the proposed approach. Using an identical
methodology to that described for the WF, the over-determined system of linear
equations for seven sets of SB-PCXI intensity data was solved using Tikhonov-
regularised QR decomposition, from which the RC’s true phase-object DF
approximation, $D_{\textrm{eff, Phase}}^{\textrm{True}}(\textbf{r};\Delta)$,
attenuation term, $I_{\textrm{ob}}(\textbf{r})$, and true attenuating-object
DF, $D_{\textrm{eff, Atten}}^{\textrm{True}}(\textbf{r};\Delta)$, signals were
calculated. The cut-off parameter was chosen to be $\rho=$ 21
$\mu\textrm{m}{{}^{2}}$ by appropriately optimising the NIQE and SNR, as
described earlier.

Figure 5 shows all of the relevant solutions and multimodal signals for the RC
sample using the above-described method; Figs. 5a–5d show the reconstructed
solutions of the system of linear equations, and Figs. 5e and 5f display the
RC’s attenuation term and true attenuating-object approximation of its
effective diffusion coefficient, respectively. Similar to the case of the WF,
the gradient term is stronger at the periphery of the entire RC sample and at
the internal fibres. Figure 5e shows the RC’s attenuation term, where the
microstructure is indistinguishable due to its attenuation-based
contrast/signal being small. Such microstructure, which is unresolved in the
attenuation term, induces measurable SAXS, and hence, can be readily resolved
in the reconstructed DF signal shown in Fig. 5f. It is important to make note
of the RC’s reconstructed $I_{\textrm{ob}}(\textbf{r})$, particularly how far
from unity it is. This indicates that the attenuating-object approximation for
the effective diffusion coefficient should be used to reconstruct a true DF
signal.

Figure 6 gives a direct comparison of the RC’s DF and attenuation
reconstructions, using (a) the most recently published rotationally-isotropic-
scatter “DF Assumed Slowly Varying” MIST approach[41] (Figs. 6a and 6b), and
(b) the current “Rapidly Varying DF Behaviour” MIST approach (Figs. 6c and
6d). Here, we only compare the present approach to the rotationally-isotropic
MIST approach[41] that considers X-ray attenuation, as that published by
Pavlov et al.[39] would erroneously reconstruct the RC’s DF signal since the
sample significantly attenuates the X-ray beam. Figures 6a and 6c show a local
reconstruction of the RC’s $D_{\textrm{eff,
Atten}}^{\textrm{True}}(\textbf{r};\Delta)$ and Figs. 6b and 6d are the
$I_{\textrm{ob}}(\textbf{r})$ reconstructions; magnified regions of each
extracted signal are also shown. By comparing the reconstructed DF signals, it
is evident that the new approach increases spatial resolution, decreases
noise, and also provides a greater subject-contrast between unresolved
microstructure. Moreover, the small fibrous network surrounding the RC pip,
otherwise known as the pericarp, is more resolvable in Fig. 6c than Fig. 6a. A
line profile across a single RC fibre, indicated by the blue asterisk in Fig.
6c, is shown in Fig. 6e. This fine feature is resolvable in the DF image
reconstructed using the current approach (denoted by the green trace), but not
when the DF is assumed to be slowly varying (denoted by the red trace). When
the DF is assumed to be slowly varying, there are supposedly Fresnel fringes,
or residual PC, at the boundaries of the traced pericarp fibre in the
recovered DF signal; this is the high-low intensity region at approximately 30
$\mu m$ and 35 $\mu m$ in Fig. 6e. When this assumption is relaxed, the
apparent DF signal induced by the strong phase effects is reconstructed
appropriately, such that the small fibre is resolvable. This line profile
indicates the evident spatial resolution difference between the two images,
which is further supported by the azimuthally averaged power spectra in Fig.
6f. These power spectra were calculated in a global region in both DF images
which encased the entire RC but not the blu tack that supports the sample (the
larger yellow box in Fig. 5f). From these power spectra, it is evident that
there is a decrease in noise, shown by the reduction in high spatial frequency
components (lower noise floor), and an increase in spatial resolution, shown
by the higher spatial-frequency position of the power-spectrum knees (denoted
by dashed vertical lines), when the DF is considered to be rapidly varying.
Although the visibility of the RC’s pip and the pericarp is low in both
reconstructions of the attenuation term, Figs. 6b and 6d, the subject-contrast
appears higher when using the previous MIST approach (when the images are
shown on the same greyscale range). However, when the greyscale range is
optimised independently for each attenuation term reconstruction, as shown in
the magnified regions for each reconstruction, the images look almost
identical.

Figure 6: Comparison of the red currant’s multimodal signals extracted using
the approach published in Alloo et al.[41], (a)–(b), and that presented within
this study, (c)–(d). Figures (a)–(b) assume the red currant’s effective
diffusion coefficient is slowly varying and Figs. (c)–(d) have no assumptions
regarding the effective diffusion coefficient. The effective diffusion
coefficient reconstructions are shown in (a) and (c), and the RC’s attenuation
term is shown in (b) and (d). The blue-asterisk line profile in (c) is shown
in (e), and (f) is the azimuthally averaged power spectra. The vertical dashed
lines in (f) denote the approximate knee of the respective power spectra. The
greyscale bars in (a) and (c) are $\times 10^{-5}$ $\mu$m.  Figure 7: Image
quality metrics, a) reciprocal of the NIQE and b) SNR, of the reconstructed
effective diffusion coefficients using Alloo et al.’s[41] approach (denoted by
the red “DF Assumed Slowly Varying” label), and that presented within this
study (denoted by the green “Rapidly Varying DF Behaviour” label). Subscript
“Global” and “Local” denote the region of interest used for calculations;
these are indicated in Fig. 3f and Fig. 5f for the wattle flower (WF) and red
currant (RC), respectively.

## Image quality of reconstructed effective diffusion coefficients

To quantitatively compare the WF’s and RC’s reconstructed DF signals using the
present approach and those calculated using the DF slowly-varying attenuating-
object approach presented in Alloo et al.[41], we investigated the SNR and
NIQE[67] for distinct regions. Two regions in each of the reconstructed DF
signals were investigated for each sample; the so-called Global region which
contained the entire sample, and a Local region which encompassed only
critical structural features of the sample. The Global and Local regions,
denoted by the large and small yellow boxes, respectively, are indicated in
Fig. 3f and Fig. 5f, for the WF and RC, respectively. The NIQE and SNR were
calculated for both of the reconstructed DF signals in each of the described
regions, for each sample. Figure 7 provides a summary of the image quality
metrics using the two MIST approaches. The red markers denote the DF signal
calculated using the previous approach[41] which assumed the DF to be slowly
varying, but considers X-ray attenuation, and the green markers denote the new
approach which considers rapidly-varying DF behaviour but only weak X-ray
attenuation. It is shown from the SNR that the new approach reconstructs
superior-quality DF images across all regions in both samples. The reciprocal
NIQEs show similar behaviour for the DF signals, with only one datapoint as an
exception, where the local reconstruction of the red current is perceived as
“better” when the DF is assumed to be slowly varying. This scoring can be
explained by the previously-mentioned point describing how the NIQE interprets
the residual PC as a sharper edge, due to its increased contrast, when in
reality it is incorrectly describing the DF signal. It may be argued that the
overall increase in image quality shown in the remaining datapoints is only
due to the filtering performed on the new DF reconstruction. However, multiple
different filters (e.g., Gaussian and median) were applied to the “DF slowly-
varying reconstruction”s, for which the NIQE and SNR were calculated, and
there was no such filtering that gave a comparable image quality to that
obtained using the present approach.

The described results for the WF and RC samples reveal that the present MIST
approach reconstructs superior DF signals for weakly- and non-negligibly-
attenuating samples, compared to our previously-published MIST approaches. We
have also demonstrated that (a) for a weakly-attenuating object the proposed
approach gives a qualitatively better PC reconstruction, and (b) for an
attenuating sample both approaches give PC reconstructions with similar image
quality. These conclusions were made based on the image quality measures of
the reconstructions, rather than the quantitativeness of the reconstructions.
Unsurprisingly, the quantitative difference between the reconstructed phase-
shift, $\phi_{\textrm{ob}}(\textbf{r})$, from the two approaches, is larger
for an attenuating object than for a weakly-attenuating object, which is
reconstructed equivalently. This is exactly what is expected, based on the
underlying assumptions of both theoretical formalisms. That is, the “DF
slowly-varying” approach[41] considers X-ray attenuation in the initial
Fokker-Planck description, whereas the present approach neglects attenuation
initially, before extending to the case of a weakly-attenuating material using
the projection approximation. It, therefore, follows that the previous
approach more accurately reconstructs the PC signal for an attenuating object,
compared to the new approach. The two MIST approaches are complementary
methods derived using their own distinct assumptions, therefore it is expected
that their respective applicability will be sample dependent. There are
regimes in which one approach is more suitable than the other. The present
approach is optimal for weakly-attenuating objects (e.g., WF), however, can
also be successfully applied to attenuating objects (e.g., RC).

## Concluding remarks

Within this paper, we developed and implemented a rapid deterministic approach
that can reconstruct high-resolution multimodal signals of samples using SB-
PCXI. The present MIST approach is not restricted to samples that have a
slowly-varying DF signal, in contrast to other MIST approaches[39, 41, 42],
and instead can be used to model rapidly-varying DF behaviour. We applied the
approach to two samples that differ in X-ray attenuation characteristics and
compared these signals to those reconstructed using two earlier variants of
MIST[39, 41]. Using the new approach, the SNR, spatial resolution, and
perceived image quality—in the majority of local and global regions-of-
interest, across the DF reconstructions for both samples—were higher.

Multimodal X-ray imaging has already proven useful in numerous applications,
and this work provides theoretical development towards reconstructing the best
possible images using an SB-PCXI technique. It furthermore assists the
translation of SB-PCXI into a user-friendly low-dose technique, as the
proposed approach requires just 4 sets of SB-PCXI data. It is computationally
efficient, requiring just 3 minutes to calculate multimodal signals for a
2100$\times$2500 pixel image using a laptop computer with an 11th Gen Intel(R)
Core(TM) i7-1165G7 2.80 GHz processor and 64 GB RAM. Moreover, the
experimental setup is simple, and the SB-PCXI technique itself has low
coherence requirements.

We anticipate that the present approach can be extended to the case of a
highly-attenuating object in which the DF is spatially rapidly varying, and
also to the case of rotationally-anisotropic position-dependent SAXS[42]. The
directional DF[43, 44] approach using MIST[42] assumes the DF to be slowly
varying. By applying the present approach to solving the directional DF
inverse problem, it may be possible to optimise, by means of increased spatial
resolution and SNR, the two-dimensional reconstructions, thereby furthering
the future goal of MIST tensor-tomography. The present MIST approach could
also be used to provide a rapidly-computed deterministic initial guess for
alternative speckle-tracking approaches that solve the inverse SB-PCXI problem
iteratively[25, 34, 30, 31]. Such iterative techniques—which have a broader
domain of applicability because they make fewer assumptions than is the case
for our work—are computationally expensive, as the multimodal signals are
reconstructed with no definite initial guess. These approaches might be more
computationally efficient if they had some “good” initial guess for the
sample’s multimodal signals, such as that provided by the present MIST
approach. The speed of convergence of subsequent iterative refinement might
significantly increase, also making these approaches more appealing for
broader adoption.

## Acknowledgements

The authors thank Sebastien Berujon, Eric Ziegler, and Emmanuel Brun for
collecting and sharing the SB-PCXI data of the red currant sample collected at
ESRF, published originally in Berujon et al.[30]. The authors are grateful for
the help provided by the beamline scientists, Chris Hall, Daniel Hausermann,
Anton Maksimenko, and Matthew Cameron, at the Imaging and Medical beamline at
the Australian Synchrotron, part of ANSTO, where the images in Figs. 3 and 4
were captured under proposal 18648. K. S. Morgan acknowledges support from the
Australian Research Council (FT18010037). This research was supported by an
AINSE Ltd. Postgraduate Research Award (PGRA) and the New Zealand Synchrotron
Group Limited’s capability fund grant. We acknowledge the University of
Canterbury for awarding a doctoral scholarship to S. J. Alloo.

## Competing interests

The authors declare no competing interests.

## Data availability

The Python3 script, with appropriate test data, is available in the open-
access repository on GitHub[72]. Further experimental data are available upon
reasonable request, please contact the corresponding author, S. J. Alloo.

## References

  * [1] David M Paganin.  Coherent X-ray Optics.  Oxford University Press, Oxford, 2006.
  * [2] Atsushi Momose, Shinya Kawamoto, Ichiro Koyama, Yoshitaka Hamaishi, Kengo Takai, and Yoshio Suzuki.  Demonstration of X-ray Talbot interferometry.  Jpn. J. Appl. Phys., 42(7B):L866–L868, 2003.
  * [3] Eckhart Förster, K Goetz, and P Zaumseil.  Double crystal diffractometry for the characterization of targets for laser fusion experiments.  Kristall und Technik, 15(8):937–945, 1980.
  * [4] Han Wen, Eric E Bennett, Monica M Hegedus, and Stefanie C Carroll.  Spatial harmonic imaging of x-ray scattering—initial results.  IEEE Trans. Med. Imaging, 27(8):997–1002, 2008.
  * [5] Kaye S Morgan, David M Paganin, and Karen KW Siu.  Quantitative single-exposure x-ray phase contrast imaging using a single attenuation grid.  Opt. Express, 19(20):19781–19789, 2011.
  * [6] A Olivo, Fulvia Arfelli, Giovanni Cantatore, Renata Longo, RH Menk, S Pani, M Prest, P Poropat, Luigi Rigon, G Tromba, and E. Vallazza.  An innovative digital imaging set-up allowing a low-dose approach to phase contrast applications in the medical field.  Med. Phys., 28(8):1610–1619, 2001.
  * [7] A Snigirev, I Snigireva, V Kohn, S Kuznetsov, and I Schelokov.  On the possibilities of x-ray phase contrast microimaging by coherent high-energy synchrotron radiation.  Rev. Sci. Instrum., 66(12):5486–5492, 1995.
  * [8] Peter Cloetens, Raymond Barrett, José Baruchel, Jean-Pierre Guigay, and Michel Schlenker.  Phase objects in synchrotron radiation hard x-ray imaging.  Journal of Physics D: Applied Physics, 29(1):133–146, 1996.
  * [9] Hye-Won Jung, Ivan Lee, Sang-Heon Lee, Kaye S Morgan, David Parsons, and Martin Donnelley.  Mucociliary transit assessment using automatic tracking in phase contrast x-ray images of live mouse nasal airways.  J. Med. Biol. Eng., 42:545–554, 2022.
  * [10] S Drevet, B Favier, B Lardy, G Gavazzi, and E Brun.  New imaging tools for mouse models of osteoarthritis.  GeroScience, 44(2):639–650, 2022.
  * [11] Gabriela Sena, Gabriel Fidalgo, Katrine Paiva, Renan Barcelos, Liebert Parreiras Nogueira, Marcos Vinícius Colaço, Marcelo Salabert Gonzalez, Patricia Azambuja, Gustavo Colaço, Helio Ricardo da Silva, Anderson Alvarenga de Moura Meneses, and Regina C Barroso.  Synchrotron X-ray biosample imaging: opportunities and challenges.  Biophys. Rev., 14:625–633, 2022.
  * [12] Lorenzo Massimi, Charlotte K Hagen, Marco Endrizzi, Peter RT Munro, Glafkos Havariyoun, PM Sam Hawker, Bennie Smit, Alberto Astolfo, Oliver J Larkin, Richard M Waltham, Z Shah, Stephen W. Duffy, Rachel L. Nelan, Anthony Peel, Tamara Suaris, J. Louise Jones, Ian G. Haig, and Alessandro Olivo.  Laboratory-based x-ray phase contrast CT technology for clinical intra-operative specimen imaging.  Medical Imaging 2019: Physics of Medical Imaging, 10948:109481R, 2019.
  * [13] Navnath S Indore, Chithra Karunakaran, and Digvir S Jayas.  Synchrotron tomography applications in agriculture and food sciences research: a review.  Plant Methods, 18:101–127, 2022.
  * [14] Andrew McNeil, Robert S Bradley, Philip J Withers, and David Penney.  Imaging fossilised spiders in amber using lab-based phase contrast X-ray tomography.  Developments in X-ray Tomography VII, 7804:440–446, 2010.
  * [15] Gregory D Edgecombe, Varpu Vahtera, Stuart R Stock, Aki Kallonen, Xianghui Xiao, Alexander Rack, and Gonzalo Giribet.  A scolopocryptopid centipede (Chilopoda: Scolopendromorpha) from Mexican amber: synchrotron microtomography and phylogenetic placement using a combined morphological and molecular data set.  Zool. J. Linn. Soc., 166(4):768–786, 2012.
  * [16] Melissa A Denecke, Wout de Nolf, Alexander Rack, Remi Tucoulou, Tonya Vitova, Gerald Falkenberg, Sousan Abolhassani, Peter Cloetens, and Bernhard Kienzler.  Speciation of actinides in granite subjected to tracer studies.  In Stepan N. Kalmykov and Melissa A. Denecke, editors, Actinide Nanoparticle Research, pages 413–435. Springer, Berlin, 2011.
  * [17] S C Mayo, P Miller, S W Wilkins, D Gao, and T E Gureyev.  Laboratory-based x-ray micro-tomography with submicron resolution.  Developments in X-ray Tomography V, 6318:439–446, 2006.
  * [18] Franz Pfeiffer, Martin Bech, Oliver Bunk, Philipp Kraft, Eric F Eikenberry, Ch Brönnimann, Christian Grünzweig, and Christian David.  Hard-X-ray dark-field imaging using a grating interferometer.  Nat. Mater., 7(2):134–137, 2008.
  * [19] E Pagot, P Cloetens, S Fiedler, A Bravin, P Coan, J Baruchel, J Härtwig, and W Thomlinson.  A method to extract quantitative information in analyzer-based x-ray phase contrast imaging.  Appl. Phys. Lett., 82(20):3421–3423, 2003.
  * [20] Miles N Wernick, Oliver Wirjadi, Dean Chapman, Zhong Zhong, Nikolas P Galatsanos, Yongyi Yang, Jovan G Brankov, Oral Oltulu, Mark A Anastasio, and Carol Muehleman.  Multiple-image radiography.  Phys. Med. Biol., 48(23):3875–3896, 2003.
  * [21] Ying Ying How and Kaye S Morgan.  Quantifying the x-ray dark-field signal in single-grid imaging.  Opt. Express, 30(7):10899–10918, 2022.
  * [22] M Endrizzi, BIS Murat, P Fromme, and A Olivo.  Edge-illumination X-ray dark-field imaging for visualising defects in composite structures.  Compos. Struct., 134:895–899, 2015.
  * [23] T E Gureyev, D M Paganin, B Arhatari, S T Taba, S Lewis, P C Brennan, and H M Quiney.  Dark-field signal extraction in propagation-based phase-contrast imaging.  Phys. Med. Biol., 65(21):215029, 2020.
  * [24] Thomas A Leatham, David M Paganin, and Kaye S Morgan.  X-ray dark-field and phase retrieval without optics, via the Fokker-Planck equation.  arXiv:2112.10999, 2021.
  * [25] Sébastien Bérujon, Eric Ziegler, Roberto Cerbino, and Luca Peverini.  Two-dimensional x-ray beam phase sensing.  Phys. Rev. Lett., 108(15):158102, 2012.
  * [26] Kaye S Morgan, David M Paganin, and Karen KW Siu.  X-ray phase imaging with a paper analyzer.  Appl. Phys. Lett., 100(12):124102, 2012.
  * [27] Joseph W Goodman.  Speckle Phenomena in Optics: Theory and Applications.  SPIE Press Bellingham, Washington, 2020.
  * [28] Marie-Christine Zdora.  State of the art of X-ray speckle-based phase-contrast and dark-field imaging.  J. Imaging, 4(5):60, 2018.
  * [29] P. C. Sabatier.  Past and future of inverse problems.  J. Math. Phys., 41(6):4082–4124, 2000.
  * [30] Sebastien Berujon and Eric Ziegler.  X-ray multimodal tomography using speckle-vector tracking.  Phys. Rev. Appl., 5(4):044014, 2016.
  * [31] Sebastien Berujon and Eric Ziegler.  Near-field speckle-scanning-based x-ray tomography.  Phys. Rev. A, 95(6):063822, 2017.
  * [32] Hongchang Wang, Yogesh Kashyap, and Kawal Sawhney.  From synchrotron radiation to lab source: advanced speckle-based X-ray imaging using abrasive paper.  Sci. Rep., 6(1):20476, 2016.
  * [33] Hongchang Wang, Yogesh Kashyap, Biao Cai, and Kawal Sawhney.  High energy X-ray phase and dark-field imaging using a random absorption mask.  Sci. Rep., 6(1):30581, 2016.
  * [34] Marie-Christine Zdora, Pierre Thibault, Tunhe Zhou, Frieder J Koch, Jenny Romell, Simone Sala, Arndt Last, Christoph Rau, and Irene Zanette.  X-ray phase-contrast imaging and metrology through unified modulated pattern analysis.  Phys. Rev. Lett., 118(20):203903, 2017.
  * [35] David M Paganin, Hélène Labriet, Emmanuel Brun, and Sebastien Berujon.  Single-image geometric-flow x-ray speckle tracking.  Phys. Rev. A, 98(5):053813, 2018.
  * [36] David M Paganin and Kaye S Morgan.  X-ray Fokker–Planck equation for paraxial imaging.  Sci. Rep., 9:17537, 2019.
  * [37] Kaye S Morgan and David M Paganin.  Applying the Fokker–Planck equation to grating-based x-ray phase and dark-field imaging.  Sci. Rep., 9:17465, 2019.
  * [38] M. R. Teague.  Deterministic phase retrieval: a Green’s function solution.  J. Opt. Soc. Am., 73(11):1434–1441, 1983.
  * [39] Konstantin M Pavlov, David M Paganin, Heyang Thomas Li, Sebastien Berujon, Hélène Rougé-Labriet, and Emmanuel Brun.  X-ray multi-modal intrinsic-speckle-tracking.  J.Opt., 22(12):125604, 2020.
  * [40] Hannes Risken.  The Fokker-Planck equation: Methods of Solution and Applications.  Springer, Berlin, 2 edition, 1989.
  * [41] Samantha J Alloo, David M Paganin, Kaye S Morgan, Marcus J Kitchen, Andrew W Stevenson, Sheridan C Mayo, Heyang T Li, Ben M Kennedy, Anton Maksimenko, Joshua C Bowden, and Konstantin M Pavlov.  Dark-field tomography of an attenuating object using intrinsic x-ray speckle tracking.  J. Med. Imaging, 9(3):031502, 2022.
  * [42] Konstantin M Pavlov, David M Paganin, Kaye S Morgan, Heyang Thomas Li, Sebastien Berujon, Laurène Quénot, and Emmanuel Brun.  Directional dark-field implicit x-ray speckle tracking using an anisotropic-diffusion Fokker-Planck equation.  Phys. Rev. A, 104(5):053505, 2021.
  * [43] Torben H Jensen, Martin Bech, Oliver Bunk, Tilman Donath, Christian David, Robert Feidenhans’l, and Franz Pfeiffer.  Directional x-ray dark-field imaging.  Phys. Med. Biol., 55(12):3317–3323, 2010.
  * [44] Torben Haugaard Jensen, Martin Bech, Irene Zanette, Timm Weitkamp, Christian David, Hans Deyhle, Simon Rutishauser, Elena Reznikova, Jürgen Mohr, Robert Feidenhans’l, and Franz Pfeiffer.  Directional x-ray dark-field imaging of strongly ordered systems.  Phys. Rev. B, 82:214103, Dec 2010.
  * [45] Avinash C Kak and Malcolm Slaney.  Principles of Computerized Tomographic Imaging.  IEEE Press, New York, 1988.
  * [46] Erin A Miller, Timothy A White, Benjamin S McDonald, and Allen Seifert.  Phase contrast x-ray imaging signatures for security applications.  IEEE Trans. Nucl. Sci., 60(1):416–422, 2013.
  * [47] Hyunwoo Lim, Jonghyeok Lee, Soohyun Lee, Hyosung Cho, Hunwoo Lee, and Duhee Jeon.  Low-density foreign body detection in food products using single-shot grid-based dark-field X-ray imaging.  J. Food Eng., 335:111189, 2022.
  * [48] Masami Ando, Naoki Sunaguchi, Daisuke Shimao, Adam Pan, Tetsuya Yuasa, Kensaku Mori, Yoshifumi Suzuki, Ge Jin, Jong-Ki Kim, Jae-Hong Lim, Seung-Jun Seo, Shu Ichihara, Norihiko Ohura, and Rajiv Gupta.  Dark-field imaging: recent developments and potential clinical applications.  Phys. Med., 32(12):1801–1812, 2016.
  * [49] Andre Yaroshenko, Katharina Hellbach, Ali Önder Yildirim, Thomas M Conlon, Isis Enlil Fernandez, Martin Bech, Astrid Velroyen, Felix G Meinel, Sigrid Auweter, Maximilian Reiser, Oliver Eickelber, and Franz Pfeiffer.  Improved in vivo assessment of pulmonary fibrosis in mice using x-ray dark-field radiography.  Sci. Rep., 5(1):17492, 2015.
  * [50] Katharina Hellbach, Andre Yaroshenko, Konstantin Willer, Thomas M Conlon, Margarita B Braunagel, Sigrid Auweter, Ali Ö Yildirim, Oliver Eickelberg, Franz Pfeiffer, Maximilian F Reiser, and Felix G Meinel.  X-ray dark-field radiography facilitates the diagnosis of pulmonary fibrosis in a mouse model.  Sci. Rep., 7(1):340, 2017.
  * [51] Simone Schleede, Felix G Meinel, Martin Bech, Julia Herzen, Klaus Achterhold, Guillaume Potdevin, Andreas Malecki, Silvia Adam-Neumair, Sven F Thieme, Fabian Bamberg, Konstantin Nikolaou, Alexander Bohla, Ali Ö Yildirim, Roderick Loewen, Martin Gifford, Ronald Ruth, Oliver Eickelberg, Maximilian F Reiser, and Franz Pfeiffer.  Emphysema diagnosis using X-ray dark-field imaging at a laser-driven compact synchrotron light source.  Proceedings of the National Academy of Sciences, 109(44):17880–17885, 2012.
  * [52] Felix G Meinel, Felix Schwab, Simone Schleede, Martin Bech, Julia Herzen, Klaus Achterhold, Sigrid Auweter, Fabian Bamberg, Ali Ö Yildirim, Alexander Bohla, Oliver Eickelber, Maximilian F Reiser, Franz Pfeiffer, and Konstantin Nikolaou.  Diagnosing and mapping pulmonary emphysema on X-ray projection images: incremental value of grating-based X-ray dark-field imaging.  PloS one, 8(3):e59526, 2013.
  * [53] Matthew R Arnison, Kieran G Larkin, Colin JR Sheppard, Nicholas I Smith, and Carol J Cogswell.  Linear phase imaging using differential interference contrast microscopy.  J. Microsc., 214(1):7–12, 2004.
  * [54] C Kottler, C David, F Pfeiffer, and O Bunk.  A two-directional approach for grating based differential phase contrast imaging using hard x-rays.  J. Microsc., 15(3):1175–1181, 2007.
  * [55] MD De Jonge, B Hornberger, C Holzner, D Legnini, D Paterson, I McNulty, C Jacobsen, and S Vogt.  Quantitative phase imaging with a scanning transmission x-ray microscope.  Physical Review Letters, 100(16):163902, 2008.
  * [56] D. Paganin, A. Barty, P. J. McMahon, and K. A. Nugent.  Quantitative phase-amplitude microscopy. III. The effects of noise.  Journal of Microscopy, 214(1):51–61, 2004.
  * [57] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery.  Numerical Recipes in FORTRAN: The Art of Scientific Computing (2nd edn).  Cambridge University Press, 1996.
  * [58] David M Paganin, Sheridan C Mayo, Tim E Gureyev, Peter R Miller, and Steve W Wilkins.  Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object.  J. Microsc., 206(1):33–40, 2002.
  * [59] Timur Eugenievich Gureyev, Andrew W Stevenson, David Maurice Paganin, T Weitkamp, A Snigirev, I Snigireva, and SW Wilkins.  Quantitative analysis of two-component samples using in-line hard X-ray images.  Journal of Synchrotron Radiation, 9(3):148–153, 2002.
  * [60] Mario Alejandro Beltran, DM Paganin, Kentaro Uesugi, and Marcus John Kitchen.  2D and 3D x-ray phase retrieval of multi-material objects using a single defocus distance.  Optics Express, 18(7):6423–6436, 2010.
  * [61] Timur E Gureyev, Sara Mohammadi, Yakov Nesterets, Christian Dullin, and Giuliana Tromba.  Accuracy and precision of reconstruction of complex refractive index in near-field single-distance propagation-based phase-contrast tomography.  J. Appl. Phys., 114(14):144906, 2013.
  * [62] Ehsan Fathi and Babak Maleki Shoja.  Deep neural networks for natural language processing.  In Venkat N. Gudivada and C.R. Rao, editors, Computational Analysis and Understanding of Natural Languages: Principles, Methods and Applications, volume 38 of Handbook of Statistics, pages 229–316. Elsevier, 2018.
  * [63] Almerico Murli, Luisa D’Amore, and Valentina De Simone.  The Wiener filter and regularization methods for image restoration problems.  Proceedings 10th International Conference on Image Analysis and Processing, pages 394–399, 1999.
  * [64] Andrey N Tikhonov and Vasiliy Y Arsenin.  Solutions of Ill-posed Problems.  Winston, Washington DC, 1977.
  * [65] Jianping Zhu.  QR factorization for the regularized least squares problem on hypercubes.  Parallel Comput., 19(8):939–948, 1993.
  * [66] Yonggi Park, Lothar Reichel, Giuseppe Rodriguez, and Xuebo Yu.  Parameter determination for Tikhonov regularization problems in general form.  J. Comput. Appl. Math., 343:12–25, 2018.
  * [67] Anish Mittal, Rajiv Soundararajan, and Alan C Bovik.  Making a “completely blind” image quality analyzer.  IEEE Signal Process Lett., 20(3):209–212, 2012.
  * [68] Norman B Nill.  Scene power spectra: the moment as an image quality merit factor.  Appl. Opt., 15(11):2846–2854, 1976.
  * [69] Andrew W Stevenson, Jeffrey C Crosbie, Christopher J Hall, Daniel Häusermann, Jayde Livingstone, and Jessica E Lye.  Quantitative characterization of the X-ray beam at the Australian Synchrotron Imaging and Medical Beamline (IMBL).  J. Synchrotron Radiat., 24(1):110–141, 2017.
  * [70] J-C Labiche, J Segura-Puchades, D Van Brussel, and J-P Moy.  FRELON camera: fast readout low noise.  ESRF Newsletter, 25:41–43, 1996.
  * [71] P-A Douissard, A Cecilia, X Rochet, X Chapel, T Martin, Th van de Kamp, L Helfen, T Baumbach, L Luquot, X Xiao, J Meinhardt, and A Rack.  A versatile indirect detector design for hard X-ray microimaging.  J. Instrum., 7(09):P09016, 2012.
  * [72] Samantha J Alloo.  Multimodal Intrinsic Speckle Tracking: Python Repository, 2022.  https://github.com/samanthaalloo/MIST_Alloo23.
  * [73] TS Imaging, 2022.  http://ts-imaging.science.unimelb.edu.au/Services/Simple/.