Applying the Fokker--Planck equation to grating-based x-ray phase and dark-field imaging

TL;DR

This paper applies the Fokker--Planck equation to model the drift and diffusion of structured illumination in grating-based x-ray phase and dark-field imaging, aiming to understand signal interactions and improve image reconstruction.

Contribution

It introduces a novel application of the Fokker--Planck equation to describe phase and dark-field signals in x-ray imaging, providing insights into signal cross-talk.

Findings

Fokker--Planck equation models drift and diffusion in structured illumination.

Enhanced understanding of cross-talk between attenuation, phase, and dark-field signals.

Potential basis for improved inverse problem solutions in x-ray imaging.

Abstract

X-ray imaging has conventionally relied upon attenuation to provide contrast. In recent years, two complementary modalities have been added; phase contrast and dark-field x-ray imaging, capturing weakly attenuating and sub-pixel sample structures respectively. These three modalities can be accessed using a crystal analyser, a grating interferometer or by looking at a directly-resolved grid, grating or speckle pattern. Grating and grid-based methods extract a differential phase signal by measuring how far a feature in the illumination has been shifted transversely due to the presence of a sample. The dark-field signal is extracted by measuring how the visibility of the structured illumination is decreased, typically due to the presence of sub-pixel structures in a sample. The strength of the dark-field signal may depend on the grating period, the pixel size and the set-up distances, and additional dark-field signal contributions may be seen as a result of strong phase effects or other factors. In this paper we show that the finite-difference form of the Fokker--Planck Equation can be applied to describe the drift (phase signal) and diffusion (dark-field signal) of the periodic or structured illumination used in phase contrast x-ray imaging with gratings, in order to better understand any cross-talk between attenuation, phase and dark-field x-ray signals. In future work, this mathematical description could be used as a basis for new approaches to the inverse problem of recovering both phase and dark-field information.

Equations

$$\frac{\partial }{\partial t}p(x,t)=-\frac{\partial }{\partial x}(D_1(x,t)p(x,t))+\frac{{\partial }^2}{\partial x^2}(D_2(x,t)p(x,t)),$$

$$\frac{\partial }{\partial z}I(x)=-\frac{1}{k}\frac{\partial }{\partial x}\left(I\frac{\partial \varphi (x)}{\partial x}\right)+\frac{{\partial }^2}{\partial x^2}(D(x,z)I(x)).$$

$$\frac{\partial }{\partial t}p(\text{x},t)=\sum _{r=1}^{\infty }\frac{(-1{)}^r}{r!}\frac{{\partial }^r}{\partial x^r}\left[D_r(\text{x},t)p(x,t)\right]=-\frac{\partial }{\partial x_i}\left[D_i(\text{x},t)p(\text{x},t)\right]+\frac{1}{2!}\frac{{\partial }^2}{\partial x_i\partial x_j}\left[D_{ij}(\text{x},t)p(\text{x},t)\right]-\cdots$$

$$\frac{\partial }{\partial z}I(\text{x},z)=-\frac{\partial }{\partial x}\left(\frac{1}{k}\frac{\partial \varphi }{\partial x}I(\text{x},z)\right)-\frac{\partial }{\partial y}\left(\frac{1}{k}\frac{\partial \varphi }{\partial y}I(\text{x},z)\right)+\frac{{\partial }^2}{\partial x^2}(D_{x}I(\text{x},z))+\frac{{\partial }^2}{\partial y^2}(D_{y}I(\text{x},z))+\frac{{\partial }^2}{\partial x\partial y}(D_{x,y}I(\text{x},z))$$

$$I=a\sin \left(\frac{x}{p}\right)+b.$$

$$I(x,z=\Delta )\approx I(x,z=0)+\Delta \left(-\frac{1}{k}\frac{\partial }{\partial x}\left(I\frac{\partial \varphi }{\partial x}\right)+\frac{{\partial }^2}{\partial x^2}(DI)\right).$$

$$I(x,z=\Delta )\approx \left(a\sin \left(\frac{x}{p}\right)+b\right)+\Delta \frac{{\partial }^2}{\partial x^2}\left(D\left(a\sin \left(\frac{x}{p}\right)+b\right)\right).$$

$$I(x,z=\Delta )=\left(1-\frac{D\Delta }{p^2}\right)a\sin \left(\frac{x}{p}\right)+b.$$

$$I(z=\Delta )={\int }_{-\infty }^{\infty }I(x^{\prime },z=0)\frac{e^{-\frac{(x-x^{\prime }{)}^2}{2\times D(x^{\prime })\Delta }}}{\sqrt{2\pi \times D(x^{\prime })\Delta }}d{x}^{\prime }.$$

$$I(z=\Delta )=a{e}^{-\frac{D\Delta }{2p^2}}\sin \left(\frac{x}{p}\right)+b.$$

$$I(x,z=\Delta )=a\sin \left(\frac{x}{p}\right)+b+\Delta \left[\frac{{\partial }^{2}D}{\partial x^2}\left(a\sin \left(\frac{x}{p}\right)+b\right)-\frac{D}{p^2}a\sin \left(\frac{x}{p}\right)+\frac{2}{p}\frac{\partial D}{\partial x}a\cos \left(\frac{x}{p}\right)\right].$$

$$I(x,z=\Delta )=a\sqrt{{\left(1+\Delta \frac{{\partial }^{2}D}{\partial x^2}-\frac{D\Delta }{p^2}\right)}^2+{\left(\frac{2\Delta }{p}\frac{\partial D}{\partial x}\right)}^2}\sin \left(\frac{x}{p}+\frac{2\Delta p\frac{\partial D}{\partial x}}{p^2+\Delta p^2\frac{{\partial }^{2}D}{\partial x^2}-D\Delta }\right)+b\left(1+\Delta \frac{{\partial }^{2}D}{\partial x^2}\right).$$

$$D(x)=((D_{right}-D_{left})/2)\text{Erf}[x]+((D_{right}+D_{left})/2)$$

