Review and experimental verification of X-ray darkfield signal interpretations with respect to quantitative isotropic and anisotropic darkfield computed tomography

Talbot(-Lau) interferometric X-ray darkfield imaging has, over the past decade, gained substantial interest for its ability to provide insights into a sample's microstructure below the imaging resolution by means of ultra small angle scattering effects. Quantitative interpretations of such images depend on models of the signal origination process that relate the observable image contrast to underlying physical processes. A review of such models is given here and their relation to the wave optical derivations by Yashiro et al. and Lynch et al. as well as to small angle X-ray scattering is discussed. Fresnel scaling is introduced to explain the characteristic distance dependence observed in cone beam geometries. Moreover, a model describing the anisotropic signals of fibrous objects is derived. The Yashiro-Lynch model is experimentally verified both in radiographic and tomographic imaging in a monochromatic synchrotron setting, considering both the effects of material and positional dependence of the resulting darkfield contrast. The effect of varying sample-detector distance on the darkfield signal is shown to be non-negligible for tomographic imaging, yet can be largely compensated for by symmetric acquisition trajectories. The derived orientation dependence of the darkfield contrast of fibrous materials both with respect to variations in autocorrelation width and scattering cross section is experimentally validated using carbon fiber reinforced rods.


Introduction
A Talbot(-Lau) grating interferometer is a specific realization of a shearing interferometer based on the Talbot-effect. It can be implemented for the X-ray spectrum and, in addition to classic attenuation contrast, further gives access to diffraction contrasts in the form of differential phase shift information and ultra small angle X-ray scattering contrast (cf. [7,8,42,64,50]). The latter is commonly referred to as "darkfield contrast" in analogy to the respective scattering contrasts in other fields of imaging and is of particular interest for its sensitivity to the unresolved substructure of the sample. The unique advantages of the grating interferometer with respect to other diffractive X-ray imaging techniques, such as scanning small angle X-ray scattering (scanning SAXS), the crystal analyzer based "diffraction enhanced imaging" or "multiple image radiography" methods (cf. [76,45,47,67]), and speckle-based imaging techniques (cf. the recent review by Zdora [75]), are its ability to directly capture planar images (in contrast to pixel or line scanning methods), its sensitivity to micrometer scaled diffractive effects within centimeter sized fields of view (in contrast to microscopy techniques) and finally its tolerance to non-monochromatic X-rays and practical feasibility in laboratory environments (as opposed to synchrotron facilities). The darkfield contrast generated by Talbot interferometers (along with attenuation and differential phase contrast) provides complementary information both in non-destructive testing and life sciences due to its sensitivity to sub-resolution structures such as micro cracks and porous or fibrous matter. The initial . Coherent X-rays are modulated by a periodic, phase-shifting or absorbing grating G 1 . Coherent diffraction within the Fresnel-regime leads to a periodic reproduction (in intervals of the Talbot distance) of the periodic amplitude modulation induced by the G 1 grating (Talbot-effect). The patterned beam can be analyzed by means of a period-matched absorption grating G 2 in front of an integrating detector. The detector pixel size is typically at least one order of magnitude larger than the grating period. By moving the analyzer grating G 2 in front of the detector about one period in multiple steps, a phase stepping curve can be acquired (cf. right panel. Stepping of G 1 will have an equivalent effect). Wavefront changes due to a sample will result in attenuation, deflection or blurring of the G 1 -induced beam modulation and results in corresponding changes of the phase stepping curve. Samples may be placed between G 1 and G 2 or in front of G 1 . Experimental realizations often use an additional absorption grating G 0 before G 1 , structuring common laboratory X-ray sources into multiple narrow slit sources of sufficient coherence (Talbot-Lau interferometer).
interest on darkfield imaging arose from early synchrotron experiments on mammography with monochromatic radiation and crystal analyzers [25,6]. With the introduction of Talbot-Lau imaging to the laboratory by Pfeiffer et al. [50,49], many more examples have been shown. Applications of X-ray darkfield imaging besides the characterization of micro-calcifications in mammography [38,15] include imaging of lungs [71,61,17,32], characterization of bone and dentin [51,43,27,26,21] as well as the analysis of general porous and fibrous materials or microscopic defects in the field of non destructive testing [54,24,16]. Other applications that have been shown include water transport in cement [52,70] or the monitoring of germinating seeds [44]. Due to its origin in the ultra small angle scattering properties of a given sample, darkfield contrast also reflects anisotropies in scattering and thereby allows to detect local orientations within fibrous materials (cf. [22,23,51,3,56,57,19]). With respect to quantitative interpretations of the obtained images, a solid understanding of the underlying contrast mechanisms is of great importance, particularly in the light of long processing chains beginning with phase stepping analysis (cf. [58,29,37,10,20] and references therein) and ending in tensor valued volume reconstructions of anisotropically scattering materials (cf. [2,35,62,68,11,9]). The aim of the present article is to provide a unified view on the numerous explanations on darkfield signal origination that have been given in previous literature, and to provide experimental support for the central results. The existing theories will further be extended to cone beam geometries using the Fresnel scaling relation. Moreover, a model for the description of generally oriented anisotropic scatterers will be derived based on the concepts presented by Yashiro et al. [74] and Lynch et al. [33]. Experimental support will be given likewise.

Talbot Interferometer
Generally, an X-ray Talbot-Lau interferometer consists of three micrometer-pitched gratings, of which the first is placed close to the X-ray source, shaping it into multiple small slit sources in order to increase coherence of the emitted radiation. It can be omitted for sufficiently coherent sources such as synchrotron radiation or microfocus X-ray tubes (Talbot interferometer, cf. Fig. 1). The second grating imposes a periodic phase or amplitude modulation. Although the coherent wavefront is subject to interference while propagating, its periodic intensity modulation is restored at specific distances characteristic to the grating pitch and X-ray energy (Talbot effect). At these distances, it can be analyzed by means of a third grating with matching periodicity in combination with an X-ray detector placed behind it (having a pixel size considerably larger than the grating period). When the structured beam is perturbed by a sample, three effects can be observed: As for classic X-ray imaging, the intensity may be diminished due to absorption. Moreover, the periodic pattern may be reduced in contrast (visibility) due to scattering, and shifted in phase due to refraction. The effects of a sample can be analyzed by comparison of the phase stepping curves for the perturbed (by the sample) and unperturbed beam.

