Intragranular strain estimation in far-field scanning X-ray diffraction using a Gaussian process

Henningsson, A.; Hendriks, J.

doi:10.1107/S1600576721005112

research papers

JOURNAL OF
APPLIED
CRYSTALLOGRAPHY

ISSN: 1600-5767

Volume 54| Part 4| August 2021| Pages 1057-1070

https://doi.org/10.1107/S1600576721005112

Open

access

Intragranular strain estimation in far-field scanning X-ray diffraction using a Gaussian process

Axel Henningsson ^a ^* and Johannes Hendriks ^b

^aDivision of Solid Mechanics, Lund University, Lund, Sweden, and ^bSchool of Engineering, The University of Newcastle, Callaghan, NSW 2308, Australia
^*Correspondence e-mail: axel.henningsson@solid.lth.se

Edited by A. Borbély, Ecole National Supérieure des Mines, Saint-Etienne, France (Received 16 February 2021; accepted 13 May 2021; online 14 June 2021)

A new method for estimation of intragranular strain fields in polycrystalline materials based on scanning three-dimensional X-ray diffraction (scanning 3DXRD) data is presented and evaluated. Given an a priori known anisotropic compliance, the regression method enforces the balance of linear and angular momentum in the linear elastic strain field reconstruction. By using a Gaussian process (GP), the presented method can yield a spatial estimate of the uncertainty of the reconstructed strain field. Furthermore, constraints on spatial smoothness can be optimized with respect to measurements through hyperparameter estimation. These three features address weaknesses discussed for previously existing scanning 3DXRD reconstruction methods and, thus, offer a more robust strain field estimation. The method is twofold validated: firstly by reconstruction from synthetic diffraction data, and secondly by reconstruction of a previously studied tin (Sn) grain embedded in a polycrystalline specimen. Comparison against reconstructions achieved by a recently proposed algebraic inversion technique is also presented. It is found that the GP regression consistently produces reconstructions with lower root-mean-square errors, mean absolute errors and maximum absolute errors across all six components of strain.

Keywords: three-dimensional X-ray diffraction (3DXRD); intragranular strain; Gaussian processes; scanning X-ray diffraction.

1. Introduction

Three-dimensional X-ray diffraction (3DXRD), as pioneered by Poulsen (2004 ) and co-workers, is a nondestructive materials probe for the study of bulk polycrystalline materials. The experimental technique is typically implemented at synchrotron facilities where access to hard X-rays (≥10 keV) facilitates the study of dense materials with sample dimensions in the millimetre range. In contrast to powder diffraction, 3DXRD enables studies on a per-grain basis, which requires that the Debye–Scherrer rings consist of a set of well defined, separable single-crystal peaks. To achieve this, the beam and sample dimensions must be selected accordingly, limiting the number of grains illuminated per detector readout. By various computer-aided algorithms (cf. Lauridsen et al., 2001 ), the single-crystal diffraction peaks can be segmented and categorized on a per-grain basis, enabling the study of individual crystals within a sample. Typical quantities retrieved from such analyses are the grain average strain and average orientation (Poulsen et al., 2001 ; Oddershede et al., 2010 ). From further analysis it also possible to retrieve an approximate grain topology map (Poulsen & Schmidt, 2003 ; Poulsen & Fu, 2003 ; Markussen et al., 2004 ; Alpers et al., 2006 ).

Reducing the X-ray beam cross section to sub-grain dimensions not only allows for the study of samples with large numbers of grains but also enables the investigation of intragranular variations. This special case of 3DXRD is commonly referred to as scanning 3DXRD since, to acquire a full data set, the narrow beam must be scanned across the sample. In this setting, it is possible to measure the diffraction signal from approximate line segments across the grains, collecting information on the intragranular structure. Any inversion procedure, in pursuit of such intragranular quantities, then poses a rich tomography problem where the ray transform typically involves higher-order tensorial fields.

Recent advances in diffraction contrast tomography (Reischig & Ludwig, 2020 ) show promising results for inversion for both orientation and strain fields in three dimensions with intragranular resolution. In scanning 3DXRD where higher angular resolution on scattering vectors is achieved at the cost of diffraction peak resolution (Nervo et al., 2014 ), multiple proposals for inversion operating solely from scattering vectors have been made. Initially, Hayashi et al. (2015 , 2017 ) proposed a method for a per-voxel strain refinement to approximate intragranular strains using scanning 3DXRD data. Unfortunately, this procedure was shown to introduce bias in the reconstruction related to strain state (Hayashi et al., 2019 ; Hektor et al., 2019 ). These obstacles were later overcome by Henningsson et al. (2020 ), who proposed an inversion method that takes the tomographic nature of the problem into account. As has been pointed out by several other authors (cf. Margulies et al., 2002 ; Lionheart & Withers, 2015 ), the sampling of strain is not uniform in 3DXRD and, as a result, some additional constraints on the reconstructed field are often desirable. Henningsson et al. (2020) proposed a simple smoothing constraint to each of the strain components with success. However, the parameter selection and the physical interpretation of these constraints are not well defined.

For powder-diffraction-type data, excellent progress to overcome the weaknesses highlighted above has been made using a Gaussian process (Hendriks et al., 2020 ). In this current work, we adapt the Gaussian process (GP) framework to scanning 3DXRD and extend it to a wider class of anisotropic materials. This framework allows for the introduction of a static equilibrium constraint, which ensures that the retrieved strain reconstruction will satisfy the balance of both angular and linear momentum. The GP naturally incorporates spatial correlation in the predicted fields via a covariance function, which, together with the equilibrium prior, replaces the need for previously used smoothing constraints. Moreover, the GP produces an estimate of the uncertainty in the reconstructed strain field, as a by-product of regression. Overall, the presented regression procedure addresses several weaknesses of previous work and provides a tool for uncertainty estimation in the reconstructed strain fields.

2. Diffraction measurements

2.1. Experimental acquisition

In scanning 3DXRD, a polycrystalline specimen is placed on a sample stage associated with an attached coordinate system ( $[\hat{{\bf x}}_{\omega}]$ , $[\hat{{\bf y}}_{\omega}]$ , $[\hat{{\bf z}}_{\omega}]$ ). The sample stage commonly has several degrees of freedom, some of which are used for initial alignment and calibration and others for data collection. Since the calibration procedure is the same for all 3DXRD-type experiments, here we only describe the degrees of freedom related to data acquisition; for details on calibration see Oddershede et al. (2010), Edmiston et al. (2011 ), Borbely et al. (2014 ) and Sharma et al. (2012 ). A fixed laboratory coordinate system ( $[\hat{{\bf x}}_{\rm l}]$ , $[\hat{{\bf y}}_{\rm l}]$ , $[\hat{{\bf z}}_{\rm l}]$ ) is introduced, which is related to the sample coordinate system through a positive rotation about $[\hat{{\bf z}}_{\rm l}]$ and a translation in the $[{{ y}}_{\rm l}]$ $[{{ z}}_{\rm l}]$ plane (Fig. 1). For a given sample position (y_l, z_l), rotation in ω is performed in discrete steps of Δω. The scattered intensity in each Δω rotation interval is generally integrated during the acquisition, resulting in a series of frames for each (y_l, z_l) position. After any necessary spatial distortion corrections have been made, the raw pixelated image stacks (y_d, z_d, ω) can be segmented into separate connected regions of diffracted intensity for which centroids and average intensities can be calculated. The resulting data set is six dimensional, with each diffraction peak average intensity and detector centroid (θ, η) mapping to a sample stage setting (y_l, z_l, ω).

Figure 1
Scanning 3DXRD experimental setup. The sample coordinate system (subscript ω) is attached to the sample turntable while the laboratory (subscript l) coordinate system is fixed in relation to the sample. The sample is rotated and translated in the y_lz_l plane across the beam to record diffraction from the full volume [modified from Henningsson et al. (2020

)].

2.2. Laue equations and scattering notation

From the diffraction peak centroids (θ, η) it is possible to compute scattering vectors, G, defined in the laboratory frame as

$[{\bf G}_{\rm l} = {{2\pi}\over{\lambda}}\left[\matrix{\cos(2\theta)-1\cr -\sin(2\theta)\sin(\eta)\cr \sin(2\theta)\cos(\eta)}\right].\eqno(1)]$

Using the notation of Poulsen (2004) and considering that the Laue equations are fulfilled during diffraction, we can also express the scattering vectors as

$[{\bf G}_{\rm l} = \boldOmega{\bf U}{\bf B}{\bf G}_{hkl},\eqno(2)]$

where Ω and U are unitary square 3 × 3 rotation matrices describing, respectively, the turntable rotation around $[\hat{{\bf z}}_{\omega}]$ and the crystal unit-cell orientation with respect to the ω-coordinate system. The matrices U and B can now be uniquely defined as the polar decomposition of their inverse product, (UB)⁻¹, in which the rows contain the real-space unit-cell lattice vectors a, b and c described in the sample ω-coordinate system:

$[({\bf U}{\bf B})^{-1} = \left[\matrix{{\bf a}^{\rm T}\cr {\bf b}^{\rm T}\cr {\bf c}^{\rm T}}\right] = \left[\matrix{a_{1}&a_{2}&a_{3}\cr b_{1}&b_{2}&b_{3}\cr c_{1}&c_{2}&c_{3}}\right].\eqno(3)]$

The integer vector G_hkl = [h k l]^T holds the Miller indices.

2.3. Grain mapping

Given a measured set of scattering vectors, the procedure known as grain mapping is concerned with finding a set of uniform crystals that explain the data. In this setting, grains are represented by their average (UB)⁻¹ matrices together with their real-space centroid coordinates. To contextualize the grain-mapping procedure, a simplified schematic of the scanning 3DXRD analysis steps is presented in Fig. 2.

Figure 2
Simplified schematic of analysis steps commonly performed on scanning 3DXRD data. From raw detector data (I), the per-peak centroids (η, θ) and average intensities are retrieved (II). The scattering vectors can then be computed (III) and input into a peak-grain mapping algorithm (IV). From the produced maps, per-grain shape reconstruction can take place (V). Finally, intragranular quantities can be sought (VI).

In essence, the grain-mapping procedure results in a map between diffraction peaks and individual average grain (UB)⁻¹ matrices and centroids. The diffraction peaks associated with a single grain can be extracted from such peak–grain maps and grain shape reconstruction can be performed by tomographic methods (cf. Poulsen & Schmidt, 2003; Alpers et al., 2006), utilizing the scattered intensity associated with each diffraction peak. The peak–grain maps also enable studies on a per-grain basis, something which simplifies analysis both conceptually and computationally. Software for performing grain mapping is freely available in the ImageD11 package (Wright, 2005 ), and further details on various algorithm options can be found in the literature (Moscicki et al., 2009 ; Oddershede et al., 2010; Edmiston et al., 2011; Sharma et al., 2012; Schmidt, 2014 ). In this paper we are concerned with reconstruction of intragranular strain, and thus we focus on the final step of analysis and proceed with the assumption that all preceding quantities have been computed.

3. Measurement model

3.1. Strain revealing transformations