$$I(x,z=\Delta )=a\sqrt{{\left(1-\frac{\Delta \left(\frac{1}{2}\text{Erf}(x)(D_{right}-D_{left})+\frac{D_{left}+D_{right}}{2}\right)}{p^2}-\frac{\Delta \left(2e^{-x^2}x(D_{right}-D_{left})\right)}{\sqrt{\pi }}\right)}^2+{\left(\frac{\Delta \left(e^{-x^2}(D_{right}-D_{left})\right)}{\sqrt{\pi }p}\right)}^2}$$

$$\times \text{ sin}\left(\frac{x}{p}-\frac{\Delta p\left(e^{-x^2}(D_{right}-D_{left})\right)}{\sqrt{\pi }\left(-\Delta \left(\frac{1}{2}\text{Erf}(x)(D_{right}-D_{left})+\frac{D_{left}+D_{right}}{2}\right)-\frac{\Delta p\left(2e^{-x^2}x(D_{right}-D_{left})\right)}{\sqrt{\pi }}+p^2\right)}\right)$$

$$+b\left(1-\frac{\Delta \left(2e^{-x^2}x(D_{right}-D_{left})\right)}{\sqrt{\pi }}\right).$$

$$I(z=\Delta )={\int }_{-\infty }^0\left(a\text{sin}\left(\frac{x^{\prime }}{p}\right)+b\right)\frac{e^{-\frac{(x-x^{\prime }{)}^2}{2D_{left}\Delta }}}{\sqrt{2\pi D_{left}\Delta }}d{x}^{\prime }+{\int }_0^{\infty }\left(a\text{sin}\left(\frac{x^{\prime }}{p}\right)+b\right)\frac{e^{-\frac{(x-x^{\prime }{)}^2}{2D_{right}\Delta }}}{\sqrt{2\pi D_{right}\Delta }}d{x}^{\prime }.$$

$$I(z=\Delta )=\frac{ia}{4}{\mathrm{e}}^{-\frac{2ipx+(D_{left}+D_{right})\Delta }{2p^2}}[-{\mathrm{e}}^{\frac{D_{left}\Delta }{2p^2}}\left(-1-\mathrm{Erf}\left(\frac{px-i{D}_{right}\Delta }{p\sqrt{2D_{right}\Delta }}\right)+{\mathrm{e}}^{\frac{2ix}{p}}\left(1+\mathrm{Erf}\left(\frac{px+i{D}_{right}\Delta }{p\sqrt{2D_{right}\Delta }}\right)\right)\right)$$

$$+{\mathrm{e}}^{\frac{D_{right}\Delta }{2p^2}}\left(1-\mathrm{Erf}\left(\frac{px-i{D}_{left}\Delta }{p\sqrt{2D_{left}\Delta }}\right)-{\mathrm{e}}^{\frac{2ix}{p}}\left(1-\mathrm{Erf}\left(\frac{px+i{D}_{left}\Delta }{p\sqrt{2D_{left}\Delta }}\right)\right)\right)]$$

$$+2b{\mathrm{e}}^{\frac{2ipx+(D_{left}+D_{right})\Delta }{2p^2}}\left(2-\mathrm{Erf}\left(\frac{x}{\sqrt{2D_{left}\Delta }}\right)+\mathrm{Erf}\left(\frac{x}{\sqrt{2D_{right}\Delta }}\right)\right).$$

$$\frac{\partial }{\partial z}I=\frac{{\partial }^2}{\partial x^2}(DI).$$

$$\frac{\partial }{\partial z}I=-\frac{D}{p^2}I.$$

$$I(x,z)=I(x,z=0)e^{-\frac{Dz}{p^2}},$$

$$I(x,z=\Delta )=a\sin \left(\frac{x}{p}\right)+b-\frac{\Delta }{k}\left[\left(a\sin \left(\frac{x}{p}\right)+b\right)\frac{{\partial }^2\varphi }{\partial x^2}+\frac{1}{p}a\cos \left(\frac{x}{p}\right)\frac{\partial \varphi }{\partial x}\right].$$

$$I(x,z=\Delta )=\left(1-\frac{\Delta }{k}\frac{{\partial }^2\varphi }{\partial x^2}\right)\left[a\sin \left(\frac{x}{p}-{\tan }^{-1}\left(\frac{\frac{\partial \varphi }{\partial x}}{p\left(\frac{k}{\Delta }-\frac{{\partial }^2\varphi }{\partial x^2}\right)}\right)\right)+b\right].$$

$$I(x,z=\Delta )=a\sin \left(\frac{x}{p}-{\tan }^{-1}\left(\frac{\Delta }{pk}\frac{\partial \varphi }{\partial x}\right)\right)+b.$$

$$I(x,z=\Delta )$$

$$I(x,z=\Delta )=\left(\frac{\Delta \left(2e^{-x^2}x({\varphi }_{right}-{\varphi }_{left})\right)}{\sqrt{\pi }k}+1\right)\left(a\sin \left(\left(\frac{x}{p}-{\tan }^{-1}\left(\frac{e^{-x^2}({\varphi }_{right}-{\varphi }_{left})}{p\sqrt{\pi }\left(\frac{k}{\Delta }+\frac{2e^{-x^2}x({\varphi }_{right}-{\varphi }_{left})}{\sqrt{\pi }}\right)}\right)\right)\right)+b\right).$$

$$I(x,z=\Delta )$$

$$I(x,z=\Delta )=MV\sin \left(\frac{2\pi (x-s)}{p}\right)+M.$$

$$\text{Stepping curve}(n)$$

$$=$$

$$I_{left}$$

Full text

1 Introduction