Isotropic Darkfield Contrast
Analog to the characterization of the modulation transfer properties of an optical system, visibility or contrast reduction of the structured beam (i.e., the reduction of the first harmonic Fourier component normalized to the zeroth as e.g. defined by [49]) can generally be modeled as Gaussian blurring of the reference sinusoid profile. The underlying conception is a convolution of the incident intensity distribution with an effective point spread function or scattering profile of a sample. Normalizing the Gaussian to the ratio of transmitted intensity t ∈ [0, 1] (i.e., also accounting for attenuation or general reduction of intensity), the following convolution kernel can be stated: with σ φ being its standard deviation in units of radians and ∆φ the phase with respect to the periodic irradiation profile (not to be confused with a scattering angle). The width σ φ will later be expressed in terms of geometric parameters based on absolute scales of the instrument. Yet to begin with, the present formulation will be convenient. By explicitly solving the convolution integral over a generic sinusoid reference and identifying the changed mean intensity and amplitude parameters o smp and a smp on the right hand side, the sample-induced visibility reduction v, defined Due to the orthogonality of the Fourier basis, this result also applies to non-sinusoidal phase stepping curves and generally describes the contrast visibility of a particular harmonic of the periodic pattern implicitly selected by the respective definition of φ with respect to spatial dimensions, which will be concretized in the following.

Linear diffusion interpretation
Given the period T G2 of the Talbot interferometer's analyzer grating, the sample induced blurring width σ φ of the first harmonic can be expressed in spatial units (indicated by the subscript x): When further arguing that the specific width σ x at the location of the analyzer grating arises due to a beam divergence caused by the sample at distance d, the respective divergence may as well be characterized by an angular standard deviation σ θ . Using the small-angle approximation ∆θ ≈ tan ∆θ = ∆x d for d ∆x, we find: Within the model of angular beam diffusion and geometric propagation onto the analyzer grating, σ θ represents an invariant property of the sample (at a given X-ray energy), while v, σ φ and σ x include properties of the instrument. This geometric interpretation has e.g. been given by Wang et al. (2009 [63]) to describe the origination of darkfield contrast, and Grünzweig et al. (2010/2013 [4, 18]) introduced this concept as "linear diffusion coefficient" further normalizing σ 2 θ to the sample thickness ∆z in order to obtain a thickness independent material constant analog to the classic "linear absorption coefficient" µ, such that: with z denoting the optical axis and x, y being the planar projection image coordinates. Anticipating the following Sections, it shall be noted here that the angular diffusion interpretation as presented here only holds in a parallel beam geometry and in the limit of sufficiently large and smooth scattering structures. The critical length scale depends on the interferometer parameters, as will be explained in the following (cf. Sections 3.2 and 3.3). Additional geometric scaling effects beyond the trigonometric relation σ x ≈ dσ θ need to be explicitly accounted for in the case of a Talbot interferometer in cone beam geometry, i.e., in the typical laboratory use case (cf. Secion 3.4). Lynch et al. (2010-2011 [74, 33]) derive the origination of darkfield contrast (or visibility reduction) starting from a first principles approach, explicitly modeling spatially varying refractive indices of grating and sample. By coherent propagation of the complex wavefront to the location of the analyzer grating within the optical Fresnel regime and subsequent Fourier analysis (analog to the effect of the analyzer grating), a complete model of the Talbot imaging process in a parallel beam geometry is obtained. While Yashiro et al. considered the sample to be placed in front of the Talbot interferometer, Lynch et al. considered the case of the sample placed between the modulator and analyzer grating. Both arrive at the same conclusions employing a statistic model of the sample's refractive properties on the scale below the system's spatial resolution. More specifically, the sample is characterized by its total phase shift Φ along the optical axis (its absorbing properties are treated separately), which varies in the perpendicular x, y plane (the imaging plane), i.e., Φ = Φ(x, y). Decomposing Φ(x, y) into low-and high-frequency (smooth and fine) components Φ s (x, y) and Φ f (x, y), the darkfield contrast can be attributed to the high frequency component Φ f . High frequency variations in absorption are neglected. The lower cutoff frequency of Φ f is determined by the spatial resolution of the imaging system, i.e., by its effective point spread width.

Ab initio derivation of visibility reduction
The final result given by both Yashiro and Lynch (although using different notations) relates the visibility v to the autocorrelation of the high frequency part Φ f of the sample's phase shifting with γ(ξ) being the normalized autocorrelation function (i.e., γ(0) = 1), σ 2 Φ f the variance of Φ f and ξ the correlation distance. γ(ξ) and σ 2 Φ f are to be understood as respective mean properties of Φ f on the spatial scale of the imaging system's point spread width at a given location (x, y).
The explicit x, y dependence has been omitted here for improved readability. The one dimensional autocorrelation parametrized by ξ is performed along the interferometer's sensitivity direction, i.e., perpendicular to the interferometer grating bars. This for most setups is the horizontal or x direction.
The correlation distance ξ is determined by the X-ray wavelength λ, the analyzer grating period T G2 and the distance d between sample and analyzer grating [74,33]: The correlation distance can thus be easily tuned by varying d, i.e., by moving the sample closer to or further from the analyzer grating. The system constants λ and T G2 (and the maximum possible distance d) determine the accessible order of magnitude of ξ, which typically ranges at the micrometer level for most X-ray Talbot interferometers. Although Yashiro and Lynch did not explicitly consider the case of d being larger than the G 1 -G 2 distance (i.e., the sample was assumed to be either directly in front of G 1 , or between G 1 and G 2 ), the employed Fresnel propagator formalism (cf. [46]) doesn't give reason to expect this situation to be fundamentally different.
The distinction between effective cross section σ 2 Φ f and autocorrelation γ(ξ) can also be interpreted as a representation of orthogonal sample properties. While γ(ξ) is characteristic to the structure perpendicular to the optical axis, σ 2 Φ f characterizes the sample along the optical axis. In particular, σ 2 Φ f will scale with sample thickness, while γ(ξ) depends on the sample's characteristic structure.