Henningsson et al. (2020) described the procedure to calculate strains in individual lattice planes from scanning 3DXRD measurements via the Bragg equations as first laid out by Poulsen et al. (2001) and Margulies et al. (2002). To enrich the framework, allow for consistent use of the Laue equations and clarify how the integration of strain can take place, here we adopt a different route, rewriting the Laue equations and performing a first-order Taylor series expansion. We start by recollecting that the 3 × 3 continuum deformation gradient tensor, F, should have the property that

$[{\bf v} = {\bf F}{\bf v}_{0},\eqno(4)]$

where v₀ is a vector in the reference configuration and v is the corresponding deformed vector. Applying this to a crystal reference unit cell (a₀, b₀, c₀) given in the sample ω-coordinate system and collecting the equation in matrix format, we find that

$[\left[\matrix{{\bf a}&{\bf b}&{\bf c}}\right] = {\bf F}\,\left[\matrix{{\bf a}_{0}&{\bf b}_{0}&{\bf c}_{0}}\right].\eqno(5)]$

With (3) this allows us to identify that

$[\eqalign{& ({\bf U}{\bf B})^{\rm -T} = {\bf F}({\bf U}_{0}{\bf B}_{0})^{\rm -T}\Leftrightarrow\cr & {\bf F}({\bf U}_{0}{\bf B}_{0})^{\rm -T}({\bf U}_{0}{\bf B}_{0})^{\rm T} = ({\bf U}{\bf B})^{\rm -T}({\bf U}_{0}{\bf B}_{0})^{\rm T}\Leftrightarrow\cr &{\bf F} = ({\bf U}{\bf B})^{\rm -T}({\bf U }_{0}{\bf B}_{0})^{\rm T},}\eqno(6)]$

where U₀ and B₀ define an undeformed crystal lattice. We can now relate the quantities involved in the Laue equations (1) to the strain tensor by considering that the infinitesimal strain tensor in the sample ω-coordinate system is defined as

$[\boldepsilon_{\omega} = \textstyle{{1}\over{2}}({\bf F}^{\rm T}+{\bf F}) -{\bf I},\eqno(7)]$

where I is the identity tensor. An introduction to elasticity theory is provided by Ottosen & Ristinmaa (2005 ). Insertion of (6) into (7) gives

$[\boldepsilon_{\omega} = \textstyle{{1}\over{2}}\left[({\bf U}_{0} {\bf B}_{0})({\bf U}{\bf B})^{-1}+({\bf U} {\bf B})^{\rm -T}({\bf U}_{0}{\bf B}_{0})^{\rm T}\right]- {\bf I}.\eqno(8)]$

The observable quantity in 3DXRD is the scattering vectors and a useful formulation must therefore relate ε_ω to G_ω, together with the known quantities U₀ and B₀. To achieve this we consider the strain in a single direction, introducing the unit vector $[\hat{\boldkappa}]$ into (8) as

$[\hat{\boldkappa}^{\rm T}\boldepsilon_{\omega}\hat{\boldkappa} = \textstyle{{1}\over{2}}\hat{\boldkappa}^{\rm T}\left[({\bf U}_{0} {\bf B}_{0})({\bf U}{\bf B})^{-1}+({\bf U} {\bf B})^{\rm -T}({\bf U}_{0}{\bf B}_{0})^{\rm T}\right]\hat{ \boldkappa}-1.\eqno(9)]$

The problem is now to select $[\hat{\boldkappa}]$ such that the right-hand side reduces to an observable quantity. From (2) we may define

$[{\bf G}_{\omega} = \boldOmega^{-1} {\bf G}_{\rm l} = {\bf U}{\bf B}{\bf G}_{hkl},\eqno(10)]$

and sample the strain parallel to this scattering vector as

$[\hat{\boldkappa} = {{{\bf G}_{\omega}}\over{||{\bf G}_{ \omega}||}} = {{{\bf U}{\bf B}{\bf G}_{hkl}}\over{|| {\bf G}_{\omega}||}}.\eqno(11)]$

Insertion into (9) now reduces (9) to

$[\eqalignno{& \hat{\boldkappa}^{\rm T}\boldepsilon_ {\omega}\hat{\boldkappa} \cr &\!\!=\!{\textstyle{{1}\over{2}}}{{{\bf G}_{hkl}^{\rm T}({\bf U} {\bf B})^{\rm T}}\over{||{\bf G}_{\omega}||}}\left[({\bf U}_{0} {\bf B}_{0})({\bf U}{\bf B})^{-1}\!+({\bf U} {\bf B})^{\rm -T}({\bf U}_{0}{\bf B}_{0})^{\rm T}\right]{{({\bf U}{\bf B}){\bf G}_{hkl}}\over{||{\bf G}_{\omega}|| }}-\!1 \cr &\!\!\displaystyle= {{1}\over{2{\bf G}^{\rm T}_{\omega}{\bf G}_{\omega}}} \left[{\bf G}_{hkl}^{\rm T}({\bf U}{\bf B})^{\rm T}({\bf U}_{0}{\bf B}_{0}){\bf G}_{hkl}+{\bf G}_{hkl}^{\rm T}({\bf U}_{0}{\bf B}_{0})^{\rm T}({\bf U}{\bf B}) {\bf G}_{hkl}\right]-\!1 \cr &\!\!=\displaystyle{{1}\over{2{\bf G}^{\rm T}_{\omega}{\bf G}_{\omega}}} \left[{\bf G}^{\rm T}_{\omega}{\bf G}^{(0)}_{\omega}+({\bf G}^{(0)}_{\omega})^{\rm T}{\bf G}_{\omega}\right]-1 = {{{\bf G}_ {\omega}^{\rm T}{\bf G}^{(0)}_{\omega}}\over{{\bf G}_{\omega}^{\rm T} {\bf G}_{\omega}}}-1,&(12)}]$

where

$[{\bf G}^{(0)}_{\omega} = \boldOmega^{-1}{\bf G}^{(0)}_{\rm l} = {\bf U}_{0}{\bf B}_{0}{\bf G}_{hkl}.\eqno(13)]$

This selection of unit vector $[\hat{{ \boldkappa}}]$ not only guarantees that ε_ω is the only unknown in (12) but further spreads the sampling of strain to all directions defined by the measured set of scattering vectors G_ω. For high X-ray energies, although not uniform, this spread is typically good (Lauridsen et al., 2001), explaining why, in general, strain reconstruction is possible in 3DXRD.

3.2. Tensorial ray transform

So far we have worked with equations (2)–(12) as if scattering occurred from a single point. This is typically the approximation made in 3DXRD when only grain average properties are required. For scanning 3DXRD, when pursuing intragranular quantities, we must consider that scattering takes place from grain sub-regions, illuminated by the narrow X-ray beam. In fact, if the scattered intensity is the same from all points within the grain, the scattering vectors known to us from the experiment are average quantities over regions, $[{\cal R}]$ , within the grain such that

$[\langle{\bf G}_{\omega}\rangle = {{1}\over{V}}\int\limits_{{\cal R}} {\bf G}_{\omega}\,{\rm d}v = {{1}\over{V}}\int\limits_{{\cal R}}{\bf U}{\bf B}{\bf G}_{hkl}\,{\rm d}v,\eqno(14)]$

where V is the total volume of $[{\cal R}]$ , dv is the differential on $[{\cal R}]$ and 〈·〉 indicates volume average. We run now the risk of invalidating our previous result (12) since the local scattering vectors G_ω = G_ω(x_ω, y_ω, z_ω) are unknown in scanning 3DXRD. To maintain a useful expression we must further transform (12) into an equation in 〈G_ω〉 rather than G_ω. However, since the strain is nonlinear in G_ω, direct volume integration of (12) is not possible. Fortunately though, we may obtain an approximation by Taylor expansion of (12) at $[{\bf G}_{\omega} = {\bf G}_{\omega}^{(0)}]$ to first order:

$[\hat{{\boldkappa}}^{\rm T}\boldepsilon_{\omega}\hat{\boldkappa}\simeq 1-{{{\bf G}_{\omega}^{\rm T}{\bf G}^{(0)}_{ \omega}}\over{({\bf G}^{(0)}_{\omega})^{\rm T}{\bf G}^{(0)}_{\omega}}}.\eqno(15)]$

By selecting a uniform reference configuration in space, integration of (15) now gives, with (14), that

$[y = {{1}\over{V}}\int\limits_{{\cal R}}\hat{\boldkappa}^{\rm T}\boldepsilon_{\omega}\hat{\boldkappa}\,{\rm d}v\simeq 1-{{ \langle{\bf G}_{\omega}\rangle^{\rm T}{\bf G}^{(0)}_{\omega}}\over{({\bf G}^{(0)}_{\omega})^{\rm T}{\bf G}^{(0)}_{\omega}}},\eqno(16)]$

where we introduce the scalar average strain measure $[y = y(\hat{\boldkappa})]$ .

Finally, in any inversion scheme where ε_ω constitute the free variables, we must be able to execute the forward model that is the integral of (16). For this purpose the direction of strain, $[\hat{\boldkappa}]$ , must be approximated. Using the already introduced assumption that G_ω varies weakly over $[{\cal R}]$ we can write

$[\hat{\boldkappa}\simeq{{\langle{\bf G}_{\omega}\rangle}\over{ ||\langle{\bf G}_{\omega}\rangle||}}.\eqno(17)]$

We note that, equally, the approximation $[\hat{\boldkappa}\simeq{\bf G}^{(0)}_{\omega}/||{\bf G} ^{(0)}_{\omega}||]$ could have been made.

In conclusion, (16) and (17) relate the measured average scattering vectors, 〈G_ω〉, to the underlying strain field, ε_ω(x_ω, y_ω, z_ω), with the strain tensor being the only involved unknown quantity.

The approximations made in (16) and (17) will give rise to an error in the integrated strain value y. The magnitude of this error will strongly depend on the spatial distribution of intragranular strain and orientation. To demonstrate that the approximations made in (16) and (17) are accurate for small strains and moderate mosaic spreads, we provide an extended analysis of this error in Appendix A. This discussion also highlights why, and when, it is possible to reconstruct intragranular strain independently of intragranular orientation in scanning 3DXRD.

3.3. Estimated uncertainty

To finalize the measurement model we introduce an additive Gaussian error e into (16), representing measurement uncertainty. Furthermore, to simplify both computation and further analytical derivations we approximate the volume integral over $[{\cal R}]$ by a corresponding line integral going through the geometrical centre of this region. Finally, we have the measurement model

$[y = {{1}\over{L}}\int\limits_{{\cal L}}\hat{\boldkappa}^{\rm T}\boldepsilon_{\omega}\hat{\boldkappa}\,{\rm d}l+e,\eqno(18)]$

where L is the length of the line segment $[{\cal L}]$ , dl is the differential on $[{\cal L}]$ and e is the additive normally distributed noise:

$[e\sim {\cal N}({\bb E}[e],{\bb C}[e,e]),\eqno(19)]$

with expectation value $[{\bb E}[e]]$ and covariance $[{\bb C}[e,e]]$ .