X-ray imaging has been widely adopted in medicine and security, providing non-invasive visualisation of internal structure at high spatial and temporal resolution. The conventional mechanism for x-ray contrast is attenuation, where dense and/or thick objects reduce the intensity of the illuminating x-ray wavefield. More recently, the capabilities of x-ray imaging have been extended in research laboratories and synchrotrons, capturing weakly-attenuating features like soft biological tissue [1] through measurements of the wavefield phase, and revealing sub-pixel structures like dentinal tubules [2] via a ‘dark-field’ signal that originates from unresolved scattering. Because an x-ray detector is typically only sensitive to wavefield intensity, specialised optical set-ups are required to convert variations in wavefield phase and dark-field into measurable variations in intensity. A mathematical model of the behaviour can then be used to retrieve the phase and dark-field information, typically from several exposures taken with different optical settings. There are a number of set-ups that have demonstrated the ability to do this, beginning with crystal interferometry [3] in the 1960s. The advent of highly coherent synchrotron x-ray sources several decades later led to propagation-based phase contrast x-ray imaging, where the coherent wavefield self-interferes during propagation to reveal phase information in addition to attenuation effects [4, 5]. More recently a number of techniques have been developed that utilise high-frequency gratings to measure attenuation, phase and dark-field information. These techniques include grating interferometry [6, 7, 8, 9, 10], edge-illumination [11] and single-grating imaging [12, 13]. A generalisation of single-grating imaging is speckle-tracking [14, 15], where the grating is substituted with a object, typically a sheet of sandpaper, that creates a reference pattern that need not be periodic, and is typically random.

This manuscript focuses on the grating and speckle-based phase contrast x-ray imaging set-ups. These set-ups can extract an attenuation image, a differential phase contrast image (or, in some cases, two complementary directions of differential phase contrast), and a dark-field image (or directional dark-field images in some cases, when applied to a sample with aligned elongated sub-pixel features). These image modalities can be understood by considering a pencil beam, or beamlet, passing through the sample. This beamlet will decrease in total intensity if the sample is attenuating. Variations in the thickness or refractive index of the sample can alter the phase of the x-ray wave, and result in a change in the direction of beamlet propagation, as seen in Fig. 1a. Therefore, the phase signal can be extracted by looking at some kind of transverse shift in the intensity profile at a distance downstream (e.g. $z=\Delta$, where $z$ is the optic axis). If the beamlet is symmetrically increased in width, as a result of sub-resolution sample structures (Fig. 1b) [16], this is described as a scatter, visibility or dark-field signal, measured by looking at a reduction in the visibility of the illumination [17]. Visibility can be defined using Michelson’s definition $V=(I_{max}-I_{min})/(I_{max}+I_{min})$, where $I_{max}$ is the local maximum intensity and $I_{min}$ is the local minimum intensity [18], so that high visibility indicates that $I_{max}$ is much larger than $I_{min}$. A strong dark-field signal will be seen as a significant decrease in illumination visibility. The dark-field signal measured by a grating-based system originates largely from small angle x-ray scattering (SAXS) or incoherent scattering that is below the resolution of the imaging system, but can also result from edges [19] or focusing [20], beam hardening [21, 22] and other effects. In most imaging applications, the sample will introduce both a transverse shift and broadening of the beamlet at a distance $z=\Delta$, as seen in Fig. 1c. The complementary transverse directions (e.g. $x$, $y$) of phase contrast and directional dark-field would be measured by looking at transverse shifts and broadening in the direction shown in the diagram and the direction perpendicular to the page.

The intensity profiles shown on the right of each panel in Fig. 1 are reminiscent of a probability density function that experiences changes in mean (Fig. 1a) and standard deviation (Fig. 1b). The sum of x-ray photons within each pixel, measured as intensity, can be modelled as a probability density function since the integrated intensity from sample exit plane to downstream plane is constant for a beam propagating in vacuum. This paper therefore studies the Fokker–Planck equation, an equation typically used to describe the time-evolution of a probability density function, and applies it to describe the propagation-evolution of a beamlet in x-ray phase/dark-field imaging [23]. In a companion paper[24], we show how the Fokker–Planck equation is an extension to the transport of intensity equation [25] used in coherent x-ray optics and also provide a derivation from first principles. Here, we begin with a short probability-based introduction to the Fokker–Planck equation (Sec. 2), then show that it can be applied to better understand how the x-ray dark-field (Sec. 3.1) and differential phase (Sec. 3.2) signals can interact, looking at additional effects seen near edges and as a result of focusing. These two sections (3.1 and 3.2) look specifically at imaging using patterned illumination that is directly resolved (e.g. single-grid imaging with an attenuating [12, 26] or phase-shifting grating [27, 28], or speckle-tracking [14, 15, 29]). Given that attenuation is a simple scaling of the intensity signal, hence a multiplicative factor at the front of the Fokker–Planck equation, we leave out this effect for simplicity. The final section, 3.4, demonstrates that the results and approach described in the earlier sections can be extended to the cases where the illumination is analysed via a second grating or mask (e.g. grating interferometry or edge-illumination imaging).

2 The X-ray Fokker–Planck Equation

The Fokker–Planck equation, shown below, describes the evolution of a probability density function $p$ with time $t$. Changes can include a time-dependent drift in both mean ($D_{1}$, a transport-type energy flow) and/or standard deviation ($D_{2}$, a diffusive-type energy flow) [30],

[TABLE]

as shown in Fig. 2 a after a given time interval $t=\Delta$.

Since the path taken by an ensemble of photons will follow a probability density function and be seen as such when those locations are registered as intensity on the detector, we can use the Fokker–Planck equation to model changes in the measured x-ray intensity. The Fokker–Planck equation can also be seen as the transport of intensity equation (TIE) with the addition of a diffusive term [24]. In the grating-based phase contrast x-ray imaging methods mentioned above, a grating will typically split the beam into an array of ‘beamlets’, each of which will pass through the sample and propagate through free space to a detector. We can begin by looking at the one-dimensional intensity distribution measured downstream of a sample that is illuminated by an array of beamlets, with a transverse shift of each beamlet seen at $z=\Delta$ when the beamlet passes through a phase-shifting sample (Fig. 1 a), and a reduction in visibility of each beamlet seen at $z=\Delta$ when the beamlet scatters at small angles from sub-pixel features in the sample (Fig. 1 b).

To move from the time-evolution of a probability function (Fig. 2a, Eqn. (1)) to the case of spatial drift and diffusion of an array of x-ray beamlets in the presence of a sample (Fig. 2b, Eqn. (2), where $I(x)$ is a sinusoid), we should replace:

• •

Time $t$ with propagation distance $z$

• •

Probability density function $p$ with x-ray intensity $I$