Relation to the linear diffusion interpretation
The geometric beam diffusion interpretation given by Wang, Bech andGrünzweig et al. (2009-2013 [63, 4, 18]) can be found as a special case of the more general result v = exp(−σ 2 Φ f (1 − γ(ξ)) given by Yashiro and Lynch when considering the limit of smooth convex (spheroid) scattering structures of diameter D ξ. The autocorrelation function γ can then be approximated by exp(− 1 2 ξ 2 (D/3) 2 ) (cf. [53] and references therein): In this limit, the linear distance dependence of the blurring width σ φ on the sample distance d as assumed in the angular diffusion interpretation is reproduced. By substitution of Eq. 11 into Eq. 6, the following identity results: The linear diffusion model becomes inaccurate as soon as γ deviates from a parabolic approximation, which either is the case for too large values of ξ/D (i.e., for too small diameters D), or when γ actually is not parabolic even for small ξ/D. The former case applies for too small scattering particles, while the latter case applies for non-smooth shaped structures. While the geometric beam diffusion interpretation was originally proposed independent of the particular imaging geometry (parallel beam versus cone beam), it's relation to the more rigorous Yashiro-Lynch theory suggests that additional geometric scaling effects as predicted by coherent optics need to be taken into account when considering laboratory instruments operating in cone beam geometry.

Magnification and Fresnel scaling in cone beam geometry
The derivations by Yashiro and Lynch, as well as the numeric simulations by Malecki et al. [36] confirming their results, were explicitly performed in a planar illumination context, which in experimental terms is commonly given only at synchrotron facilities. Laboratory setups operate at moderate distances from an X-ray point source as compared to the extent of the illuminated field of view and therefore exhibit non-negligible changes in geometric magnification for different positions along the optical axis. However, by means of the Fresnel scaling theorem, all results obtained for plane wave illumination within the Fresnel approximation can be directly transferred to a cone beam or point source illumination scenario by geometric scaling of all dimensions, i.e. [46]: with M being the geometric magnification factor defined by the ratio of source-detector distance SDD (assumed to be equivalent to the source-G 2 distance) and source-object distance SOD. The optical axis is oriented along the z-axis, x and y are the orthogonal image plane coordinates. I (SOD) is the intensity distribution at the detector (placed directly behind the analyzer grating G 2 at distance d from the sample) in the case of a source placed at distance SOD. I (∞) is the intensity pattern for the case of an infinitely distant source (i.e., the case of plane wave illumination). The Fresnel scaling theorem essentially states that classic geometric magnification based on the intercept theorem applies also to coherent wave propagation within the Fresnel regime, affecting both the plane perpendicular to the optical axis as well as propagation distances along the optical axis. While this scaling relation is taken into account in the design of laboratory type Talbot-Lau interferometers by scaling the position and period of the modulating G 1 grating appropriately with respect to the analyzer grating G 2 [64,1,12], it has, to the author's knowledge, not yet been regarded in models of the distance dependence of the visibility or darkfield contrast.
When consequently applying the above scaling relation to the correlation distance ξ (cf. Eq. 9), the effective correlation distance in terms of the actual sample dimensions becomes, for a cone beam (CB) geometry, ) .
This result can be either interpreted as a consequence of scaled sampling from the sample's phase Φ(x, y) according to the Fresnel scaling theorem, or, more intuitively, as the consequence of the geometric magnification of the sample at the location of the analyzer grating, which serves as the reference scale for all other (downscaled) planes along the optical axis (cf. the definition of M , Eq. 14). This also has implications on the angular beam diffusion model (cf. Sections 3.1 and 3.3) for the origination of darkfield contrast. By its conception of diffused incident beams, geometric scaling was assumed to arise, independent of the actual system geometry, only due to the trigonometric relation between diffusion angle σ θ and sample-detector distance d. When considering the beam diffusion model as an approximation of the Yashiro-Lynch theory in the case of large sample structures in relation to ξ, as was discussed in Section 3.3, the following modified relation is found though: i.e., σ θ ∝ d M (d) rather than ∝ d. Consequences of the specific form of ξ CB are its symmetry about d = SDD/2 and its maximum at d = SDD/2. Experimental evidence can be found e.g. in the Thesis of M. Chabior [5], who indeed finds, for a paper sample in a cone beam setup, a deviation from the quadratic distance dependence of ln(v) as expected by the linear diffusion model (cf. Section 3.1). The severe deviations found in the same work for aluminum and graphite samples can on the other hand be explained by an additional violation of the implicit D ξ assumption (cf. Section 3.3). Eq. 15 further reproduces the symmetry in the distance dependence about the midpoint between source and G 2 that was experimentally shown e.g. by Chabior [5] and Prade et al. [53]. This symmetry has previously been associated with a geometric argument given by Donath et al. [12], which was also adopted by Strobl  [60]. Figure 2 reproduces the data published by Chabior [5] (using an interferometer designed to 45keV with an analyzer grating period T G2 of 4µm, i.e., λ/T G2 ≈ 6.89 × 10 −6 ) and compares it to the Yashiro-Lynch model under consideration of Fresnel scaling (cf. Eq. 15) and the generic autocorrelation model also used in [74], with L being a characteristic correlation length and the Hurst exponent H characterizing the shape of the autocorrelation function and in particular allowing for both convex (sharp) and concave (smooth) autocorrelation profiles. Besides the correlation distance ξ, also the effective integration area of a detector pixel with respect to the sample dimensions is subject to geometric scaling. As the autocorrelation function γ(ξ CB ) represents the mean statistical properties of the sample's phase shifting properties below the imaging resolution (cf. Section 3.2), γ itself will be subject to change whenever the the imaging resolution (i.e., the pixel integration area) crosses characteristic lengths scales of the sample's structure. This effect was experimentally observed by Koenig et al. [31], who further consistently find a dependence also on the X-ray focal spot size, which likewise affects the imaging resolution.