The measurement noise is assumed to be zero mean ( $[{\bb E}[e_{i}] = 0]$ ) and independent ( $[{\bb C}[e_{i},e_{j}] = 0]$ ) with the variance selected in accordance with previous work (Borbely et al., 2014; Henningsson et al., 2020),

$[{\bb C}[e_{i},e_{i}] = \left({{\partial y}\over{\partial r}}\right) ^{-2},\eqno(20)]$

where r = r(θ) is the radial detector coordinate and the indices i and j indicate unique measurements. Other estimations of $[{\bb C}[e_{i},e_{i}]]$ are possible. Importantly, though, the variance should depend on the scattering angle 2θ, as, for a 2D detector with uniform pixel size, the measurement uncertainty increases with decreasing scattering angle.

4. Regression procedure

Equation (18) is a ray transform that contains information on the average directional strain for a region within the grain. The problem to reconstruct the full strain tensor field from a series of such measurements is therefore tomographic in nature, and the measurements y are highly spatially entangled as the regions $[{\cal L}]$ will intersect in general. A collection of N measurements,

$[{\bf y} = \left[\matrix{y_{1}&y_{2}&...&y_{j}&...&y_{N}}\right]^{\rm T},\eqno(21)]$

could represent the second member of a linear equation system where (18) is used to form a system matrix and a vector of unknown strains defined on some finite basis. This has been described by Henningsson et al. (2020) for a voxel basis, using a weighted least-squares (WLSQ) approach to retrieve the strain field. As we will discuss in Section 4.1, in this work we adapt these ideas to a Gaussian process framework, not solving for a deterministic strain field but instead calculating the probability distribution of strain at each spatial coordinate, revealing a distribution over strain tensor functions.

Before proceeding any further, it is useful to introduce a vector notation along with some geometrical quantities related to the integration path $[{\cal L}]$ (Fig. 3).

Figure 3
A single crystal under elastic deformation illuminated by an X-ray beam. Scattering takes place along the illuminated region $[{\cal L} = {\cal L}_{1}+{\cal L}_{2}]$ .

Since ε^T = ε we can uniquely represent the strain tensor field in sparser format by introducing the column vector

$[\bar{\boldepsilon}({\bf x}) = \left[\matrix{\epsilon_{xx}({\bf x})&\epsilon_{yy}({\bf x})&\epsilon_{zz}({\bf x})& \epsilon_{xy}({\bf x})&\epsilon_{xz}({\bf x})&\epsilon_{yz}({\bf x})}\right]^{\rm T}.\eqno(22)]$

To represent the tensor product $[\hat{\boldkappa}^{\rm T}\boldepsilon_{\omega}\hat{\boldkappa}]$ involved in (18) using $[\bar{\boldepsilon}]$ we seek the corresponding vector $[\bar{\boldkappa}]$ such that the equality $[\bar{\boldkappa}^{\rm T}\bar{\boldepsilon} = \hat{\boldkappa}^{\rm T}\boldepsilon_{\omega}\hat{\boldkappa}]$ holds true. We find by expansion that

$[\bar{\boldkappa} = \left[\matrix{\kappa_{x}^{2}&\kappa_{y}^{2}&\kappa_ {z}^{2}&2\kappa_{x}\kappa_{y}&2\kappa_{x}\kappa_{z}&2\kappa_{y}\kappa_{z}}\right]^{\rm T}.\eqno(23)]$

Next, denoting the intersection points between the X-ray beam and the grain boundary by p₀, p₁,…, p_M and letting the Euclidean length of these illuminated regions be labelled L_i = ||p_i − p_i+1||₂, we find, for measurement number j, that

$[y_{j} = e_{j}+\sum^{i = M-1}_{i = 0}{{1}\over{L_{i}}}\int\limits^{L_{i}}_{0}\bar{{\boldkappa}}^{\rm T}\bar{\boldepsilon}({\bf p}_{i}+\hat{{\bf n }}s)\,{\rm d}s = e_{j}+{\cal M}_{j}\bar{\boldepsilon},\eqno(24)]$

where the symbol $[{\cal M}_{j}]$ is shorthand for the integral operator corresponding to measurement number j, s is a scalar, $[\hat{{\bf n}}]$ is a unit vector along the X-ray beam and $[\bar{\boldepsilon} = \bar{\boldepsilon}({\bf p}_{i}+ \hat{{\bf n}}s)]$ is a function over a compact support in the grain volume. Considering the full measurement set y defined in (21), we introduce a compact notation,

$[{\bf y} = { {\cal M}\!\!\!\!\!\!\!{\cal M}}\bar{\boldepsilon}+{\bf e},\eqno(25)]$

where $[{\cal M}\!\!\!\!\!\!\!{\cal M}]$ and e are column vectors formed in analogy with (21).

4.1. Gaussian process regression

A Gaussian process is any stochastic process in which all subsets of a generated stochastic sequence of measurements form multivariate normal distributions (Rasmussen, 2003 ). The regression procedure associated with a Gaussian process, known as Gaussian process regression, can be described in terms of basic statistical theorems and quantities. The central idea is to exploit the fact that linear operators acting on normally distributed variables form again normal distributions. The goal is to arrive at the distribution of the Gaussian process that, for some spatial function f(x), describes the probability of finding a value f at coordinate x together with the covariance of f(x) with other spatial locations f(x′).

In the scanning 3DXRD case, we consider the measurement series, y, generated by some underlying strain tensor field, $[\bar{\boldepsilon}({\bf x})]$ , and seek to calculate at each spatial coordinate, x, the probability distribution $[p[{\bar{\boldepsilon}}({\bf x})|{\bf y}]]$ , i.e. the probability of finding a specified strain tensor $[\bar{\boldepsilon}]$ at x given the measurements y. As we will show, if we assume a Gaussian process prior and Gaussian noise, this probability distribution is multivariate normal, and the covariance of strain at any two points, $[{\bb C}[\bar{\boldepsilon} = \bar{\boldepsilon}( {\bf x}),\bar{\boldepsilon}^{\prime} = \bar{\boldepsilon}({\bf x}^{\prime})]]$ , together with the strain expectation value, $[{\bb E}[\bar{\boldepsilon}({\bf x})]]$ , will be revealed by the regression.

If it is assumed that $[\bar{\boldepsilon}({\bf x})]$ is normally distributed,

$[\bar{\boldepsilon}({\bf x})\sim{\cal N}({\bb E} [\bar{\boldepsilon}],{\bb C}[\bar{\boldepsilon},\bar{ \boldepsilon}^{\prime}]),\eqno(26)]$

it follows directly that y is multivariate normal,

$[{\bf y}\sim{\cal N}({\bb E}[{\bf y}],{\bb C} [{\bf y},{\bf y}]),\eqno(27)]$

since it is a linear combination of the independent normal distributions $[\bar{\boldepsilon}({\bf x})]$ and e. Considering then the joint distribution of $[\bar{\boldepsilon}({\bf x})]$ and y we can calculate

$[\eqalignno{& \left[\matrix{\bar{\boldepsilon}\cr {\bf y}}\right]\cr & \sim{\cal N}\left( \left[\matrix{{\bf I }\cr {\cal M}\!\!\!\!\!\!\!{\cal M}}\right]{\bb E}[\bar{\boldepsilon}], \left[\matrix{{\bb C}[\bar{\boldepsilon},\bar{\boldepsilon}^{\prime}]&{\bb C}[\bar{\boldepsilon},\bar{{ \boldepsilon}}^{\prime}]{\cal M}\!\!\!\!\!\!\!{\cal M}^{\rm T}\cr {\cal M}\!\!\!\!\!\!\!{\cal M}{\bb C}[\bar{\boldepsilon},\bar{ \boldepsilon}^{\prime}]&{\cal M}\!\!\!\!\!\!\!{\cal M}{\bb C}[\bar{ \boldepsilon},\bar{\boldepsilon}^{\prime}]{\cal M}\!\!\!\!\!\!\!{\cal M}^{\rm T}+{\bb C}[{\bf e},{\bf e}]}\right]\right),\cr &&(28)}]$

where I is an identity operator and we use the fact that y is a linear transformation of two normally distributed variables $[\bar{\boldepsilon}({\bf x})]$ and e. The joint probability of (28) now gives us the sought distribution, $[p[\bar{\boldepsilon}({\bf x})|{\bf y}]]$ , which is again normal. Its variance and expectation value can be found by writing out (28) in analytical exponent form, with fixed y, and completing the exponent square. The closed-form solution can be obtained as

$[\eqalign{{\bb E}[\bar{\boldepsilon}|{\bf y}] & = {\bb E}[\bar{\boldepsilon}]+{\bb C}[\bar{\boldepsilon},\bar{\boldepsilon}^{\prime}]{\cal M}\!\!\!\!\!\!\!{\cal M}^{\rm T}\big({\cal M}\!\!\!\!\!\!\!{\cal M}{\bb C}[\bar{\boldepsilon},\bar{ \boldepsilon}^{\prime}]{\cal M}\!\!\!\!\!\!\!{\cal M}^{\rm T}\cr & \quad+{\bb C} [{\bf e},{\bf e}]\big)^{-1}\left({\bf y} -{\bb E} [{\bf y}]\right),\cr {\bb C}[\bar{\boldepsilon},\bar{\boldepsilon }^{\prime}|{\bf y}] & = {\bb C}[\bar{\boldepsilon},\bar{ \boldepsilon}^{\prime}]-{\bb C}[\bar{\boldepsilon},\bar{ \boldepsilon}^{\prime}] {\cal M}\!\!\!\!\!\!\!{\cal M}^{\rm T} \big( {\cal M}\!\!\!\!\!\!\!{\cal M}{\bb C}[\bar{\boldepsilon},\bar{ \boldepsilon}^{\prime}]{\cal M}\!\!\!\!\!\!\!{\cal M}^{\rm T} \cr & \quad+{\bb C} [{\bf e},{\bf e}]\big)^{-1}{\cal M}\!\!\!\!\!\!\!{\cal M}{\bb C} [\bar{\boldepsilon},\bar{\boldepsilon}^{\prime}].}\eqno(29)]$

Before any approximate or analytical solutions to the involved transformations of $[{\bb C}[\bar{\boldepsilon},\bar{\boldepsilon}^{\prime}]]$ by $[{\cal M}\!\!\!\!\!\!\!{\cal M}]$ can be given, it remains first to specify the prior distribution of $[\bar{\boldepsilon}({\bf x})]$ .

4.2. Equilibrium prior

Since the closed-form solution of (29) requires only that $[\bar{\boldepsilon}({\bf x})]$ is normal, we are free to incorporate prior knowledge on $[\bar{\boldepsilon}({\bf x})]$ by making a parametrization of $[\bar{\boldepsilon}({\bf x})]$ as linear transformations of some other underlying normal distributions. Since $[\bar{\boldepsilon}({\bf x})]$ represents a linear elastic strain field and the scanning 3DXRD experiment is assumed to take place on a sample at rest, we expect that the accompanying stress field $[\bar{{\boldsigma}}]$ will be in static equilibrium. This can be expressed as a linear map

$[\bar{\boldepsilon}({\bf x}) = {\bf H}\bar{{ \boldsigma}}({\bf x}),\eqno(30)]$