(In the case shown in Fig. 2, the centred probability density function $p$ is replaced by an array of x-ray intensity beamlets $I$, described by a sinusoid)

• •

Drift velocity $D_{1}$ with the change in angle of x-ray wavefield propagation, determined by the sample-induced phase shift $\frac{1}{k}\frac{\partial\phi}{\partial x}$

• •

Diffusion coefficient $D_{2}$ with a value $D(z)$ that characterises the length scale of unresolved x-ray scattering at distance $z$.

This gives the ‘X-ray Fokker–Planck equation’ [24],

[TABLE]

Note that the conventional Fokker–Planck model will describe a diffusion that increases the spread of light in a non-linear manner with propagation. This is not consistent with the linear increase in spread with propagation distance seen when a scattered wave-field propagates through free space. For this reason, $D(x,z)$ is set as a function of $z$, so that the diffusion width can increase linearly with propagation distance $z$ (see Paganin & Morgan for further justification[24]).

In order to fully describe directional, three-dimensional or more complex behaviour, we would need an equation that includes more spatial dimensions and terms. This is a good time to note that the Fokker–Planck equation is the truncated form of the more general Kramers–Moyal equation [23, 30]. This more general equation includes third-order and higher-order terms, and describes changes in $n$ spatial dimensions, with ${\bf{x}}=(x_{1},x_{2},x_{3}...)$ as Cartesian coordinates, as below,

[TABLE]

For the case of two dimensions and truncation immediately before the ‘$-...$’ shown above, we have the anisotropic-diffusion Fokker-Planck equation, which well-describes a two-dimensional x-ray beamlet landing on a two-dimensional detector,

[TABLE]

(absorbing any multiplying fractions into $D_{x},D_{y}$ and $D_{xy}$). This kind of diffusion in multiple directions is ideal for modelling a directional and/or irregularly shaped x-ray dark-field signal[31], and is an avenue for future work. The parameters $D_{x},D_{y}$ and $D_{xy}$ describe a diffusion profile that has the shape of an ellipse, and hence requires three parameters to describe it fully. These three parameters may be mapped one-to-one to the semi-major axis, the semi-minor axis, and the angular orientation of the elliptical diffusion profile. Note also, that if the infinitely-many terms in the Kramers–Moyal equation are truncated to no higher than second-order spatial derivatives, the result is the Fokker–Planck Equation. As shown by Pawula, one should either truncate to second-order or not at all, since ‘if one assumes that terms above a given order are zero in the Kramers–Moyal expansion, this assumption implies that all terms above second order are zero’[32].

3 Application of the Fokker–Planck model to describe phase and dark-field effects

To model x-ray grating-based phase and dark-field imaging, the key method considered in the present paper, we can set the intensity that illuminates the sample to a sinusoid. This can describe the intensity that would be seen downstream of a grating or grid, and can be specified for period $p$, amplitude $a$ and mean intensity $b$,

[TABLE]

Note that we have omitted the factor of $2\pi$ inside the sinusoid that would be required to describe a period of exactly $p$. This choice was made to keep the equations as simple as possible without neglecting any effects. The general expression given in Eqn. 5 means the illumination without any sample in the beam would have a mean intensity of $b$ and visibility of $a/b$ (e.g. the grey curve in Fig. 3 a). Given that a Fourier series could be used to form any other illumination (e.g. periodic ‘square’ illumination sampled an order of magnitude more finely than the grid period [26] or a speckle pattern generated by a piece of sandpaper [14, 15, 29]), the choice of a sinusoidal illumination is generally useful. In the illustrations given in Figs. 3-6, we have sketched a sinusoidal illumination that is likely seen at a detector downstream of an absorption or phase grid, effectively a normalised intensity (reaching from 0 to 1) that has been smoothed to a more realistic visibility of 0.5 so that the sinusoid is described by $a\approx\frac{1}{4}$ and $b\approx\frac{1}{2}$. A visibility that is less than 1 is typically seen due to limited coherence, some penumbral blur and/or sampling of relatively fine grids (e.g. 7 pixels per grid period) that are chosen to maximise the resolution of the retrieved images [12].

To examine measurements at a given detector position, we can take the x-ray Fokker–Planck equation (Eqn. (2)) with the finite-difference approximation for propagation to $z=\Delta$ (functional dependencies, where obvious, have been dropped for simplicity and clarity),

[TABLE]

In the following sections we will consider how these two steps (Eqns. (5) and (6)) can be applied to the dark-field/diffusion term (Section 3.1) and the phase/drift term (Section 3.2) to show how the Fokker–Planck equation is relevant to grating-based x-ray imaging and to understand how spatial variations in the phase and dark-field may introduce signal variations that could be ‘falsely’ interpreted as attenuation, phase or dark-field signals. The Supplementary material includes python code for interactive visualisation of the less-intuitive mathematical expressions from Section 3, which can be used to more intuitively understand the effects described herein (also found at https://github.com/KayeMorgan7/XrayFokkerPlanckSim).

3.1 Dark-field / Scattering effects

We look first at an object with scattering properties, but no phase-shifting properties ($\partial\phi/\partial x=0$), which will leave only the second term on the right side of the x-ray Fokker–Planck equation (Eqn. (2)). We can insert the sinusoidal illumination $I(x,z)$ (Eqn. (5)) and take the finite-difference approximation for Eqn. (6) to leave,

[TABLE]

3.1.1 Uniform scattering

If we consider an area where the diffusion properties do not vary quickly with $x$ (so that $D$ can be moved outside of $\frac{\partial^{2}}{\partial x^{2}}$), we can easily differentiate the sinusoid twice with respect to $x$, and are left with

[TABLE]

Equation 8 describes a decrease in the amplitude of the sinusoid, and hence the visibility of the sinusoid, as seen in Fig. 3 a), which is indeed the signal which is collected in dark-field imaging. A more significant effect will be seen for a larger diffusion coefficient $D$, larger propagation distance $\Delta$ or a smaller grid period $p$, as expected.