Relation to Small Angle Scattering
Small angle scattering (SAS) techniques, or more specifically, small angle X-ray scattering (SAXS), sense a sample's differential scattering cross section (its scattering angle distribution) usually by means of illuminating it with a collimated X-ray pencil beam. A planar detector array placed behind the sample at a sufficient distance captures the scattered radiation. By application of fundamental principles from optics and Fourier analysis, the observable two-dimensional diffraction pattern is commonly shown to correspond to the Fourier transform of the auto correlation function of the sample's scattering length density (cf. e.g. the textbook by da Sivia [59]). As the latter is, on a larger scale, related to the sample's mass density, the observed pattern can be, more qualitatively, understood as the Fourier transform of the autocorrelation function of a sample's micro-or mesoscopic structure. The resolution and range of sampled correlation distances depends on the sensed angular range, with smaller scattering angles corresponding to larger correlation distances, i.e., to larger structure scales.
While SAXS commonly addresses length scales in the 10 1 to 10 2 nanometer range, the signal captured by X-ray Talbot interferometers corresponds to ultra small angle X-ray scattering (USAXS) typically related to lengths scales on the micrometer scale (cf. Eqs. 8 and 9). An important practical difference between SAXS and USAXS with respect to data analysis and interpretation emerges from the relation to the incident radiation: while the larger scattering angles in SAXS generate signals well separated from the incident radiation, USAXS analyses require explicit consideration of contributions of the original beam.
Following the convolution concept of darkfield contrast origination, Modregger et al. infer the effective point spread function (PSF) of a sample by means of actual deconvolution of the phase stepping curves (PSC) acquired with and without sample [41]. The result is directly interpreted, in analogy to classic SAXS, as differential scattering cross section and Fourier transform of the sample's autocorrelation function [40,39].
Another analogy between small angle scattering theory and darkfield imaging was drawn by Strobl [60] and Prade et al. [53], who model an incoherent superposition of scattered and unscattered fractions of the incident radiation in order to arrive at an expression equivalent to the wave optical results by Yashiro and Lynch (cf. Eq. 8) in parallel beam geometry. An experimental comparison of the latter model (and consequently also of the Yashiro-Lynch model) to SAXS is given by Gkoumas et al. [14], who investigate contributions of the structure factor to the autocorrelation function -as expected in classic scattering experiments -for dense sphere suspensions in a parallel beam synchrotron setting.
The PSF deconvolution approach is, in terms of Fourier analysis, equivalent to the evaluation of higher order Fourier components of the PSCs, of which commonly only the zeroth (mean attenuation) and the first (mean attenuation and visibility reduction) are evaluated. I.e., when directly interpreting the point spread function (or convolution kernel) as a SAXS pattern, the respective Fourier coefficients q m of the phase stepping curves with and without sample are expected to be related by an autocorrelation function γ (mξ): with m enumerating the harmonics. The prime explicitly distinguishes it from the autocorrelation function γ as defined by Yashiro and Lynch (and Strobl and Prade), although γ and γ are expected to be closely related. When directly comparing γ with the Yashiro-Lynch theory for a thin sample and m = 1, we find (assuming 1 i.e., the Yashiro-Lynch theory is qualitatively consistent with the expectations from classic scattering theory, yet differs in detail. First, SAXS theory is inherently founded on a single scattering assumption, i.e., a thin sample assumption. This gives rise to the linear vs. exponential relation with the autocorrelation function. The constant offset (1 − σ 2 Φ f ) corresponds, within the model employed by Strobl and Prade et al. [60,53], to the unscattered fraction of the incident radiation.

Lambert-Beer relation for tomographic imaging
With respect to tomographic imaging, linearity of the given contrast modality is a necessary prerequisite. I.e., the signal generated by a stack of samples must equal the sum of their isolated signals. Analogously, the signal response is required to be proportional to the thickness of a homogeneous sample. In classic X-ray computed tomography, this is provided by the well known Lambert-Beer law of intensity attenuation.
An equivalent relation is also found for the darkfield contrast modality: most directly, it can be inferred from the Gaussian convolution model given previously in Section 3. Under the assumption that multiple samples (sections of a larger sample), characterized by their effective point spread widths σ φ,i , don't interact, their combined effect will be a series of convolutions (as direct consequence of their individual effects). As convolutions of Gaussian kernels are additive in their variances, the combined signal is then given by: This argument has been employed e.g. by Wang et al. [63] (in analogy to Khelashvili et al. [30]) and Bech et al. [4] to derive the feasibility of darkfield tomography. The relation can further be generalized to non-Gaussian convolution kernels by means of the convolution theorem, as has been done by Modregger et al. [40,39]: with denoting the convolution operation and F the Fourier transformation. q is the conjugate variable to ∆φ after transformation. g(∆φ) is the effective point spread kernel found for a sample, with g i (∆φ) being individual contributions. The prime on v (q) indicates the missing normalization to v (0). The Gaussian convolution model is in fact a special case of this more general formulation. The frequency (q) dependence of v is usually not explicitly addressed in the majority of articles on the subject, as v is commonly explicitly defined as the visibility of the first harmonic of the grating period, i.e., v = v (q 1 )/v (q 0 ) (cf. Pfeiffer et al. (2008) [49]). The linear superposition of exponential arguments is, independent of the above considerations, also found by Yashiro et al. [74] and Lynch et al. [33] within their wave optical derivations of darkfield contrast origination (cf. Section 3.2).

Anisotropic Darkfield Contrast
The darkfield contrast is an oriented effect. As has been described in the previous sections, it arises from refractive effects below the spatial resolution of the imaging system and can be understood in terms of the scattering cross section and autocorrelation function of the sample's substructure. So far, the actual direction of autocorrelation was not explicitly discussed and rather implicitly determined by the orientation of the interferometer gratings. In fact, the scalar visibility contrast found in a particular experiment is to be understood as a feature of a higher dimensional quantity, which may or may not be isotropic. For samples exhibiting structural anisotropy in the plane perpendicular to the optical axis (i.e., parallel to the interferometer gratings), the darkfield signal will vary when rotating the sample (or the interferometer) about the optical axis due to the variation in characteristic length scales along the instrument's direction of sensitivity.
This anisotropy has initially been experimentally demonstrated by Wen et al. [66], Jensen et al. [22,23], Revol et al. [55], Potdevin et al. [51] and Schaff et al. [57], who considered planar samples exhibiting highly ordered fibrous structures perpendicular to the optical axis, such as carbon fibers, wood fibers, dentinal tubules or trabecular bone. Technically, these cases all address distinctive variations of the autocorrelation width of long fibers with respect to their orientation. Yashiro et al. [72] considered the case of moderately anisotropic structures within an extended sample and investigated the individual effects of autocorrelation width, scattering cross section and autocorrelation shape (by both rotating the sample and sampling multiple correlation lengths, cf. Sections 3.2 and 3.5). A systematic experimental observation of darkfield signal variation in dependence of the 3D orientation -i.e., also considering inclinations with respect to the optical axis -of highly oriented fibers was presented by Bayer et al. [3].
The theoretic modeling of these anisotropy effects has been treated phenomenologically using sinusoids reproducing the periodicity, phase and an assessment of the degree of anisotropy of the signal in dependence of a single orientation angle. This representation straight forwardly enables planar vector radiographs, augmenting 2D X-ray images with directional information on the unresolved substructure in the detection plane. Mathematically, it is equivalent to a linear approximation of the directional dependencies (as has also been pointed out by Jensen et al. [23]). Malecki et al. [34] further investigated the sinusoid approximation in numerical simulations of fibrous structures parallel to the detection plane. Current approaches extending the concept of directional imaging to tomography on volumetric samples (cf. [2,35,62,68,11,9]) derive from these planar signal models and thus implicitly presume that the autocorrelation width perpendicular to the optical axis is the principal origin of orientation dependency of the darkfield signal. A first heuristic model of darkfield anisotropy extending beyond autocorrelation effects and also accounting optical axis projected density optical axis projected density, small pixel optical axis projected density, small pixel Figure 3: Illustration of the projection of a Gaussian mass density distribution along the optical axis (left). The standard deviation of the projection will typically be larger than the distribution's standard deviation parallel to the projection plane (both marked in orange for comparison, cf. Eqs. [24][25]. When considering the case of pixels smaller than an elongated object (with regions outside that pixel shaded in gray), the projected volume contributing to an individual pixel will be confined by the pixel itself (center and right). Variations in the total projection (hatched) affecting the scattering cross section captured by that pixel are then dominated by the extent of the mass distribution along the optical axis (indicated in green). Cf. Eq. 29 vs. Eq. 30.
for variations in scattering cross section has recently been presented parallel to the present work by Felsner et al. [13].
In the following, a model of general darkfield anisotropy for arbitrarily oriented scatterers shall be derived based on the wave optical considerations given by Yashiro and Lynch, thereby implicitly accounting both for autocorrelation and scattering cross section dependencies.

Anisotropic Gaussian density model of anisotropic scatterers
In preparation to advances in X-ray darkfield tensor tomography, the anisotropy properties with respect to general sample orientations shall be reviewed based on darkfield signal origination discussed previously in Section 3.2. The anisotropy of arbitrary elongated structures shall be modeled by means of an anisotropic Gaussian mass (or electron) density distribution ρ( r) = ρ 0 e − 1 2 rT r (22) with r denoting a point in three dimensional space and being a symmetric, positive definite tensor characterized by positive eigenvalues (σ −2 1 , σ −2 2 , σ −2 3 , representing the inverse variances of ρ( r)) and a unitary rotation matrix R (with the supscript T denoting the transpose operation). R defines the orientation of the Gaussian ellipsoid object characterized by three orthogonal standard deviations σ 1 , σ 2 and σ 3 . The Gaussian mass density ρ( r) may be interpreted as description of a stochastic ensemble of smaller structures, or, more abstractly, as a generating structure to approximate the phase shifting properties along the optical axis (Eqs. [24][25] and resulting autocorrelation function (Eq. 32) of convex anisotropic structures irrespective of the precise microscopic material distribution (e.g., a cylindrical fiber). It may in this respect also be interpreted as a generalization of the Gaussian approximation of the autocorrelation properties of spherical objects (cf. Section 3.3).

Scattering cross section
The scattering cross section, as has been discussed in Section 3.2, is governed by the variance of phase front fluctuations σ 2 Φ f below the scale of the imaging system's spatial resolution. The total phase shift Φ(x, y) along the optical axis (here: z) caused by the presumed Gaussian density distribution is given by: with and Φ 0 being a constant factor resulting from the Gaussian density distribution's maximum value ρ 0 and the material's refractive index for the considered wavelength. The difference between T xx and T xx − Tzz is illustrated in Figure 3 (left). When defining the local mean variance used by Yashiro and Lynch et al. [74,33] within the discussion of darkfield origination (cf. Section 3.2) by means of a Gaussian weighting kernel accounting for the imaging point spread function's width σ PSF and assuming this point spread to be larger than the considered structure sizes σ i , i.e., σ PSF /σ i 1, we find: and can thus finally state: for objects smaller than the pixel size.
√ T zz −1 was identified here as the standard deviation along the optical axis and is therefore in more general terms described by √n Tn −1 for arbitrary orientationsn of the optical axis ( n = 1). The determinant det(T ) = (σ 1 σ 2 σ 3 ) −2 is invariant under rotations R and therefore, as Φ 0 , an invariant factor with respect to the orientation dependence of σ 2 Φ f (provided the density distribution fits into the considered integration range defined by σ PSF ).
The following conclusions can be drawn from these results on σ 2 Φ f . It is proportional to the extent √n Tn −1 of the considered structure along the optical axis, and with det(T ) −1 = σ 1 σ 2 σ 3 further proportional to its volume. In the eigenbasis of T , the relation σ 2 Φ f ∝ det(T )nTn −1 simplifies to σ 2 Φ f ∝ (σ 1 σ 2 )σ 2 3 : the variance σ 2 Φ f is directly proportional to the variance σ 2 3 of the underlying density distribution along the optical axis and scales with its cross sectional area, which is proportional to σ 1 σ 2 .
In the limit of objects extending way beyond the size of a detector pixel (e.g., long fibers), the cross sectional area of the object's projection onto a pixel becomes almost independent of its orientation, as the considered area is rather confined by the detection area of the pixel itself, which is proportional to σ 2 PSF (cf. Figure 3). In these cases, the orientation dependence of the scattering cross section affecting that pixel is thus, given the previous observations, dominated by the variance of the density distribution along the optical axis, while the cross sectional area is bounded by σ 2 PSF , as opposed to the extent of the material distribution: for objects larger than the pixel size.

Autocorrelation function
The normalized planar autocorrelation γ 2D ( ξ) corresponding to Φ(x, y) is, based on the additivity of variances for convolutions of Gaussians, given by assuming ξ to lie in the x-y plane perpendicular to the optical axis (z). When assuming that all relative rotations of interferometer and object are accounted for in R and thus in T and T Φ , the autocorrelation direction may without loss of generality be defined parallel to the x axis, i.e., ξ = (ξ, 0, 0), such that As the mass density distribution characterized by T by design already represents all structures within the optical path relevant to the considered detector pixel, no further averaging of autocorrelation properties over the pixel area is required. The autocorrelation width moreover is an immanent property of the mass density distribution, and thus in contrast to the scattering cross section not limited by the pixel extent σ PSF .
For isotropic density distributions (σ i = σ), the off-diagonals of T and T Φ will be zero. The autocorrelation then reduces to γ(ξ) = exp(− 1 4 ξ 2 σ 2 ), corresponding directly to the 3D autocorrelation function of an isotropic Gaussian density distribution without explicit averaging over the optical axis. In order to provide further intuition to the modeled Gaussian density distribution, this result may be compared to the approximate autocorrelation function for spheres of diameter D (γ(ξ) ≈ exp(− 1 2 ξ 2 (D/3) 2 ), cf. Eq. 10)). I.e., an isotropic Gaussian mass density distribution with standard deviation σ can be interpreted as modeling a sphere of diameter D ≈ 4σ.

Model of darkfield contrast anisotropy
Substituting the above results (Eqs. 29-32) on scattering cross section and autocorrelation into Eq. 8, the fringe contrast visibility is given by: which for scatterers, such as long fibers, larger than the typical integration width (e.g., pixel size) characterized by σ PSF , modifies to (34) with T characterizing the inverse variances of the considered scatterers' density distribution (Eq. 23), n the unit vector along the optical path, ξ the oriented autocorrelation distance and T Φ (Eq. 25) describing the inverse variances of the scatterers' total phase shift along the optical axis (Eq. 24), based on their density distribution characterized by T . The subscripts "sp" and "lp" indicate the "small pixel" and "large pixel" cases as compared to the extent of the considered scatterers characterized by T .
When defining, without loss of generality,n parallel to the z-axis and ξ parallel to the x-axis (assuming all rotations to be considered within T , cf. Eq. 23), and further using the first order approximation 1 − exp(−u) ≈ u for small u, then the following simplified proportionality relations for the expected orientation dependence of the darkfield contrast µ DF = − ln(v) result for the large pixel (lp) and small pixel (sp) cases respectively:

Experimental verification
With respect to applications in quantitative darkfield tomography, the Yashiro-Lynch results are explicitly reproduced in a tomography setting using a phantom containing spherules of various diameters. Furthermore, the derived anisotropy properties are reproduced using a carbon fiber phantom, with particular focus on the role of the scattering cross section in addition to the autocorrelation width. The experiments have been performed at the ID19 beamline of the European Synchrotron Radiation Facility (ESRF) using the Talbot-Interferometer by Weitkamp et al. [65]. A π-shifting phase grating G 1 of 4.8µm period and an absorbing analyzer grating G 2 of 2.4µm period were used in a monochromatic 35keV setting. Images were sampled at an effective pixel size of 61.6µm.

Feature size and positional dependence of darkfield contrast
While the dependence of the system's correlation length ξ on the distance d between sample and analyzer grating G 2 is beneficial with respect to quantitative material characterization analog to other scattering techniques, the resulting distance dependence of the darkfield signal is an undesired perturbation in the context of tomographic darkfield imaging. For centimeter scaled samples, the sample extent is typically not negligible anymore with respect to its mean distance d 0 from the analyzer grating, wherefore the distance dependence is explicitly examined using a rotating circular arrangement of differently sized spherules and compared to the theoretical expectations discussed in Sections 3.1-3.3. Figure 4 shows a sketch of the experimental setup and a summary of the respective results. The distance and sphere diameter dependent signals (Fig. 4 left) extracted from projections of the rotating sample cylinder shall be in particular compared to the explicit expressions for µ DF and the autocorrelation function γ D (ξ) for spherical particles of diameter D given by Lynch et al. [33] and Yashiro et al. [72] respectively: Constants of the material and experiment will be specified later. Beforehand, the specific methodology of the present data analysis shall be further outlined: By Eq. 11, the darkfield signal is expected to be approximately quadratic in d. The data is thus, after subtraction of the background darkfield signal caused by the sample container (cf. Fig. 4, upper left), considered in the square root domain, i.e. transformed to σ φ = √ 2µ DF = −2 ln(v). It is then fitted to the first order Taylor expansion  Figure 4: Experimental data (left) on the dependence of darkfield contrast both on the distance d between sample and analyzer grating G 2 and on the structure size. Sinusoidal distance variations for multiple samples were realized by rotating a cylindrical container of radius r ≈ 1 cm comprising nine capillary tubes arranged equidistantly along the perimeter, eight of which are filled with spherules of diameters D ranging between 0.25 µm to 80 µm (cf. sketch on the upper right). The ninth empty one serves as reference for the background signal generated by the sample container, which is approximated using a cos 2 function due to the expected 180°symmetry. Darkfield signals are evaluated for all viewing angles which allow for unobstructed views on individual capillaries (hence the uneven sampling pattern), and the background signal is subtracted prior to further processing. An example projection image is shown on the right. As the distance dependence is expected to be linearly approximable for σ φ = √ 2µ DF (cf. Eqs. 3 and 11), the signals are evaluated accordingly and approximated by a sinusoid corresponding to the variation in sample-G 2 distance. Quantitative comparisons of both the spherule size dependence and the distance dependence of the darkfield signal with the respective models given in Eqs. 37, 41 (with Eq. 38) and 43 are shown on the lower right. Deviations of µ DF (D, d 0 ) from the expected model (in particular at D = 10 µm) can be related to unaccounted variations in the volume fraction of spheres, cf. Figures 5-6. The d-dependence in cases of D < ξ (particularly pronounced for D = 0.25 µm) can be attributed to the non-vanishing autocorrelation function of the sphere arrangement rather than the spheres themselves.
of the unknown actual function σ φ (D, d) about the mean sample-G 2 -distance d 0 , with The parameters r and ω − ω 0 denote the distance and rotation phase of the considered sample capillary with respect to the rotational axis, and d 0 + ∆d(ω) describes the resulting orientation dependent sample-G 2 -distance d. Deviations from α D = 1/d 0 correspond to deviations from the assumptions of the linear diffusion model (Eqs. 5, 6 and 11), as can be easily verified by comparison of Eq. 39 with the Taylor expansion of the respective model σ φ (D, d) ∝ d σ θ (D). While the mean darkfield signal µ DF (D, d 0 ) = 1 2 σ φ (D, d 0 ) 2 contained in Eq. 39 can now be directly compared to Eq. 37 given an estimate of the involved material constants (which will be provided later), further transformations are required for a comparison of α D with the d(ω)-dependent autocorrelation function γ D (ξ(d)) (Eq. 38). By substituting the relations and into Eq. 40, the following correspondence results at d = d 0 : which can be compared to the observed values of α D determined using Eq. 39. The above expression may either be evaluated using Eq. 38, or may further be approximated using the Gaussian autocorrelation model given in Eqs. 10-11. Eq. 41 then simplifies to In the limit of large spheres, i.e., D → ∞, the following limit as expected by the linear diffusion model is approached both for Eq. 41 with Eq. 38 as well as its approximation in Eq. 43: A comparison of the experimental data ( Fig. 4 left) with Eqs. 37, 41 (using Eq. 38) and 43 is shown on the bottom right of Fig. 4. In contrast to Eq. 37, Eqs. 41-43 do have only one single free parameter, d 0 . A least squares fit yields and thus ξ 0 ≈ 1.5 µm (46) which is in excellent agreement with the actual experimental configuration of d 0 = (0.1 ± 0.01) m. Note that the autocorrelation of a given finite shape (as the spheres) inherently becomes zero for correlation distances larger than the extent (diameter) of that shape, wherefore the submicrometer sized glass spherules are not described by Eqs. 38 and 41. That they nevertheless exhibit a distance dependence (cf. Fig. 4) can rather be attributed to the fact that the dense packing of spheres of 0.25 to 1.2 µm diameter implies that the autocorrelation of the sphere arrangement, i.e., the structure factor (see also Section 3.5 and reference [14]), becomes non-negligible at the given mean correlation distance ξ 0 of 1.5 µm.
Finally, with respect to the quantitative comparison of the actual darkfield signals µ DF (D) at d 0 to the theoretic prediction as given by Eq. 37, the following constants are used: with λ being the X-ray wavelength (corresponding to 35keV), T G2 the analyzer grating period, d the sample-G 2 distance, f the volume fill factor (assuming dense sphere packing), ∆z the average sample thickness and r e ρ e,PMMA the scattering length density for PMMA. 1 ∆χ is the refractive index, whose imaginary part is, in the present case, negligible with respect to its absolute magnitude. Note the difference in notation with respect to [33] regarding µ DF which here for consistency refers to the actual darkfield signal as opposed to a thickness normalized darkfield coefficient. Thickness is here explicitly accounted for by ∆z.

Quantitative darkfield tomography
Tomographic reconstruction, i.e., the transformation of projections of an object to a volume representation of that object, generally presumes that all projections correspond to linear combinations of a static set of scalar volume elements (voxels), which are to be reconstructed. For samples rotating about an axis parallel to the interferometer's gratings, as required for tomographic imaging, the distance dependence of the darkfield contrast inevitably also leads to an apparent orientation dependence of the signal originating from isotropic scatterers away from the rotational axis. As has been shown previously, this distance dependence is not generally negligible and further not independent of the actual sample properties, which are inherently unknown prior to reconstruction. Explicit modeling of this dependence in the context of general iterative tomographic reconstruction methods would therefore be a non-trivial option. The effects of distance dependence can however be eliminated to a large degree by symmetric acquisition of projections from opposing directions, by means of performing full 360°scans: Due to the linearity -and thus additivity -of the darkfield contrast, individual volume elements (voxels) may be analyzed isolated without loss of generality. When considering an individual voxel at a distance ∆d from the rotational axis (at distance d 0 ) with respect to the optical path, the respective signals from opposing object orientations are expected to exhibit the following relations (to first order in σ φ (d)): with α characterizing the distance dependence as discussed in the previous Section. The mean signal from opposing orientations is thus with α 2 ∆d 2 ω representing the relative error with respect to µ DF (d 0 ). Given that α ≤ d −1 0 is expected (cf. Eq. 44 and Fig. 4, lower right), that error is not larger than ∆d 2 ω /d 2 0 , which can be realistically kept below 5% even for sample diameters (and thus values of ∆d) of almost up to d 0 /2. Figure 5 shows tomographic reconstructions of the data previously presented in Figure 4. The effects of darkfield distance dependence on tomographic reconstructions are illustrated by means of reconstructions from two complementary sets of projections covering an angular range of only 180°e ach.  Figure 6: Evaluation of the darkfield contrast after tomographic reconstruction (Fig. 5). In order to eliminate the volume fractions f , the ratio µ DF /µ (darkfield over absorption) is considered. For a given material, variations in µ DF /µ are now expected to arise solely due to variations in the sphere diameter D. The respective model according to Lynch (cf. Eqs. 37,47) is compared to the data, assuming d 0 = 10.7cm as found previously (Eq. 45). The absolute scale is fitted to the data. The good agreement of theory and experiment support the validity of the Yashiro-Lynch model of darkfield contrast origination with respect to the predicted structure size dependence. Further, the quantitative effects of darkfield distance dependence on short-scan tomographies are shown as well (see also Fig. 5).
A quantitative analysis of the respective gray values is given in Figure 6. By considering the absorption normalized darkfield signal µ DF /µ, effects of varying packing density are eliminated, so that variations in darkfield contrast can be expected to solely arise from differences in material (Glass vs. PMMA) and structure size (spherule diameter). The data for PMMA is compared to Eq. 37 (normalized by the absorption coefficient for PMMA at 35keV) at the correlation distance ξ 0 ≈ 1.5µm corresponding to the mean sample-G 2 distance d 0 ≈ 10.7 cm as found previously (Eq. 45), i.e., the distance between the axis of rotation and the analyzer grating G 2 : where the proportionality constant fitting the experimental data is found to be The theoretical expectation (cf. Eqs. 37 and 47) evaluates to based on the constants provided in the previous section and µ (35keV) PMMA ∆z = 0.31 cm −1 , as found in the NIST Xcom database assuming a mass density of 1.2 g cm 3 for PMMA. I.e., in addition to the distance dependence discussed previously, also the size dependence is perfectly consistent with the models given by Yashiro and Lynch, and can further be verified within tomographic reconstructions. The absolute quantitative scale found in the experiment agrees with the theoretic model within a margin of 25%.

Anisotropic darkfield contrast
In order to verify the expected orientation dependence of the darkfield contrast both with respect to effects of varying autocorrelation width and scattering cross section, a sample consisting of long carbon fibers has been imaged at a multitude of orientations in analogy to the experiment performed by Bayer et al. [3]. In contrast to the latter experiment, the sample is explicitly chosen smaller than the field of view. The carbon fibers are embedded within three fiber reinforced plastic rods of about 1cm length and 2mm diameter, and extend over the full length of the rods. The fiber sample is attached to a polygonic sample cage (made of UV resin) by means of an acrylic stand and hot glue. The cage allows to easily vary the inclination of the carbon fibers with respect to the rotational axis of a tomography setup, while the acrylic stand centers the sample in the polygonic cage. Although the support structures are made of non-scattering materials, a small contribution to the darkfield contrast is generated by their edges. Figure 7 shows selected examples. In order to keep their impact on the following analyses minimal, the signal of the sample support structures has been masked where possible (i.e., outside of the sample's silhouette).
The acquired images can be analyzed in two ways: most obviously, the average darkfield signal per detector pixel over the area of the sample silhouette may be considered, yielding a signal corresponding to fibers much longer than the pixel size (analogous to [3]). In order to instead reproduce the case of fibers fully contained within a single integrating pixel, the phase stepping curves' complex amplitudes as well as their mean transmission are averaged over the full detector area prior to the evaluation of visibility and its negative logarithm (the darkfield signal). The result is then equivalent to that of a larger integrating detector.
For the comparison of the rotation series acquired at varying fiber inclinations to the anisotropy model derived in Section 4.4, the following mass distribution model of the rotating fibers is used (cf. also Eq. 23):  Fig. 8. The rods are inclined about 65°(top row) and 36°(bottom row) with respect to the rotational axis. The grating sensitivity is parallel to the horizontal image axis. White arrows indicate examples of darkfield signals originating from the sample support structure, which has been masked outside of the sample silhouette. In the center column, the projected carbon fiber orientation is perpendicular to the sensitivity axis. The darkfield contrast is maximal in this case, and for the top row exceeds the chosen color scale ranging from 0 (black) to 0.7 (white). The center and outermost columns show the isolated effects of varying scattering cross section (per pixel) and varying autocorrelation width respectively. Quantitative results are shown in Figure 8.  In the former case, the fibers are always fully contained within the integration area. Four inclination angles of the fibers with respect to the rotational axis are considered (cf. legend), and a full rotation over 360°is performed at each inclination. Experimental data and theoretical model are shown in the left and right column respectively. Diameter/length aspect ratios σ D : σ L of 1:4 and 1:2 have been assumed for the model data (cf. Eqs. [48][49][50] shown on the right hand side for the small and large pixel case respectively in order to approximate the experimental observations on the left hand side. with σ D and σ L denoting standard deviations characterizing the effective (with respect to darkfield effects) mean diameter and length of the fiber bundle and ω and θ describing its orientation in terms of rotation and inclination as sketched in Figure 8.
For the first case considering small detector pixels (indicated by the subscript "sp"), the following is expected from the previous theoretic derivations (cf. Eq. 36): For the case of a single large detector ("lp") integrating over the full extent of all fibers, is expected in contrast (cf. Eq. 35). Figure 8 shows a respective comparison of experimental data and theoretic model. While the individual carbon fibers contained in the considered sample are expected to have an aspect ratio of about 10 3 (ca. 10 −2 m length at ca. 10 −5 m diameter), the observed signal is best reproduced with aspect ratios σ L /σ D of 4 and 2 in the small and large pixel case respectively. The general discrepancy between the extreme aspect ratio of individual fibers and the deduced aspect ratios of the fiber ensemble is expected to arise from a finite distribution width of fiber orientations. Such variations in orientation are generated whenever fibers are bent or not perfectly aligned parallel, which is especially expected among the three separate rods constituting the sample. Similarly, the reduced aspect ratio found for the case of the integrating detector (as compared to the small pixel case) might be attributed to the larger ensemble of fibers considered simultaneously in that case.
While the present data doesn't allow further microscopic analyses of the observed aspect ratios, the observed orientation dependence with respect to rotations and inclinations is in good agreement with the theoretic expectation despite the considerable number of first order approximations that have been made towards the derivation of Eqs. 49-50 (cf. Section 4). First of all, both the dependence on changes in scattering cross section and in autocorrelation width are reproduced, although the integrating pixel case appears to be more susceptible to imperfections of the sample. The scattering cross section dependence is found to be -as expected by Eq. 49 -considerably more pronounced in the case of objects exceeding the pixels' integration area. Moreover, the narrowly peaked rotation angle dependence for strongly inclined fibers is reproduced, which can be attributed to the influence of the off-diagonal term T 2 xz /T zz originating from inclinations of the anisotropic mass distribution with respect to the optical axis (cf. Figure 3 left).

Conclusions
The origination of darkfield contrast has been discussed extensively in the past from various points of view. The existing variety of approaches can be shown to be largely consistent with the wave optical derivations by Yashiro and Lynch, which on the other hand provide the crucial link to Fresnel optics, which allows for a well founded extension of the existing models to cone beam geometries based on the Fresnel scaling relation. Following the argumentation of Yashiro and Lynch, the results can furthermore be extended to anisotropic scatterers. A complete model considering both the anisotropy of the autocorrelation width as well as the scattering cross section is derived, as required for the description of arbitrarily oriented fibers (as opposed to anisotropy considerations solely within a planar sample perpendicular to the optical axis). All results are supported by experiments, with particular focus on the demands of quantitative tomography. Additional effects relevant to quantitative darkfield interpretations that have not been explicitly considered here, yet shall not be left unmentioned, include influences of beam polychromaticity [28,48,73], as well as discontinuities at material boundaries and higher order optical effects affecting the periodicity of the Talbot pattern [69].