where H is an anisotropic compliance matrix that is orientation dependent, $[{\bf H} = {\bf H}({\bf U})\simeq{\bf H}({\bf U}_{0})]$ . The set of analytical functions $[\bar{{\boldsigma}}({\bf x})]$ that satisfy balance of angular and linear momentum are known as the Beltrami stress functions. These may be described as a linear map

$[\bar{{\boldsigma}}({\bf x}) = {{\cal B}\!\!\!\!{\cal B}}\bar{ {\boldPhi}}({\bf x}),\eqno(31)]$

where $[\bar{{\boldPhi}}({\bf x})]$ is a column vector holding six Beltrami stress functions, which are required to be twice differentiable, and

$[{\cal B}\!\!\!\!{\cal B} \!=\! \left[\matrix{\!\!0&\displaystyle{{\partial^{2}} \over {\partial z^{2}}}\!\!& \displaystyle{{\partial^{2}} \over {\partial y^{2}}}\!\!&0\!\!&0\!\!&\displaystyle-2{{\partial^{2}} \over {\partial y \partial z}}\!\!\cr \displaystyle{{\partial^{2}} \over {\partial z^{2}}}\!\!&0\!\!&\displaystyle{{\partial^{2}} \over {\partial x^{2}}}\!\!&0\!\!&\displaystyle-2 {{\partial^{2}} \over {\partial x\partial y}}\!\!&0\!\!\cr \!\!\displaystyle{{\partial^{2}} \over {\partial y^{2}}}\!\!&\displaystyle{{\partial^{2}} \over {\partial x^{2}}}&0&\displaystyle-2 {{\partial^{2}} \over {\partial x\partial y}}\!\!&0\!\!&0\!\!\cr \!\!0\!\!&0\!\!&\displaystyle-{{\partial^{2}} \over {\partial x\partial y}}\!\!&\displaystyle-{{\partial^{2}} \over {\partial z ^{2}}}\!\!&\displaystyle{{\partial^{2}} \over {\partial y\partial z}}\!\!&\displaystyle{{\partial^{2}} \over {\partial x \partial z}}\!\!\cr\!\! \displaystyle-{{\partial^{2}} \over {\partial y\partial z}}\!\!&0\!\!&0\!\!&\displaystyle{{\partial^{2}} \over {\partial x \partial z}}\!\!&\displaystyle{{\partial^{2}} \over {\partial x\partial y}}\!\!&\displaystyle-{{\partial^{2}} \over { \partial x^{2}}}\!\!\cr\!\! 0\!\!&\displaystyle-{{\partial^{2}} \over {\partial x\partial z}}\!\!&0\!\!&\displaystyle{{\partial^{2}} \over {\partial y \partial z}}\!\!&\displaystyle-{{\partial^{2}} \over {\partial y^{2}}}\!\!&\displaystyle{{\partial^{2}} \over {\partial x \partial y}}\!\!}\right].\eqno(32)]$

We have then

$[\bar{\boldepsilon}({\bf x}) = {\bf H}{\cal B}\!\!\!\!{\cal B}\bar{{\boldPhi}}({\bf x}),\eqno(33)]$

and must now make an assumption on the distribution of $[\bar{\boldPhi}({\bf x})]$ . Without any further prior knowledge we select a zero-mean normal distribution with

$[{\bb E}[\bar{\boldPhi}] = \left[\matrix{0 \cr 0\cr 0\cr 0\cr 0\cr 0}\right],\quad{\bb C}[\bar{\boldPhi},\bar{\boldPhi} ^{\prime}] = \left[\matrix{k_{1}&0&0&0&0&0\cr 0&k_{2}&0&0&0&0\cr 0&0&k_{3}&0&0&0\cr 0&0&0&k_{4}&0&0\cr 0&0&0&0&k_{5}&0\cr 0&0&0&0&0&k_{6}}\right],\eqno(34)]$

where the covariance functions $[k_{i} = k_{i}({\bf x},{\bf x}^{\prime})]$ describe the spatial correlation of the field. In this work, we have used the stationary squared-exponential covariance function,

$[k_{i}({\bf x},{\bf x}^{\prime}) = \sigma_{i}^{2}\exp\left({ {-{\bf r}^{\rm T}{\bf r}}\over{2{\bf l}_{i}^{\rm T}{\bf l}_{i}}} \right),\quad{\bf r} = {\bf x}-{\bf x}^{\prime},\quad {\bf l}_{i} = \left[\matrix{l_{ix}&l_{iy}&l_{iz}}\right]^{\rm T},\eqno(35)]$

introducing a smoothness assumption into the strain field reconstruction. The unknown hyperparameters defined by $[{\bf l}_{i}]$ and σ_i are thus in total 6 × 4 = 24 in our case. These variables will be estimated through an initial optimization process known as hyperparameter optimization; we will return to how this is done later. First we highlight that the zero-mean prior assumption on the Beltrami stress functions, $[\bar{\boldPhi}({\bf x})]$ , does not imply that the posterior distribution of strain, $[\bar{\boldepsilon}({\bf x})]$ , will be zero mean. This is realized upon examination of equation (29), which shows that a prior mean of $[{\bb E}[\bar{\boldepsilon}] = 0]$ does not imply that the conditional posterior $[{\bb E}[\bar{\boldepsilon}|{\bf y}]]$ will be zero. Other selections for the prior mean are possible; however, when such additional prior information is unknown, a zero-mean selection is preferable for simplicity.

In total, these selections impose that (i) the strain field is in a point-wise static equilibrium and (ii) the strain field has a local spatial correlation to neighbouring points. The resulting prior distribution of strain is

$[\bar{\boldepsilon}\sim{\cal N}({\bf 0},{\bf H }{\cal B}\!\!\!\!{\cal B}{\bb C}[\bar{\boldPhi},\bar{\boldPhi}^{\prime}]{\bf {\cal B}}^{\rm T}{\bf H}^{\rm T}).\eqno(36)]$

4.3. Equilibrium posterior distribution

With the prior information of equilibrium and spatial correlation now encoded into the strain field we can insert

$[{\bb C}[\bar{\boldepsilon},\bar{\boldepsilon}^{\prime}] = {\bf H}{\cal B}\!\!\!\!{\cal B}{\bb C}[\bar{\boldPhi},\bar{ \boldPhi}^{\prime}]{\cal B}\!\!\!\!{\cal B}^{\rm T}{\bf H}^{\rm T}\eqno(37)]$

into equation (29) to arrive at a final expression in which only the hyperparameters remain to be estimated. The covariance between measurements takes on the form

$[{\bb C}[{\bf y},{\bf y}] = {\cal M}\!\!\!\!\!\!\!{\cal M}{\bf H}{\cal B}\!\!\!\!{\cal B}{\bb C}[\bar{\boldPhi},\bar{\boldPhi}^{\prime}]{\cal B}\!\!\!\!{\cal B}^{\rm T}{\bf H}^{\rm T}{\cal M}\!\!\!\!\!\!\!{\cal M}^{\rm T},\eqno(38)]$

which involves, through the mappings $[{\cal M}\!\!\!\!\!\!\!{\cal M}]$ , a double integral over the two times partially differentiated squared exponential in (35). The solution to this double line integral is intractable, although some work has been done to show that for l_x = l_y = l_z it can be analytically reduced to a single integral (Hendriks, Gregg et al., 2019 ). However, the numerical integration remains too computationally costly for practical use. This motivates the use of an approximation scheme on a reduced basis for which closed-form solutions to all involved quantities of (29) are again recovered (Jidling et al., 2018 ).

4.4. Finite basis approximations

Decomposing (35) onto a Fourier basis,

$[\eqalignno{\varphi_{ik}({\bf x}) & = {{1}\over{L_{x}L_{y}L_{z}}}\sin[\lambda_{xik}(x+L_ {x})]\sin[\lambda_{yik}(y+L_{y})]\cr & \quad\times\sin[\lambda_{zik}(z+L_{z})],&(39)}]$

where the scalars λ and L are the frequencies and phases, respectively, we find that

$[k_{i}({\bf x},{\bf x}^{\prime})\simeq\textstyle\sum\limits^{k = m}_{k = 1}\varphi_{ ik}({\bf x})s_{ik}\varphi_{ik}({\bf x}^{\prime}) = {\boldvarphi}^{\rm T}_{i}{\bf s}_{i}{\boldvarphi}^{\prime}_{i}.\eqno(40)]$

s_i is a diagonal matrix of basis coefficients, s_ik, which are the spectral densities of (35). Specifically it is possible to show (Solin & Särkkä, 2020 ) that the kth spectral density is

$[s_{ik} = \sigma^{2}_{i}(2\pi)^{{{2}/{3}}}l_{ix}l_{iy}l_{iz}\exp\left[\textstyle-{{ 1} \over {2}}(l^{2}_{ix}\lambda^{2}_{xik}+l^{2}_{iy}\lambda^{2}_{yik}+l^{2}_{iz} \lambda^{2}_{zik})\right].\eqno(41)]$

With the vector notation

$[{\boldphi} = \left[\matrix{\boldvarphi_{1}\cr \boldvarphi_{2}\cr \boldvarphi_{3}\cr \boldvarphi_{4}\cr \boldvarphi_{5}\cr \boldvarphi_{6}}\right],\quad{\bf S} = \left[\matrix{ {\bf s}_{1}&{\bf 0}&{\bf 0}&{\bf 0}&{\bf 0} &{\bf 0}\cr {\bf 0}&{\bf s}_{2}&{\bf 0}&{\bf 0}&{\bf 0} &{\bf 0}\cr {\bf 0}&{\bf 0}&{\bf s}_{3}&{\bf 0}&{\bf 0} &{\bf 0}\cr {\bf 0}&{\bf 0}&{\bf 0}&{\bf s}_{4}&{\bf 0} &{\bf 0}\cr {\bf 0}&{\bf 0}&{\bf 0}&{\bf 0}&{\bf s}_{5} &{\bf 0}\cr {\bf 0}&{\bf 0}&{\bf 0}&{\bf 0}&{\bf 0}& {\bf {s}_{6}}}\right],\eqno(42)]$

where 0 is a matrix of zeros, we find the approximate covariance

$[{\bb C}[\bar{\boldPhi},\bar{\boldPhi}^{\prime}] = {\boldphi}^{\rm T}{\bf S}{\boldphi}^{\prime}.\eqno(43)]$

Insertion of (43) into (37) now yields

$[{\bb C}[\bar{\boldepsilon},\bar{\boldepsilon}^{\prime}] = {\bf H}{\cal B}\!\!\!\!{\cal B}{\boldphi}^{\rm T}{\bf S} {\boldphi}^{\prime}{\cal B}\!\!\!\!{\cal B}^{\rm T}{\bf H}^{\rm T}.\eqno(44)]$

Introducing the quantities

$[\boldphi_{\epsilon} = {\bf H}{\cal B}\!\!\!\!{\cal B}{\boldphi}^{\rm T},\quad\boldphi_{y} = {\cal M}\!\!\!\!\!\!\!{\cal M}{\boldphi }_{\epsilon},\eqno(45)]$