Evaluating a direct solution like this can encounter problems if the amplitude of the sinusoidal illumination becomes negative. The problem arises because a linear approximation has been made with propagation distance $\Delta$ (Eqn. (6)), and the linearisation eventually reaches zero and then becomes negative. Equation (8) will therefore be valid only when $D\Delta<p^{2}$. The case when the intensity sinusoid will reach an amplitude of zero is at $D\Delta/p^{2}=1$, when the dark-field blurring is over a length scale at the detector $L$ that is equal to the period of the sinusoid $p$. This naturally provides a definition for the ‘blur width’ $L$, $L=\sqrt{D\Delta}$, matching the definition found in the companion paper[24]. Note that if we incorporate $D(z)=D_{0}\Delta$, then $L=\sqrt{D_{0}}\Delta$.

To avoid the limited applicability in any numerical simulations, an alternative solution to the diffusion equation is via a linear integral transform (which can loosely be considered as convolution with a position-dependent point spread function, where the point spread function in our case is a Gaussian where the standard deviation grows with $\sqrt{D\Delta}$)

[TABLE]

For sinusoidal illumination, this gives,

[TABLE]

Note the standard deviation of $\sqrt{D\Delta}$ used in Eqn. 9 could be scaled with a factor of $\sqrt{2}$ to better match the blurring width used in the finite-difference approximation so that Eqns. 8 and 10 match at low propagation distance $z=\Delta$, but the approach here is sufficient for the order of magnitude effects in modelling the diffusive behaviour.

3.1.2 Scattering edge effects

If the diffusion property $D$ does change quickly with $x$, such that $D$ cannot be moved outside $\frac{\partial^{2}}{\partial x^{2}}$ in Eqn. (7), then we start to see more interesting behaviour. This case is most likely to correspond to an edge of a strongly-scattering region. We can expand the brackets from Eqn. (7) to get

[TABLE]

We can better understand the effects by using a trigonometric identity to express Eqn. (11) as a single sinusoid (i.e. $A\sin(x)+B\cos(x)=\sqrt{A^{2}+B^{2}}\sin\left(x+\arctan(B/A)\right)$, noting that here we also take a small angle approximation for $\arctan$):

[TABLE]

This is the diffusion effect we saw in Eqn. (8), with some additional effects near an edge, where $\frac{\partial^{2}D}{\partial x^{2}}$ is significant, as shown in Fig. 3 b. In any edge regions, we see that the amplitude of oscillations in the periodic illumination can now be increased or decreased and that the local mean intensity ($b$) can be increased or decreased, since it is multiplied by $\left(1+\Delta\frac{\partial^{2}D}{\partial x^{2}}\right)$. Because the local first and second spatial derivatives can change quickly (e.g. on the edge of a scattering object), this is a local effect that can result in an illumination that is no longer well-described by a simple sinusoid.

This can be explored further by looking at two neighbouring regions with differing scattering properties, setting

[TABLE]

so that $D=D_{right}$ far to the left of the origin, and $D=D_{left}$ to the right of the origin, with some smooth transition between, described by the Erf function (i.e. the special function known as the ‘error function’).

This gives:

[TABLE]

As mentioned in the final paragraph of section 3.1.1, the applicability of a solution found using a finite-difference approximation (Eqn. (6)) is limited to the case where the area of diffusion is smaller than the period of the illumination pattern. Using the alternative approach given in Eqn. (9), we can evaluate first the diffusion on the left of the origin, and then the diffusion effect on the right, with the intensity given by Eqn. (5) ($I(x^{\prime})=a\sin{\left(\frac{x^{\prime}}{p}\right)+b}$):

[TABLE]

This evaluates to:

[TABLE]

Equation 3.1.2 can correctly describe effects seen at an edge without the problems associated with the potential for a ‘negative’ sinusoidal solution given in Eqn. (3.1.2). This means that the diffusion effect can be modelled when it extends over several periods of the illumination. Several cases are shown in Fig. 4 b, varying the illumination period $p$ only, and these cases help to develop an understanding for the ‘diffusion edge effect’. The python code in the Supplementary Material implements Eqn. 3.1.2, providing a GUI to explore this effect (also at https://github.com/KayeMorgan7/XrayFokkerPlanckSim). These results indicate that the strongest effect is the variations in the mean value (the last term in Eqn. (3.1.2)), resulting in a signal, shaped like the second derivative of the diffusion coefficient, that sits ‘on top’ of the sinusoidal illumination. This origin of this intensity variation can be described by considering the areas where ‘additional scattered light’ and ‘less scattered light’ are observed, as shown in Fig. 4 and described in the caption.

Note, however, that with this solution (Eqn. (3.1.2)) the dark-field signal ($D_{left},D_{right}$) cannot be set to zero as this will cause problems associated with dividing by zero, but can be set to a very small number when required.

3.1.3 Modelling scattering using a dark-field coefficient

It is of interest to have a way to characterise the diffusion characteristic that is consistent with the Fokker–Planck model. To reach this, we can look at the original differential equation (Eqn. (2)) without phase effects, $\phi$=0. This gives the diffusion equation:

[TABLE]

If we have any periodic illumination (Eqn. (5)), and we consider an area over which the diffusion/scattering is not changing significantly with $x$ (so that $D$ can be moved outside of $\frac{\partial^{2}}{\partial x^{2}}$), then Eqn. (17) must satisfy:

[TABLE]

This differential equation can be solved as

[TABLE]

which agrees with the ‘linear dark-field co-efficient’ of $D/p^{2}$ seen in Lynch et al. [33], and used in dark-field computed tomography [34].

3.2 Phase effects

Looking at a sample with phase properties, $\phi\neq$ 0, independently of dark-field effects (setting $D=0$), from Eqn. (2), we can approximate the intensity observed at distance $z=\Delta$ by the finite-difference approximation (Eqn. (6)), insert sinusoidal illumination (Eqn. (5)) and evaluate the derivative with respect to $x$ to give

[TABLE]

Using the trigonometric identity for the addition of two sinusoidal signals of the same argument (provided immediately before Eqn. (12)), and leaving out terms of higher order $\Delta$ (since we linearised with $\Delta$ in Eqn. (6)), Eqn. (20) simplifies to:

[TABLE]

3.2.1 Prism effect, seen for uniform phase gradient $\frac{\partial\phi}{\partial x}$

If we are in the slowly-changing bulk of the sample, where $\frac{\partial\phi}{\partial x}$ is significant, but $\frac{\partial^{2}\phi}{\partial x^{2}}$ is small and can be neglected, Eqn. (21) simplifies to:

[TABLE]

This means that if the phase shift is linear (so $\frac{\partial^{2}\phi}{\partial x^{2}}$ is zero), then we see only the prism effect, where the grid pattern shifts transversely. The shift in the relative position is what we would expect (see Morgan et al.[13]), namely $\frac{\Delta}{k}\frac{\partial\phi}{\partial x}$), and there is no change in the mean or visibility of the signal. This is the ‘well-behaved’ case, where phase effects are not introducing any signals that look like attenuation or dark-field.