we finally arrive at the approximate posterior mean and covariance of strain using (29):

$[\eqalign{{\bb E}[\bar{\boldepsilon}|{\bf y}]& = {\bb E}[\bar{\boldepsilon}]+\boldphi_{\epsilon} {\bf S}\boldphi_{y}^{\rm T}\left(\boldphi_{y} {\bf S}\boldphi_{y}^{\rm T}+{\bb C}[{\bf e},{\bf e}]\right)^{-1}\left({\bf y}-{\bb E}[{\bf y}]\right),\cr {\bb C}[\bar{\boldepsilon},\bar{{\boldepsilon }}|{\bf y}] & = \boldphi_{\epsilon}{\bf S}\boldphi _{\epsilon}^{\rm T}-\boldphi_{\epsilon}{\bf S}\boldphi_{y }^{\rm T}\left(\boldphi_{y}{\bf S}\boldphi_{y}^{\rm T}+ {\bb C}[{\bf e},{\bf e}]\right)^{-1}\boldphi_{y} {\bf S}\boldphi^{\rm T}_{\epsilon}.}\eqno(46)]$

The computational complexity can be further reduced by algebraically rearranging this equation to avoid forming the covariance matrices (Rasmussen, 2003), resulting in

$[\eqalign{{\bb E}[\bar{\boldepsilon}|{\bf y}] & = {\bb E}[\bar{\boldepsilon}]+\boldphi_{\epsilon}\big(\boldphi_{y}^{\rm T}{\bb C}[{\bf e},{\bf e}]^{-1} \boldphi_{y}\cr & \quad+{\bf S}^{-1}\big)^{-1} \boldphi_{y}^{\rm T }{\bb C}[{\bf e},{\bf e}]^{-1}\left({\bf y}-{\bb E}[{\bf y}]\right),\cr {\bb C}[\bar{\boldepsilon},\bar{{\boldepsilon }}|{\bf y}] & = \boldphi_{\epsilon}\left(\boldphi_{y}^{\rm T}{\bb C}[{\bf e},{\bf e}]^{-1}\boldphi_{y}+ {\bf S}^{-1}\right)^{-1}\boldphi_{\epsilon}^{\rm T}.}\eqno(47)]$

Here, the inverses S⁻¹ and $[{\bb C}[{\bf e},{\bf e}]^{-1}]$ can be trivially computed, as the matrices are diagonal. For m < N, this reduces the computational complexity to $[{\cal O}(Nm^{2})]$ from $[{\cal O}(N^{3})]$ required for the inverse in (29) and (46). A numerically stable and efficient algorithm for solving these equations using the QR decomposition is given by Hendriks, Wensrich et al. (2019 ), together with analytical expressions for the various integral mappings $[{\cal M}\!\!\!\!\!\!\!{\cal M}]$ . We note here that, although the introduced Fourier basis in (39) is defined over all space, the support of the reconstructed field in (47) is for all practical purposes that of the grain volume. This follows from the fact that the mappings executed through $[{\cal M}\!\!\!\!\!\!\!{\cal M}]$ are only performed over the grain, as indicated in (24), and requires that the period of the lowest frequency basis included is larger than the grain volume.

As m → ∞ the approximate solution (47) approaches the exact solution (29) (Solin & Särkkä, 2020). In practice, however, we must select a finite m, leading to (35) being used in approximate form. To direct the selection of frequencies λ_xik, λ_yik and λ_zik in (40) use can be made of (41). In this work, we have selected the basis frequencies on an equidistant grid in (λ_xik, λ_yik, λ_zik) space such that the spectral densities were above a minimum threshold, i.e. we aim to achieve a desired coverage of the spectral density function. Specifically, we select

$[\eqalign{&\lambda_{xik} = \Delta\lambda_{xki}g_{xki},\quad L_{x } = {{\pi}\over{2\Delta\lambda_{xki}}},\quad\Delta\lambda_{xki} = {{\nu}\over{l_{ix }R}},\cr & \lambda_{yik} = \Delta\lambda_{yki}g_{yki},\quad L_{y} = {{\pi}\over {2\Delta\lambda_{yki}}},\quad\Delta\lambda_{yki} = {{\nu}\over{l_{iy}R}},\cr & \lambda_{zik} = \Delta\lambda_{zki}g_{zki},\quad L_{x} ={{\pi} \over{2\Delta\lambda_{zki}}},\quad\Delta\lambda_{zki} = {{\nu}\over{l_{iz}R}},\cr & g_{xki}^{2}+g_{yki}^{2}+g_{zki}^{2}\leq R^{2},}\eqno(48)]$

where (g_xki, g_yki, g_zki) are positive integers such that (Δλ_xkig_xki, Δλ_ykig_yki, Δλ_zkig_zki) defines equidistant grid points excluding the origin, and ν controls the desired coverage of the spectral density.

To see how ν controls this coverage, we use equation (48) to write the spectral density in (41) as a function of ν, giving

$[\eqalignno{ s_{ik}&\displaystyle = \sigma^{2}_{i}(2\pi)^{{{2} /{3}}}l_{ix}l_{iy}l_{iz}\exp\left[-{{\nu^{2}} \over {2R^{2}}}(g_{xki}^{2}+g_{yki}^ {2}+g_{zki}^{2})\right]\cr &\displaystyle\geq\sigma^{2}_{i}(2\pi)^{{{2} / {3}}}l_{ix}l_{iy}l_{iz}\exp \left(-{{\nu^{2}}/ {2}}\right),&(49)}]$

where the inequality holds because the maximum value of (g_xki, g_yki, g_zki) is R². Hence, we can see that ν controls the minimum spectral density, or alternatively we could view it as controlling the proportion of the volume under the spectral density function we wish the basis functions to cover. Taking this view, ν = 1 gives ∼68%, ν = 2 gives ∼95% and ν = 3 gives ∼99.7% volume coverage. In this work, we use ν = 3.5, corresponding to approximately 0.9996% coverage of the volume under the spectral density function.

Continuing with this reasoning, we can view R as governing the resolution with which the basis functions cover the spectral density. Whilst larger R will result in a better approximation to the covariance function it also increases the computational cost and, in general, will have diminishing returns in terms of error reduction. A suggestion is to increase R, subject to computational limits, whilst observing a substantial reduction in residuals or improvement in the out-of-sample log likelihood – described in detail in the next section. For both the simulation and real data experiments in this work we have used R = 5, which results in a total of m = 38 used basis functions for each of the six covariance functions, k_i(x, x′), i = 1, 2,…, 6. Increasing R beyond this was found to give minimal improvement.

To complete the regression scheme, we now discuss the selection of the hyperparameters l_ix, l_iy, l_iz and σ_i, which at this stage are the only unknowns in the formulation.

4.5. Hyperparameter selection

The hyperparameters, l_ix, l_iy, l_iz and σ_i, for the posterior conditional distribution can be determined through optimization (Rasmussen, 2003). Typically, this is done by either maximizing the log marginal likelihood or using a cross-validation approach and maximizing the `out-of-sample' log likelihood, i.e. the likelihood of observing a set of measurements not used in the regression, $[\tilde{{\bf y}}]$ . Following the work by Gregg et al. (2020 ), which demonstrates that maximizing the out-of-sample log likelihood yields better results for line integral measurements, we determine the hyperparameters by solving

$[\eqalignno{\Theta_{*}& = \mathop{{\rm {arg \,max}}}\limits_{\Theta}\log p_ {\Theta}(\tilde{{\bf y}}|{\bf y}) = \mathop{\rm {arg \,max}}\limits_{\Theta} -0.5\log\det{\bb C}[\tilde{ {\bf y}},\tilde{{\bf y}}|{\bf y}] \cr & \quad-0.5\left(\tilde{{\bf y}}-{\bb E}[\tilde{{\bf y}}| {\bf y}]\right)^{\rm T}{\bb C}[\tilde{{\bf y}},\tilde{{\bf y}}| {\bf y}]^{-1}\left(\tilde{{\bf y}}-{\bb E}[\tilde{{\bf y}}| {\bf y}]\right).&(50)}]$

where Θ is a vector holding the hyperparameters introduced in (35) and $[\log p_ {\Theta}(\tilde{{\bf y}}|{\bf y})]$ is the out-of-sample log likelihood. By extension of (47), we have that

$[\eqalign{{\bb E}[\bar{\tilde{{\bf y}}}| {\bf y}] & = {\bb E}[\tilde{{\bf y}}]+\boldphi_{\tilde{y}}\big(\boldphi_{y}^{\rm T}{\bb C}[{\bf e},{\bf e}]^{-1} \boldphi_{y}\cr & \quad+{\bf S}^{-1}\big)^{-1} \boldphi_{y}^{\rm T} {\bb C}[{\bf e},{\bf e}]^{-1}\left({\bf y}-{\bb E} [{\bf y}]\right),\cr {\bb C}[\tilde{{\bf y}},\tilde{{\bf y}}| {\bf y}] & = \boldphi_{\tilde{y}}\left(\boldphi_{y}^{\rm T} {\bb C}[{\bf e},{\bf e}]^{-1}\boldphi_{y}+ {\bf S}^{-1}\right)^{-1}\boldphi_{\tilde{y}}^{\rm T}+{\bb C} [{\bf e},{\bf e}].}\eqno(51)]$

Note that it is not essential that a global optimum is found in this procedure; in fact, in many cases, setting the hyperparameters to some reasonable fixed values will produce excellent reconstructions. In the case of scanning 3DXRD we have found that setting the hyperparameters uniformly to the grain diameter gives reasonable results and can serve as a good initial guess for optimization.

5. Validation

To validate the presented regression method we have generated simulated scanning 3DXRD data using a previously developed algorithm (Henningsson, 2019 ). This tool has been used with success in the past (cf. Hektor et al., 2019; Henningsson et al., 2020) and can provide an understanding of the limitations and benefits of scanning 3DXRD reconstruction methods. Briefly, the simulation input is specified as a set of cubic single-crystal voxels featuring individual strains and orientations together with an experimental setup. We refer the reader to Henningsson (2019) for additional details on the simulation algorithm, with an undocumented implementation available via https://github.com/FABLE-3DXRD/S3DXRD/. Strain reconstructions from generated diffraction data were compared with ground-truth input strain as well as an additional reconstruction method described by Henningsson et al. (2020). This reconstruction method, previously referred to as `algebraic strain refinement' (ASR), uses a voxel basis for strain reconstruction and can, in short, be described as solving a global WLSQ problem. This least-squares approach operates from the same average directional strain data as the presented GP method.

5.1. Single-crystal simulation test case

Diffraction from a tin (Sn) grain subject to a nonuniform strain tensor field has been simulated for the nonconvex grain topology depicted in Fig. 4.

Figure 4
Grain topology input for diffraction simulation coloured by corresponding input Euler angle field in units of degrees. The top row represents central cuts through the 3D renderings below, as indicated by the red lines.

The grain was assigned an orientation field by introducing linear gradients in the three Euler (Bunge notation) angles, φ₁, Φ, φ₂, as