3.2.2 Focusing effect, seen for $\frac{\partial^{2}\phi}{\partial x^{2}}\neq 0$

If the phase gradient is not uniform, and varies with position $x$, then we will have more interesting effects. We can consider each of the terms in Eqn. (21) to understand the physical meaning. With soft x-rays and a propagation distance in the near field, $\frac{\Delta}{k}$ will be very small (e.g. $10^{-11}$). If the second phase derivative is large, i.e. approaching the order of $\frac{k}{\Delta}$, then the multiplicative term at the front of the sinusoid in Eqn. (21) can be slightly bigger or smaller than 1, describing either some amplification or deamplification of the entire illumination. Whether amplification or deamplification is seen will depend on the sign of the second derivative of the phase, which describes either defocusing or focusing of the illumination. We also see that if the second phase derivative approaches the order of $\frac{k}{\Delta}$, then the denominator inside the $\tan^{-1}$ in Eqn. (21) will become slightly bigger or smaller than $\frac{k}{\Delta}$, describing either a squeezing or a spreading out of the intensity peaks, matching with the defocusing/focusing behaviour.

To untangle the effect of phase gradients on the period of the illumination, from the effects of phase gradients on the shift of the illumination (i.e. the $x$-dependent terms within the sinusoid), we can model the case of some ‘phase curvature’ using a Taylor polynomial of order two, $\phi=Ax^{2}+Bx+C$, so that $\frac{\partial\phi}{\partial x}=2Ax+B$ and $\frac{\partial^{2}\phi}{\partial x^{2}}=2A$. The values for $A$, $B$ and $C$ could be chosen to describe the phase gradient and curvature of any local area of the sample. With this phase, Eqn. (21) becomes:

[TABLE]

with the simplification taken for small shifts in phase of the sinusoid, so that $\tan(x)\approx x$.

We see that a linear change in phase (modelled by $A=0$, $B\neq 0$) will shift the relative position of the peaks of the periodic illumination by $\frac{B\Delta}{k}$, acting as a prism without changing the period of the illumination, as seen in Eqn. (22), and on the left third of Fig. (6).

If the phase shift incurred by the sample is not just linear, but there is such curvature that a quadratic term is required to adequately describe the local phase, then we will see additional effects. This case corresponds to when $A$ is significant, approaching $\frac{k}{2\Delta}$, so that the term $1-\frac{2A\Delta}{k}$ moves away from 1. Since this term appears in three places, a large $A$ will 1) change the mean illumination, 2) alter the period and 3) introduce an additional shift in the position of the periodic illumination. If the phase curvature is positive (seen for $A>0$ in our quadratic model) then the beamlets will spread out, seen by a decreased mean illumination and a period that is multiplied by a number slightly greater than 1 and so locally increased. The shift in the position of the illumination will also be slightly reduced towards the centre of the focus. If the phase curvature in negative ($A<0$), the illumination will be focused, increasing the mean intensity and decreasing the period of the periodic illumination, as shown in Fig. 6. These observations are consistent with preservation of photons. While this kind of focusing will not introduce any visibility changes when the grid pattern is directly resolved, it can change the visibility measured in grating interferometry [35, Fig.1] (See section 3.4.1).

3.2.3 Phase edge, seen for $\frac{\partial^{2}\phi}{\partial x^{2}}\neq 0$

A phase edge is also an interesting feature to study, particularly given the bright/dark contrast seen in propagation-based phase contrast imaging when a wave passes across an edge and self-interferes. If we model a phase edge in the same way as Eqn. (13) models a scattering edge (replacing $D$ with $\phi$) and use this model in our expression for phase effects, Eqn. (21), we get

[TABLE]

When plotted, as expected, a bright/dark contrast fringe overlies the illumination, in a similar manner to the dark-field edge. The python code in the Supplementary Material implements Eqn. 24, providing a GUI to explore this effect (also at https://github.com/KayeMorgan7/XrayFokkerPlanckSim). This can result in local changes to the visibility, as seen in Fig. 6.

3.3 Combined effects

While the phase effects and dark-field effects are most easily understood in the previous two sections, where each is isolated, we can also write an equation to describe the combined effects. Since the shifting and diffusive terms in the finite-difference pproximation for the Fokker–Planck equation (Eqn. (6)) add together, we can describe the resulting intensity as the sum of Eqn. (12) and Eqn. (21), subtracting the illuminating intensity (Eqn. (5)) since it is in both equations and hence repeated. This gives:

[TABLE]

3.4 Extension of the Fokker–Planck model to Grating Interferometry and Edge Illumination

The results presented in the preceding sections can be extended to include the effect of an analyser grating or detector mask, as would be used in a grating interferometer or edge illumination set-up respectively. The use of a second grating/mask has the advantage that larger pixels can be used, which is typically more efficient and allows the imaging of larger samples than imaging where the illumination pattern is directly resolved. To incorporate the second grating/mask, we can simulate a periodic absorption immediately before the detector to ‘analyse’ the illumination we have studied so far.

Here we show the implementation of this, using a general signal that describes any of the cases we have explored (Eqn. (3.1.2), Eqn. (22), Eqn. (23), Eqn. (24) etc.), where the sinusoidal illumination arriving at the periodic absorber has been modified by a factor of $M$ in mean value, a factor of $V$ in visibility, and has been shifted by $s$, so can be described as

[TABLE]

Note that we neglect a local change in period here, due to the increased complexity, but the same approach can be applied to study this effect. Also note that we have explicitly included a factor of $2\pi$ inside the sinusoid in Eqn. 27 that is not included earlier, so that it is possible to integrate the unresolved intensity oscillations over the period of that intensity oscillation, which will now correspond exactly to $p$. We can model the absorption by the analyser grating or detector mask by multiplying the illumination (Eqn. (26)) by a sinusoid that spatially oscillates between 0 (total absorption) and 1 (total transmission) (i.e. $\frac{1}{2}\sin\left(\frac{2\pi(n-x)}{p}\right)+\frac{1}{2}$).

3.4.1 Grating Interferometry

In grating interferometry, the intensity is recorded during a ‘stepping scan’, where an absorption grating steps transversely across the illumination. This can be implemented by convolving the ‘absorption’ sinusoid with the modified illumination sinusoid (Eqn. (26)) over the period of the illumination, $p$, providing the infinitely sampled measured stepping curve as a function of absorption grating position $n$ (which should be sampled from $n=0$ to $n=p$),

[TABLE]

This stepping curve has a mean value of $Mp/2$, a visibility of $V/2$ and is shifted by $s$. This means that the results shown above for a directly-resolved grid can also be applied in the case of an interferometer with just a simple multiplicative scaling of $p/2$ for the mean (attenuation), $1/2$ for the visibility (dark-field) and no adjustment for the shift.

Note that if the period of oscillations is locally varying (seen when $\frac{\partial^{2}\phi}{\partial x^{2}}\neq 0$, Sections 3.2.2 & 3.2.3), we would need to scale the period $p$ by a factor of $f$ in the general expression for a modified sinusoid, Eqn. (26). Carrying this into Eqn. (27) produces a non-trivial solution that depends on the nature of $\frac{\partial^{2}\phi}{\partial x^{2}}$, in a similar way to the results seen in Sections 3.2.2 and 3.2.3, where we had to consider specific cases of $\frac{\partial^{2}\phi}{\partial x^{2}}$. Including spatial variations in illumination period would be important when considering how a strong phase object measured using interferometry may produce a ‘pseudo’ dark-field (e.g. [35, Fig. 1]).

3.4.2 Edge Illumination

Edge illumination imaging has several variants, for example, depending on whether there is strong scattering expected from the sample[11, 17]. With the aim only to illustrate the applicability of the Fokker–Planck approach, here we will consider just the simple case where scattering is weak. In this case, the absorption grating / detector mask, which has a period matching the illumination period, should be aligned so that the peak of the sinusoidal illumination first coincides with the left edge of the each absorption line (exposure 1, $I_{left}$), then coincides with the right edge of the absorption line (exposure 2, $I_{right}$),

[TABLE]

The shift ($s$) and attenuation ($M$) signals could then be extracted by adding and subtracting Eqns. (29) and (30), giving consistent results as in the scenario described above for grating interferometry.

Following these demonstrations of applicability, we leave further applications of the Fokker–Planck equation to Grating Interferometry and Edge Illumination to future work. Possibilities include integrating over a non-integer number of grating lines and investigating the effects of a locally-modified period.

4 Discussion

This work serves as a starting point for future applications of the Fokker–Planck and Kramers–Moyal equations in x-ray imaging techniques that capture phase and dark-field effects. We begin the discussion by mentioning a number of possibilities for future work. A clear next step in any grating-based or speckle-based approach is to utilise this model as a starting point for phase/dark-field retrieval, or simply as a basis for lensing and edge-based ‘corrections’ to retrieved attenuation, phase and dark-field signals (perhaps as a second or iterative step after the initial extraction of the three signals by conventional retrieval methods). As noted in Section 3.4.1, the Fokker–Planck model could be applied to provide an analytical model for specific cases seen in grating interferometry where phase edge[19] and lensing effects (simulated and measured in Wolf et al.[20]) result in an illumination period that is locally-varying, leading to a spatially-varying ‘mis-match’ between that illumination and the second/analyser grating[35].

Corrections for dark-field-edge, phase-edge and lensing effects are particularly of interest when the size of the features within the sample are comparable to the period of the illumination, and the effects are typically more dominant. The Fokker–Planck model described in this paper could better shed light on the transition between two regimes, the regime where sample features are large enough to create a spatially-resolved phase signal (e.g. see differential phase images of individual air sacs in the lung in Morgan et al. [36] or high-resolution lung CT in Koenig et al. [35, Fig.5]) and the regime where sample features are not spatially resolved and are instead seen as a dark-field signal (e.g. see dark-field images of the lungs[37] or [35, Fig.2]). The Fokker–Planck model naturally lends itself to a parameter-based definition for these regimes, differentiating between edge effects and bulk effects, as determined by the spatial resolution of the imaging set-up and the illumination period. This direction of future work also lends itself to an investigation into the optimum illumination pattern period to be sensitive to a given dark-field signal for a given wavelength with a given imaging set-up. This links to the correlation length, $\lambda\Delta/p$ seen in grating interferometry [38], which matches the expression seen when the ‘smear width’ $L$ is set to the illumination period $p$, zeroing the visibility of the resulting illumination at the detector (also see Fig. 3 c) of the companion paper on the X-ray Fokker–Planck equation [24]).

There are also natural extensions to the Fokker–Planck model. The first extension would be to describe more detailed[39] and size-related[33, 40] scattering functions. A dark-field signal from a non-uniform collection of overlaid and unresolved sample features could be modelled by the addition of extra diffusive terms ($\partial^{2}(D_{m}I)/\partial x^{2}$, each with a different $D_{m}$) to the Fokker–Planck equation (Eqn. (2)). Since the diffusion equation is energy-conserving, the addition of these terms would require no other adjustments. A dark-field signal from a collection of extended unresolved sample features (e.g. carbon fibres [41, 42] or bone microstructure [31, 43]) could be modelled by a two-dimensional Fokker–Planck equation (a truncation of the Kramers–Moyal equation), which provides three parameters to describe the two-dimensional scattering ellipse utilised in directional dark-field imaging [44] (which is sometimes modelled as scattering that oscillates in magnitude sinusoidally with angle [31, 45]).