$[\varphi_{1} = \Phi = \varphi_{2} = {{\pi}\over{180}}\left(45+{{x}\over{130v}}+{ {z}\over{24v}}\right),\eqno(52)]$

where v = 5 µm was the used voxel size and the grain origin was set at the grain centroid in the xy plane and at the bottom edge of the grain in z (Fig. 4). The maximum grain size in each dimension x, y and z was 26, 26 and 13 voxels, respectively.

The strain field was defined by a set of Maxwell stress functions, which are a subset of the more general class of Beltrami stress functions,

$[\bar{\boldPhi} = \left[\matrix{A(x,y,z)&B(x,y,z)&C(x,y,z)&0&0&0}\right]^{\rm T}.\eqno(53)]$

To achieve a relatively simple, but not trivial, strain field the functions A, B and C were selected as a cubic polynomial:

$[\eqalignno{ A = B = C & = \rho_{1}(x-t_{x})^{3}+\rho_{2}(y-t_{y})^{3}+\rho_{3}(z-t_{z})^{3 }\cr & \quad +\rho_{4}(x-t_{x})(y-t_{y})(z-t_{z}).&(54)}]$

The stress was converted to strain by the elastic compliance matrix C, as

$[\left[\matrix{\epsilon_{xx}\cr \epsilon_{yy}\cr \epsilon_{zz}\cr \epsilon_{xy}\cr \epsilon_{xz}\cr \epsilon_{yz}}\right] = {\bf C}^{-1}{\cal B}\!\!\!\!{\cal B}\bar{ \boldPhi} = {\bf C}^{-1}\left[\matrix{6\rho_{3}(z-t_{z})+6\rho_ {2}(y-t_{y})\cr 6\rho_{1}(x-t_{x})+6\rho_{3}(z-t_{z})\cr 6\rho_{2}(y-t_{y})+6\rho_{1}(x-t_{x})\cr \rho_{4}(t_{z}-z)\cr \rho_{4}(t_{y}-y)\cr \rho_{4}(t_{x}-x)}\right].\eqno(55)]$

Numerical values of the constants ρ₁, ρ₂, ρ₃, ρ₄, t_x, t_y and t_z are presented in Table 1.

Table 1
Strain field parameters for diffraction simulation in units of µm

ρ₁	ρ₂	ρ₃	ρ₄	t_x	t_y	t_z
100	100	100	1000	10	10	0

The elasticity matrix for single-crystal tin was taken from Darbandi et al. (2013 ) (Table 2) and converted from Voigt notation to the used strain vector notation.

Table 2
Elasticity constants for single-crystal tin in units of GPa converted from Voigt notation as given by Darbandi et al. (2013)

C₁₁	C₂₂	C₃₃	C₄₄	C₅₅	C₆₆	C₁₂	C₁₃	C₂₃
72.3	72.3	88.4	48.0	44.0	44.0	59.4	35.8	35.8

Parameters presented in Table 3 were used to define the experimental setup of the simulation.

Table 3
Experimental parameters used in single-grain simulation, corresponding to the results presented in Fig. 5

Wavelength	0.22 Å
Sample-to-detector distance	163 mm
Detector pixel size	50 × 50 µm
Detector dimensions	2048 × 2048 pixels
Beam size	5 × 5 µm
ω rotation interval	[0, 180°]
Δω step length	1°
Maximum grain size in x	130 µm
Maximum grain size in y	130 µm
Maximum grain size in z	65 µm

The unit cell in Table 4 was used to define a strain-free lattice state.

Table 4
Relaxed reference lattice parameters

a	b	c	α	β	γ
5.81127 Å	5.81127 Å	3.17320 Å	90.0°	90.0°	90.0°

The generated diffraction patterns were analysed on a per-z-slice basis using ImageD11 (Wright, 2005) to compute scattering vectors and average crystal orientations for each z slice. The grain shape was then reconstructed on the basis of the normalized diffraction peak intensities using filtered backprojection (Poulsen & Schmidt, 2003). Next, the diffraction data were converted to average directional strains, as described in Section 3, and input into the WLSQ and GP reconstruction methods. The final reconstructed strain tensor fields are illustrated together with simulation ground-truth and residual fields in Fig. 5. The corresponding root-mean-square errors (RMSEs), mean absolute errors (MAEs) and maximum absolute errors for the residual fields are given in Table 5.

Table 5
Root-mean-square errors, mean absolute errors and maximum absolute errors for the residual fields presented in Fig. 5

The result of the Gaussian process regression is compared with the weighted least-squares fit (WLSQ). Values are unitless (strain) and on the same scale (10⁻⁴) as in Fig. 5.

	RMSE		MAE		Maximum absolute error
Strain	GP	WLSQ	GP	WLSQ	GP	WLSQ
ε_xx	1.322	2.076	1.101	1.737	2.791	7.42
ε_yy	1.042	1.371	0.846	0.999	3.856	6.094
ε_zz	0.887	1.489	0.769	1.157	1.914	6.736
ε_xy	1.122	1.511	0.955	1.172	2.778	6.306
ε_xz	0.24	1.04	0.198	0.798	1.506	4.85
ε_yz	0.48	0.958	0.399	0.742	1.34	4.652

Figure 5
3D rendering of strain reconstructions for WLSQ and GP regression approaches. The top row defines the simulation ground truth as described in equation (55)

, with each column featuring a different strain component. The surface of the voxelated grain is presented, together with a pulled-out interior spherical cut centred at the grain centroid with a diameter of 50 µm. The corresponding coordinate systems are depicted in the bottom left of the figure. Three separate colormaps have been assigned to enhance contrast for the various fields. However, units of strain remain the same across plots (×10⁻⁴). The residual field is defined as the difference between the ground truth and the reconstructed strain field.

Hyperparameters were optimized using the L-BFGS-B algorithm, as implemented in SciPy (Jones et al., 2001 ), with a maximum of ten line-search steps per iteration. Gradients were computed using automatic differentiation as implemented in PyTorch (Paszke et al., 2019 ). In the first optimization iteration all hyperparameters were uniformly set to the grain radius. The convergence of the optimization is displayed in Fig. 6. The smoothness constraints for the WLSQ in the xy plane were set to 2.5 × 10⁻⁴, limiting the maximum absolute difference in each strain tensor component between two neighbouring voxels [further details are provided by Henningsson et al. (2020)].

Figure 6
Negative cross-validation log likelihood reduction during hyperparameter optimization for the simulated grain presented in Figs. 4

and 5

. Optimization was conducted using the L-BFGS-B algorithm as implemented in SciPy with a maximum of ten line-search steps per iteration. Gradients were computed using automatic differentiation as implemented in PyTorch.

To assess how well the two methods (WSLQ and GP) utilize data, the MAE and RMSE of the reconstructed strain fields, as a function of the number of input measurement integrals, has been investigated. By measurements we here refer to the integral values, y_j, as defined in (24), together with their associated vectors ( $[{\bf p}_{0}]$ , $[{\hat{\boldkappa}}]$ , $[{\hat{\bf n}}]$ ). Measurements were permuted randomly and input into the WLSQ and GP reconstruction in initial sample sizes of 1, 2, 3, 4 and 5%, after which the sample size was increased in steps of 5% as indicated by the markers in Fig. 7. Since the GP hyperparameter optimization is a nonconvex problem, the quality of any found local minima may vary between runs, and a better local minimum is not guaranteed with a larger measurement set owing to the different topology of the cost function. Thus, in order not to obscure the convergence rate of the GP method, we have selected to present results using fixed optimized hyperparameters found using 10% of the measurements as well as for non-optimized hyperparameters, set uniformly to the grain diameter. The resulting MAE and RMSE for the reconstructed residual fields were computed and averages over the six strain components were formed. The performance as a function of input measurements can be assessed by visual inspection of Fig. 7.

Figure 7
Average root-mean-square error (squares) and mean absolute error (stars) for the simulated grain presented in Figs. 4

and 5

as a function of used percentage of measurements. The performance of the Gaussian process regression (red and blue filled lines) is compared with that of the weighted least squares (black dashed lines). The RMSE and MAE were computed from the residual fields and averaged over the six reconstructed strain components to produce a scalar measure per reconstruction. Each point in the plot corresponds to a full 3D strain reconstruction using a random subset of the measured data.

5.2. Embedded tin grain

To further compare the GP and WLSQ reconstruction methods, analysis of a previously studied columnar tin grain has been included. This additional analysis further serves to show that the presented method is computationally feasible for state-of-the-art scanning 3DXRD data sets. Including hyperparameter optimization, the GP reconstruction was performed on a single CPU (Intel i7-8700K CPU @ 3.70 GHz) in 18 min and 9 s. As mentioned in Section 4.4, the computational complexity scales as $[{\cal O}(Nm^{2})]$ , where m is the number of basis functions and N the number of measurements. The corresponding runtime using fixed precomputed hyperparameters was 3.5 s. The data for this example from Hektor et al. (2019) and the input experimental parameters are identical to those presented in Table 3 except for the beam size, which was 0.25 µm. Similarly, the relaxed lattice state was as defined in Table 4. In the original experiment, the X-ray beam was scanned across the xy plane, producing a space-filling map of measurements. However, owing to time constraints, the data were collected for every second z layer, as seen in the rightmost column of Fig. 8. The reader is referred to the original publication (Hektor et al., 2019) for further information on the experimental setup, sample and preliminary data analysis.

Figure 8
Reconstructed strain field using WLSQ (left column) and the GP method (middle column) of a columnar tin grain embedded within a polycrystalline sample. The rightmost column depicts the estimated uncertainty of the GP reconstruction. The 3D surface of the voxelated grain is presented together with a pull-out enlarged interior spherical cut with its centre at the grain centroid and a radius of 1 µm. Two separate colormaps have been assigned to enhance contrast for the various fields. However, the units of strain remain the same across plots (×10⁻⁴).

As the GP method uses a nonlocal basis representation of the strain field, as defined in equation (39), interpolation between measured slices is an automatic feature of the method. For the WLSQ method, although some interpolation scheme could be selected, we have chosen to present the raw reconstructions. This also highlights the added benefit of the selected basis for the GP method. Hyperparameter optimization and smoothness constraints for the WLSQ method were applied and selected as in Section 5.1.

6. Discussion

Comparison of the true and predicted fields in Fig. 5 for the two methods indicates that the reconstructions captured well the simulated input strain state. For all strain components in Table 5, both the RMSE and MAE are of the order of the expected experimentally limited strain resolution (10⁻⁴). We note, however, that the GP has consistently lower RMSE, MAE and maximum absolute errors in comparison with the WLSQ. The enhanced performance is attributed to the joint effect of the equilibrium prior, optimized correlation kernel and nonlocal basis selection.

The results of Table 5 indicate that, in general, the strain tensor z components enjoy more accurate reconstructions than the xy components. This observation is in line with previous work (Margulies et al., 2002; Lionheart & Withers, 2015; Henningsson et al., 2020) and is explained by the nonuniform sampling of strain taking place in scanning 3DXRD. The GP regression quantifies this phenomenon via the reconstructed standard deviation fields (Fig. 5, bottom row). Indeed the uncertainty in the predicted mean is elevated for the xx and yy components, which show similar patterns to the residual fields.

On the performance of the two methods, Fig. 7 indicates that fewer measurements are needed for the GP compared with the WLSQ approach whilst achieving a more accurate result. Little reduction in the RMSE and MAE is seen for the GP after about 50% of the measurements have been introduced (about 20% for the optimized GP version). This could imply that it is possible to retrieve approximations to the full strain tensor field from reduced scanning 3DXRD data sets. This could be attractive as scanning 3DXRD typically has time-consuming measurement sequences. From Fig. 7 it is also clear that the final errors in reconstruction will be nonzero. This is so because the error in reconstruction is made up of both bias and variance. While the variance can be reduced by adding more measurements, the bias is due to systematic errors arising from incorrect model assumptions such as the line integral approximation, the truncated covariance basis series expansion, the Taylor series expansion related to the strain measure, the directional approximation of $[{\hat{\boldkappa}}]$ and possibly further unknown sources. Since the bias cannot be removed by adding more measurements, the reconstruction error will face a lower nonzero bound.

It is evident that the reconstructed fields have maximum uncertainties at the boundary of the grain, as can be seen from the cutout spheres of Figs. 5 and 8. The elevated standard deviation at the grain surface is explained by the tomographic measurement procedure, which has an increasing measurement density towards the grain centroid. Furthermore, as measurements do not exist outside of the grain, points lying on the grain surface will, in some sense, have a reduced number of points that they are correlated with. In addition to these effects, the selected line beam approximation may have an impact on the grain boundary errors. If the full 3D profile of the beam had been used instead, a higher number of scans that partially graze the grain boundary could have been included in the analysis, thus increasing the measurement density at the boundary. In the current model, if a scan has a geometric centre that does not intersect the grain, it has no impact on the reconstruction, even though the full 3D beam may have some overlap with the grain. The main challenge with using a full 3D beam profile, rather than the line approximation, is to maintain analytical expressions during integration of the partial derivatives of the basis functions over the illuminated domain.

The predicted strain field of the columnar tin grain of Fig. 8 shows similar patterns for the two regression methods. The uncertainty is again seen to be reduced on the interior of the grain, and the posterior standard deviation is of the order of the experimental strain resolution of 10⁻⁴. This validates the applicability of GP regression on real state-of-the-art scanning 3DXRD synchrotron data.

6.1. Outlook

Two future potential improvements to strain predictions should be mentioned. Firstly, the selection of covariance function, although restricted to give a positive definite covariance matrix, is not unique; other selections may outperform the squared-exponential kernel used here. Secondly, for polycrystalline samples, additional prior knowledge of grain boundary strain could be extracted by considering the total sample grain map and that tractions must cancel on the interfaces [i.e. incorporating and extending the work of Hendriks, Gregg et al. (2019)]. Two challenges with this exist: (i) the uncertainty in reconstructed grain shapes leading to uncertainty in the interface normal and (ii) uncertainty in the per-point grain orientation leading to uncertainty in the grain compliance. The first of these challenges may be addressed by using near-field techniques (Viganò et al., 2016 ) in conjunction with scanning 3DXRD to achieve higher-resolution grain maps.

7. Conclusions

Intragranular strain estimation from scanning 3DXRD data using a Gaussian process is shown to provide a new and effective strain reconstruction method. By selecting a continuous differentiable Fourier basis for the Beltrami stress functions, a static equilibrium prior can be incorporated into the reconstruction, guaranteeing that the predicted strain field will satisfy the balance of both angular and linear momentum. The regression procedure results in a per-point estimated mean strain and per-point standard deviations, providing new means of estimating the per-point uncertainty of the reconstruction. Furthermore, the proposed method incorporates the spatial structure of the strain field by making use of a generic covariance function, optimized by maximizing the out-of-sample log likelihood. With the introduction of these three features, the equilibrium prior, the per-point uncertainty quantification and the optimized spatial smoothness constraints, the proposed regression method addresses weaknesses discussed in previously proposed reconstruction methods. Specifically, in comparison with a previously proposed weighted least-squares approach, it is found, from numerical simulations, that the Gaussian process regression consistently produces reconstructions with lower root-mean-square errors, mean absolute errors and maximum absolute errors across strain components. Moreover, it is shown that the reconstruction error as a function of the number of available measurements is reduced for the Gaussian process.

APPENDIX A

Error related to measurement approximations

To demonstrate the accuray of the approximations made in (15), we investigate the error associated with the fact that both strain and crystal orientation may vary along the ray path. To do this we must consider that the strain computed from (15) is further assigned to planes with approximate normals given by (17). Thus, there is a twofold error source to capture in the following analysis, arising partly from the integrated Taylor series expansion,

$[y = {{1}\over{V}}\int\limits_{{\cal L}}\hat{\boldkappa}^{\rm T}\boldepsilon_{\omega}\hat{\boldkappa}\,{\rm d}v\simeq 1-{{ \langle{\bf G}_{\omega}\rangle^{\rm T}{\bf G}^{(0)}_{\omega}}\over{({\bf G}^{(0)}_{\omega})^{\rm T}{\bf G}^{(0)}_{\omega}}},\eqno(56)]$

and partly from assigning the average strain value, y, to an incorrect plane normal,

$[{\hat{\boldkappa}}\simeq{{\langle{\bf G}_{\omega}\rangle}/{ ||\langle{\bf G}_{\omega}\rangle||}},\eqno(57)]$

which in reality is not fixed but warps across the crystal [ $[{\hat{\boldkappa}} = {\hat{\boldkappa}}({\bf x})]$ ].

To compute the error in strain we consider first the true average strain for a single line-integral measurement, y_true, existing in a fixed direction, $[{\hat{\boldkappa}}]$ :

$[y_{\rm true} = {{1}\over{L}}\int\limits_{{\cal L}}{\hat{\boldkappa}}^{\rm T} \boldepsilon_{\omega}{\hat{\boldkappa}}\,{\rm d}s = {{1}\over{L}}{\hat{ \boldkappa}}^{\rm T}\left(\int\limits_{{\cal L}}\boldepsilon_{\omega}\,{\rm d}s\right){ \hat{\boldkappa}} = {\hat{\boldkappa}}^{\rm T}\langle\boldepsilon_{ \omega}\rangle{\hat{\boldkappa}},\eqno(58)]$

where the average strain tensor 〈ε_ω〉 is unknown from the experiment. We stress that we are interested in the true strain for a fixed $[{\hat{\boldkappa}}]$ since this is the normal that the approximation of (56) will eventually be assigned to. The sought absolute error now becomes

$[e_{y} = y-y_{\rm true} = 1-{{\langle{\bf G}_{\omega}\rangle^{\rm T} {\bf G}^{(0)}_{\omega}}\over{({\bf G}^{(0)}_{\omega})^{\rm T}{\bf G}^{(0)}_{ \omega}}}-{\hat{\boldkappa}}^{\rm T}\langle\boldepsilon_{\omega} \rangle{\hat{\boldkappa}}.\eqno(59)]$

For a given set of Miller planes, strain field, Euler angle field, reference unit cell and integration domain $[{\cal L}]$ , equation (59) can be evaluated. To do so, two integrations must be performed, yielding individually 〈G_ω〉 and 〈ε_ω〉. In the following we attempt to characterize (59) for a fixed reference unit cell (Table 4) while letting the remaining parameters vary according to a stochastic model described below. The goal is to study the distribution of the absolute errors as a function of increasing levels of intragranular strain and mosaicity to understand the limitations of the proposed measurement approximations.

We consider a spherical grain of fixed radius R₀ = 1.0 centred at the origin and define a random integration domain, $[{\cal L}]$ , as

$[\eqalign{ & {\bf x} = {\bf p}_{0}+s{\hat {\bf n}},\quad s\in[0,L],\cr & {\bf p}_{0} = R_{0}\left[\matrix{\cos(a_{1})\sin(a_{2})\cr \sin(a_{1})\sin(a_{2})\cr \cos(a_{2})}\right],\quad a_{1}\sim{\cal U}(0,2\pi),\quad a_{2}\sim {\cal U}(0,\pi),\cr & {\hat{\bf n}} = \left[\matrix{\cos(b_{1})\sin(b_{2})\cr \sin(b_{1})\sin(b_{2})\cr \cos(b_{2})}\right],\quad b_{1}\sim{\cal U}(0,2\pi),\quad b_{2}\sim {\cal U}(0,\pi), }\eqno(60)]$

where L is determined by the sphere line intersection, p₀ is a random uniform point on the sphere surface and $[{\hat{\bf n}}]$ is a unit vector also drawn from a random uniform distribution denoted $[{\cal U}(\cdot,\cdot)]$ . We further define G_hkl as

$[{\bf G}_{hkl} = [h,k,l]^{\rm T},\quad h,k,l\sim{\cal U}(-7,7),\eqno(61)]$

where the distribution has been limited to the interval [−7, 7] to represent typical scanning 3DXRD data sets. Pseudo-random strain and orientation fields are introduced by constructing random samples of superimposed Fourier waves:

$[\eqalign{&\epsilon_{xx}({\bf x}) = {{A_{\rm s}}\over{f^{\rm max} _{0}}}f_{0}({\bf x}),\quad\epsilon_{yy}({\bf x}) = {{A_{\rm s}}\over{f^ {\rm max}_{1}}}f_{1}({\bf x}),\quad\epsilon_{zz}({\bf x}) = {{A_{\rm s }}\over{f^{\rm max}_{2}}}f_{2}({\bf x}),\cr & \epsilon_{xy}({\bf x}) = {{A_{\rm s}}\over{f^{\rm max}_{3}}}f_{3}({\bf x}),\quad\epsilon_{xz}({\bf x}) = {{A_{\rm s}}\over {f^{\rm max}_{4}}}f _{4}({\bf x}),\quad\epsilon_{yz}({\bf x}) = {{A_{\rm s}}\over{f^{\rm max}_ {5}}}f_{5}({\bf x}),\cr & \varphi_{1}({\bf x}) = {{A_{\rm e}}\over{f^{\rm max}_{6}}}f_{6}({\bf x}),\quad\Phi({\bf x}) = {{A_{\rm e}}\over{f^{\rm max}_{7}}}f_{7}({\bf x}),\quad\varphi_{2}({\bf x}) = {{A_{\rm e}}\over{f^{\rm max}_{8}}}f_{ 8}({\bf x}),\cr & f_{k}({\bf x}) = \!\textstyle\sum\limits^{n = 25}_{i = 1}\!c_{i}\sin\left[\,f_{ix}(x +p_{ix})\right]\sin\left[\,f_{iy}(y+p_{iy})\right]\sin\left [\,f_{i}z(z+p_{iz}) \right],\cr & c_{i} \sim{\cal U}(-1,1),\quad f_{ix},f_{iy},f_{iz}\sim {\cal U}\left({{10}\over{R_{0}}},{{1}\over{2R_{0}}}\right),\cr & p_{ix},p_{iy},p_{iz}\sim{\cal U}(-R_{0},R_{0}),\,\, f^{\rm max} _{k} = \mathop{\rm {arg \, max}}\limits_{{\bf x}}|f_{k}({\bf x})|,\,\, {\bf x}^{\rm T}{\bf x}\,\lt\,R_{0}^{2},}\eqno(62)]$