The Fokker–Planck model shown here could be extended to incorporate real-world effects seen at synchrotron and laboratory x-ray sources, including limited spatial resolution, extended sources, polychromatic illumination and magnification. The spatial resolution limit due to the detector characteristics, including pixel size, could be modelled by convolving the general expressions derived above (e.g. Eqn. (25)) with a Gaussian of specified width. The spatial resolution associated with an x-ray imaging system will also be limited by the size of the x-ray source, with an extended source ‘blurring’ the observed image. In the Fokker–Planck model, an extended source could be modelled by summing the images created by an array of many point sources, produced by shifting the intensity measured at $z=\Delta$ slightly in order to model the grid projection seen from each point on the source. An alternative would be to replace $D(x,y)$ with $D(x,y)+d$ (or equivalently, add another diffusion term), where $d$ models the extra blur that results from source-size blurring (see p. 13 of Paganin & Morgan [24]). The effect of an extended source is also of particular relevance in the measured x-ray dark-field signal [35, Sec.5]. In the direction of less blurring, a more finely-resolved structured illumination behind a rectangular (or some other-shaped [46]) grid could be modelled using a Fourier sum of the sinusoidal illumination model utilised in this paper.

Any x-ray source that is not a synchrotron will typically produce illumination that is sufficiently polychromatic and divergent that the associated effects must be included in modelling. A divergent source could be modelled according to the Fresnel scaling theorem [47], by setting the propagation distance $z$ to $z/M$, where $M$ is the magnification, while adjusting the native pixel size accordingly. This is also key in x-ray dark-field imaging by interferometry [48] [35, Sec.4]). An image captured with polychromatic illumination could be modelled as the weighted sum of the images seen at each wavelength within the illuminating spectrum [49]. This can be implemented by integrating the Fokker–Planck equation (Eqn. (6)) over wavelength, after multiplying Eqn. (6) with a wavelength-dependent weighting factor that describes the x-ray spectrum and which could also incorporate detector sensitivity.

A sum of many slightly-different images into a single exposure can be used to describe a reduction in fringe or patterned illumination visibility for a number of cases. In addition to the cases of an extended source and a polychromatic source, reduced visibility is seen when using a diffuser. A diffuser is typically a sheet of highly-scattering material, like sandpaper, which is placed in the x-ray beam upstream of the sample and rotated at high speed, such that many positions are summed over the time of a typical x-ray exposure. This will reduce the visibility of propagation-based x-ray phase contrast fringes, as each diffuser position slightly distorts the fringes, and the sum over all diffuser positions results in a fringe that is ‘blurred’ [50]. A diffuser will not change the spatially coherent width/area, but will reduce the magnitude of the second order degree of coherence [51]. In x-ray SAXS, the wavefield is locally distorted by scattering from sub-pixel structures in the sample, locally reducing the spatially coherent width/area, and resulting in a reduction in the visibility of fine features like the illuminating grid pattern. In the case of a diffuser, there is a temporal sum of many slightly distorted images. In the case of x-ray SAXS, there is a temporal sum of many photons, where some are slightly scattered.

Another aspect of dark-field is worth discussing. Several papers have reported on an enhancement of dark-field in the vicinity of the boundary of the geometric shadow of a paraxially-illuminated compact object [52, 53, 19]. Rather than being dark-field due solely to unresolved micro-structure, as is the case for points that lie inside the geometric shadow, the dark-field signal emanating from the edge of the object is augmented by the Young–Maggi–Rubinowicz boundary diffraction wave[54]. This phenomenon, whereby illuminated objects have their projected edges acting as a secondary source of (scattered) boundary waves, was first noted in a qualitative study by Young in 1802[55]. Subsequent contributions to the underpinning wave theory are due to Maggi[56], Rubinowicz[57] and Miyamoto & Wolf[58, 59]. The same phenomenon may also be understood from the perspective of complex-ray theory. This ray-based perspective underpins Keller’s geometric theory of diffraction[60], in which Fresnel diffraction patterns may be understood as the interference between (i) the complex-ray disturbance transmitted by an object according to geometrical optics, and (ii) the complex-ray disturbance scattered from the edges of the object. A third means of understanding the dark-field-enhancement along the boundary of the geometric shadow of an illuminated object, is from the asymptotic expansion (short-wavelength limit) of the two-dimensional diffraction integral associated with a paraxially-illuminated compact object under the projection approximation[47], obtained by applying the method of stationary phase for double integrals. Here, the boundary wave corresponds to critical points of the second kind: see e.g. Sec. 3.3.3 of the text by Mandel & Wolf[30]. For applications of the boundary-wave concept to coherent X-ray optics, see e.g. Morgan et al.[26, 61].

5 Conclusion

This work has shown that the Fokker–Planck equation can be applied to phase contrast and dark-field x-ray imaging captured using a single-grid set-up, a grating interferometer, an edge-illumination set-up or a speckle-based set-up. The Fokker–Planck model was applied to provide an analytic expression for the additional signals seen as a result of edges within the sample and lensing by the sample. The authors hope that this work will provide a basis for further applications of the Fokker–Planck model in x-ray phase and dark-field imaging.

Acknowledgements

We acknowledge useful discussions with Florian Schaff, Mario Beltran, Carsten Detlefs, Timur Gureyev, Alexander Kozlov, Thomas Leatham and Tim Petersen. We acknowledge the following funding sources: ARC Future Fellowship FT180100374; Veski Victorian Postdoctoral Research Fellowship (VPRF); German Excellence Initiative and European Union Seventh Framework Program (291763). The python GUI is built on the slider demo from matplotlib, available at https://matplotlib.org/3.1.1/gallery/widgets/slider_demo.html.

Author contributions statement

K.S.M. and D.M.P. worked on this paper in close collaboration. K.S.M. prepared all figures and performed all mathematical calculations, arising from discussions between both authors. The paper was mainly written by K.S.M., with input from D.M.P.

Additional information

Competing interests: The authors declare no competing interests.