where the two scale parameters A_s and A_e regulate the maximum difference between any two points within the field and allow for the strain field (scaled by A_s) and orientation field (scaled by A_e) to vary on different scales simultaneously. The necessary normalizing factors $[f^{\rm max}_{k}]$ were computed by sampling the fields on equidistant grids of ∼1000 points and selecting the maximum absolute value. Typical fields generated by the model can be seen in Fig. 9.

Figure 9
Typical Euler angle (top row) and strain (bottom and middle rows) fields generated by the stochastic model presented in equation (62)

. The presented fields exist on a spherical domain for which a central cut slice has been presented above (z = 0). The maximum difference parameters of the field were A_e = 1.4 and A_s = 75 × 10⁻⁴.

Using the above stochastic model we have performed repetitive sampling of 1000 line integral measurements for each of 100 grain states while letting A_s and A_e successively increase between samples. The results of this analysis are presented in the histograms of Fig. 10, where each histogram corresponds to a total of 100 000 data points (100 × 1000) and a fixed maximum field variation (A_s, A_e). The average strain 〈ε_ω〉 and diffraction vector 〈G_ω〉 involved in (59) were computed by first-order numerical integration using a total of 20 integration points along each domain $[{\cal L}]$ . The reference orientation matrix, U₀, was computed by averaging over ∼1000 equally spaced points of the sphere.

Figure 10
Absolute error (59)

computed for the stochastic model defined through equations (60)

, (61)

and (62)

. For each histogram 1000 random line integral measurements have been sampled from each of 100 spherical crystal states, resulting in a total of 100 000 data points per histogram. The maximum field difference in Euler angles and strain (A_e and A_s) increases from bottom to top and left to right, respectively, as indicated by the figure labels. (Note that the maximum counts of the bottom-left plot have been clipped in order to facilitate equal axes between subplots while maintaining good visibility of the histograms.)

The results of Fig. 10 show that the error in (59) increases with the heterogeneity of both orientation and strain state. For small strains (≤50 × 10⁻⁴) and moderate mosaic spreads (≤1°) the largest errors are a few times that of the experimental resolution limit (10⁻⁴) and the bulk (>95%) of measurements are contained within [±10⁻⁴] units of strain. For samples featuring larger mosaicity (>1.0°) and strain variation (>75 × 10⁻⁴) the approximation starts to break down. At such elevated levels of deformation, however, one has to first consider if the small-strain approximation made in (7) is still valid.

Specifically, for the synthetic data set presented in this paper (Figs. 4 and 5) we conclude that the input strain and orientation fields will give rise to a negligible error, e_y < 10⁻⁴. Likewise for the tin grain presented in Fig. 8, on the basis of the observed diffraction peak spread in ω, the mosaic spread is <1.0° and thus e_y is negligible.

Acknowledgements

The authors are grateful for the beamtime provided by the ESRF, beamline ID11, where the diffraction data were collected (Hektor et al., 2019). The authors also thank Stephen Hall for valuable input on the manuscript during the final stages of writing the paper.

References

Alpers, A., Poulsen, H. F., Knudsen, E. & Herman, G. T. (2006). J. Appl. Cryst. 39, 582–588. Web of Science CrossRef CAS IUCr Journals Google Scholar
Borbely, A., Renversade, L., Kenesei, P. & Wright, J. (2014). J. Appl. Cryst. 47, 1042–1053. Web of Science CrossRef CAS IUCr Journals Google Scholar
Darbandi, P., Bieler, T., Pourboghrat, F. & Lee, T. (2013). J. Electron. Mater. 42, 201–214. Web of Science CrossRef CAS Google Scholar
Edmiston, J. K., Barton, N. R., Bernier, J. V., Johnson, G. C. & Steigmann, D. J. (2011). J. Appl. Cryst. 44, 299–312. Web of Science CrossRef CAS IUCr Journals Google Scholar
Gregg, A., Hendriks, J., Wensrich, C. & O'Dell, N. (2020). Nucl. Instrum. Methods Phys. Res. B, 480, 67–77. Web of Science CrossRef CAS Google Scholar
Hayashi, Y., Hirose, Y. & Seno, Y. (2015). J. Appl. Cryst. 48, 1094–1101. Web of Science CrossRef CAS IUCr Journals Google Scholar
Hayashi, Y., Setoyama, D., Hirose, Y., Yoshida, T. & Kimura, H. (2019). Science, 366, 1492–1496. Web of Science CrossRef CAS PubMed Google Scholar
Hayashi, Y., Setoyama, D. & Seno, Y. (2017). Mater. Sci. Forum, 905, 157–164. CrossRef Google Scholar
Hektor, J., Hall, S. A., Henningsson, N. A., Engqvist, J., Ristinmaa, M., Lenrick, F. & Wright, J. P. (2019). Materials, 12, 446. Web of Science CrossRef Google Scholar
Hendriks, J., Gregg, A., Wensrich, C. & Wills, A. (2019). Strain, 55, e12325. Web of Science CrossRef Google Scholar
Hendriks, J. N., Wensrich, C. M. & Wills, A. (2020). Strain, 56, e12341. Web of Science CrossRef Google Scholar
Hendriks, J., Wensrich, C., Wills, A., Luzin, V. & Gregg, A. (2019). Nucl. Instrum. Methods Phys. Res. B, 444, 80–90. Web of Science CrossRef CAS Google Scholar
Henningsson, A. (2019). Student paper, Lund University, Sweden, https://lup.lub.lu.se/student-papers/record/8972668. Google Scholar
Henningsson, N. A., Hall, S. A., Wright, J. P. & Hektor, J. (2020). J. Appl. Cryst. 53, 314–325. Web of Science CrossRef CAS IUCr Journals Google Scholar
Jidling, C., Hendriks, J., Wahlström, N., Gregg, A., Schön, T. B., Wensrich, C. & Wills, A. (2018). Nucl. Instrum. Methods Phys. Res. B, 436, 141–155. Web of Science CrossRef CAS Google Scholar
Jones, E., Oliphant, T., Peterson, P., et al. (2001). SciPy, https://www.scipy.org/. Google Scholar
Lauridsen, E. M., Schmidt, S., Suter, R. M. & Poulsen, H. F. (2001). J. Appl. Cryst. 34, 744–750. Web of Science CrossRef CAS IUCr Journals Google Scholar
Lionheart, W. R. B. & Withers, P. J. (2015). Inverse Probl. 31, 045005. Web of Science CrossRef Google Scholar
Margulies, L., Lorentzen, T., Poulsen, H. & Leffers, T. (2002). Acta Mater. 50, 1771–1779. Web of Science CrossRef CAS Google Scholar
Markussen, T., Fu, X., Margulies, L., Lauridsen, E. M., Nielsen, S. F., Schmidt, S. & Poulsen, H. F. (2004). J. Appl. Cryst. 37, 96–102. Web of Science CrossRef CAS IUCr Journals Google Scholar
Moscicki, M., Kenesei, P., Wright, J., Pinto, H., Lippmann, T., Borbély, A. & Pyzalla, A. (2009). Mater. Sci. Eng. A, 524, 64–68. Web of Science CrossRef Google Scholar
Nervo, L., King, A., Wright, J. P., Ludwig, W., Reischig, P., Quinta da Fonseca, J. & Preuss, M. (2014). J. Appl. Cryst. 47, 1402–1416. Web of Science CrossRef CAS IUCr Journals Google Scholar
Oddershede, J., Schmidt, S., Poulsen, H. F., Sørensen, H. O., Wright, J. & Reimers, W. (2010). J. Appl. Cryst. 43, 539–549. Web of Science CrossRef CAS IUCr Journals Google Scholar
Ottosen, N. S. & Ristinmaa, M. (2005). The Mechanics of Constitutive Modeling. Kidlington: Elsevier. Google Scholar
Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J. & Chintala, S. (2019). Advances in Neural Information Processing Systems 32, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. dAlché-Buc, E. Fox & R. Garnett, pp. 8024–8035. Curran Associates. https://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf. Google Scholar
Poulsen, H. (2004). PhD thesis, Risø National Laboratory, Roskilde, Denmark. Google Scholar
Poulsen, H. F. & Fu, X. (2003). J. Appl. Cryst. 36, 1062–1068. Web of Science CrossRef CAS IUCr Journals Google Scholar
Poulsen, H. F., Nielsen, S. F., Lauridsen, E. M., Schmidt, S., Suter, R. M., Lienert, U., Margulies, L., Lorentzen, T. & Juul Jensen, D. (2001). J. Appl. Cryst. 34, 751–756. Web of Science CrossRef CAS IUCr Journals Google Scholar
Poulsen, H. F. & Schmidt, S. (2003). J. Appl. Cryst. 36, 319–325. Web of Science CrossRef CAS IUCr Journals Google Scholar
Rasmussen, C. E. (2003). Advanced Lectures on Machine Learning, Lecture Notes in Computer Science, Vol. 3176, pp. 63–71. Berlin, Heidelberg: Springer. Google Scholar
Reischig, P. & Ludwig, W. (2020). Curr. Opin. Solid State Mater. Sci. 24, 100851. Web of Science CrossRef Google Scholar
Schmidt, S. (2014). J. Appl. Cryst. 47, 276–284. Web of Science CrossRef CAS IUCr Journals Google Scholar
Sharma, H., Huizenga, R. M. & Offerman, S. E. (2012). J. Appl. Cryst. 45, 693–704. Web of Science CrossRef CAS IUCr Journals Google Scholar
Solin, A. & Särkkä, S. (2020). Stat. Comput. 30, 419–446. Web of Science CrossRef Google Scholar
Viganò, N., Tanguy, A., Hallais, S., Dimanov, A., Bornert, M., Batenburg, K. J. & Ludwig, W. (2016). Sci. Rep. 6, 20618. Web of Science PubMed Google Scholar
Wright, J. (2005). ImageD11, https://github.com/FABLE-3DXRD/ImageD11/. Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.

JOURNAL OF
APPLIED
CRYSTALLOGRAPHY

ISSN: 1600-5767

Volume 54| Part 4| August 2021| Pages 1057-1070

https://doi.org/10.1107/S1600576721005112

Open

access

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Search IUCr Journals		doi		Advanced search
Author		volume	page

research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Intragranular strain estimation in far-field scanning X-ray diffraction using a Gaussian process

1. Introduction

2. Diffraction measurements

2.1. Experimental acquisition

2.2. Laue equations and scattering notation

2.3. Grain mapping

3. Measurement model

3.1. Strain revealing transformations

3.2. Tensorial ray transform

3.3. Estimated uncertainty

4. Regression procedure

4.1. Gaussian process regression

4.2. Equilibrium prior

4.3. Equilibrium posterior distribution

4.4. Finite basis approximations

4.5. Hyperparameter selection

5. Validation

5.1. Single-crystal simulation test case

5.2. Embedded tin grain

6. Discussion

6.1. Outlook

7. Conclusions

APPENDIX A

Error related to measurement approximations

Acknowledgements

References

research papers