An efficient system matrix factorization method for scanning diffraction based strain tensor tomography

Henningsson, A.; Hall, S.A.

doi:10.1107/S2053273323008136

research papers

FOUNDATIONS
ADVANCES

ISSN: 2053-2733

Volume 79| Part 6| November 2023| Pages 542-549

https://doi.org/10.1107/S2053273323008136

Open

access

An efficient system matrix factorization method for scanning diffraction based strain tensor tomography

Axel Henningsson ^a ^* and Stephen A. Hall ^a

^aDivision of Solid Mechanics, Lund University, Ole Römersväg 1, Lund, Sweden
^*Correspondence e-mail: axel.henningsson@solid.lth.se

Edited by D. A. Keen, STFC Rutherford Appleton Laboratory, United Kingdom (Received 2 February 2023; accepted 18 September 2023; online 29 September 2023)

Diffraction-based tomographic strain tensor reconstruction problems in which a strain tensor field is determined from measurements made in different crystallographic directions are considered in the context of sparse matrix algebra. Previous work has shown that the estimation of the crystal elastic strain field can be cast as a linear regression problem featuring a computationally involved assembly of a system matrix forward operator. This operator models the perturbation in diffraction signal as a function of spatial strain tensor state. The structure of this system matrix is analysed and a block-partitioned factorization is derived that reveals the forward operator as a sum of weighted scalar projection operators. Moreover, the factorization method is generalized for another diffraction model in which strain and orientation are coupled and can be reconstructed jointly. The proposed block-partitioned factorization method provides a bridge to classical absorption tomography and allows exploitation of standard tomographic ray-tracing libraries for implementation of the forward operator and its adjoint. Consequently, RAM-efficient, GPU-accelerated, on-the-fly strain/orientation tensor reconstruction is made possible, paving the way for higher spatial resolution studies of intragranular deformation.

Keywords: X-ray diffraction; strain tensor; tomography; diffraction imaging.

1. Introduction

Diffraction-based strain tomography is an experimental technique deployed for estimation of the six-component elastic strain tensor field, $[\boldvarepsilon({\bf x})]$ , within the bulk of polycrystalline aggregates. Whether used with X-rays (Hektor et al., 2019 ; Korsunsky et al., 2005 ; Lionheart & Withers, 2015 ) or neutrons (Hendriks et al., 2020 ), the method offers a unique possibility to probe the internal heterogeneity of the strain in dense materials in a non-destructive way. In essence, the measured diffraction signal from the specimen can be reduced to average strains along line integral domains across a sample volume. Each of these scalar strain measures, $[\gamma_{j}]$ , can be associated to a spatial sampling direction, $[{\boldkappa}_{j}]$ , which in general varies between measurements. Considering a set of $[j =1,\ldots,m]$ such measurements, it is possible to construct a global linear system of equations,

$[{\bf A}{\bf s} ={\boldgamma},\eqno(1)]$

where $[{\bf s}]$ holds the basis coefficients of a decomposed strain tensor field and $[{\boldgamma}]$ is a vector with all measurements [the formation of $[{\boldgamma}]$ from raw diffraction images and the decoupling of the crystal strain from orientation are discussed by Henningsson & Hendriks (2021 )]. The rows of the system matrix, $[{\bf A}]$ , are required to contain the integral weights of the strain tensor basis functions combined with non-linear combinations of the components of $[{\boldkappa}_{j}]$ . Using measurements only from a single axis of rotation, there exist no known, closed-form, direct back-projection algorithms to recover the strain tensor field $[{\bf s}]$ . This motivates the need for iterative solvers, which may seem to require the assembly of $[{\bf A}]$ . Unfortunately, the storage of $[{\bf A}]$ can be very RAM inefficient and the assembly routines needed to construct $[{\bf A}]$ involve ray-tracing through the strain tensor volume. Indeed, direct storage of the forward operator is unfeasible for high-resolution scalar, absorption-based, tomography. Considering a fixed resolution, strain tomography of symmetric strain tensor fields offers no relaxation in this respect as the number of non-zeros in the system matrix, $[{\bf A}]$ , is a factor six greater compared with scalar tomography. As a result, several existing reconstruction methods have been cast in settings with few strain tensor basis functions limiting the achievable reconstruction resolution (Henningsson et al., 2020 ; Henningsson & Hendriks, 2021; Hendriks et al., 2020). Similarly, in scanning 3D X-ray diffraction (scanning-3DXRD) microscopy applications, it is common to collapse the rich 2D pixel intensity distribution of the recorded diffraction peaks to a single centre of gravity prior to the pursuit of strain reconstruction (Hayashi et al., 2015 ).

We present a system matrix factorization for strain tensor tomography in which the forward operator, $[{\bf A}]$ , can be implemented as a weighted sum of scalar forward projections as

$[{\bf A}{\bf s} = \textstyle\sum\limits^{i = 6}_{i = 1}{\bf S}_{i} {\cal P}{\bf a}_{i},\eqno(2)]$

where $[{\cal P}]$ is a scalar forward projection operator, $[{\bf a}_{1},\ldots,{\bf a}_{6}]$ are the six individual components of the strain tensor field and $[{\bf S}_{1},\ldots,{\bf S}_{6}]$ are diagonal weight matrices. This factorization allows for RAM-efficient, on-the-fly implementations to be easily achieved with existing tomographic libraries (for scalar projection). Additionally, our proposed factorization allows for access to GPU-accelerated implementations commonly deployed in scalar tomography to facilitate large sparse iterative solvers (Palenstijn et al., 2011 ; van Aarle et al., 2015 , 2016 ).

For illustrative purposes we have selected to present our derivations in the context of strain reconstruction and for the experimental setup of scanning-3DXRD. The methodology is, however, also applicable for other neutron and X-ray scanning diffraction experiments given that a fixed axis of rotation is used and that the diffraction peak centre-of-mass positions can be accurately measured (typically in far-field geometry). The key ingredient in our derivation is the linearity of the diffraction model, which allows us to rearrange the order of the involved operators. In contrast to the far-field diffraction setting considered in this paper, near-field diffraction methods (Reischig & Ludwig, 2020 ) model the full detector intensity distribution of the diffraction peaks rather than the peak centroid positions. As a result, the forward operator in near-field diffraction models depends non-linearly on the intragranular strain and orientation. To highlight that our factorization method is applicable to multiple models, as long as they fall within the class of far-field diffraction, we derive and demonstrate (in Appendix C) a factorization similar to that of equation (2) for a previously suggested diffraction model that features coupling between the intragranular strain and orientation. This factorization enables efficient reconstruction of the full intragranular deformation field, including both strain and orientations.

2. Per-ray factorization

Given an unknown, symmetric, second-order strain tensor field,

$[{\boldvarepsilon}({\bf x}) = \left[\matrix{\varepsilon_{1}({\bf x})&\varepsilon_{4}({\bf x})&\varepsilon_{5}({\bf x})\cr \varepsilon_{4}({\bf x})&\varepsilon_{2}({\bf x})&\varepsilon_{6 }({\bf x})\cr \varepsilon_{5}({\bf x})&\varepsilon_{6}({\bf x})&\varepsilon_{3}({\bf x})}\right],\eqno(3)]$

defined on a 3D spatial domain, $[{\bf x}= [x\,y\,z]^{\rm T}]$ , we shall consider measurements of the average strain, $[\gamma_{j}]$ , on the line integral domain $[{\cal R}_{j}]$ as

$[\gamma_{j} ={{1}\over{L_{j}}}\int\limits_{{\cal R}_{j}}{\boldkappa}_{j}^{\rm T}\boldvarepsilon({\bf x}){\boldkappa}_{j}\,{\rm d}{\bf x},\eqno(4)]$

where $[{\boldkappa}_{j} = [\kappa_{1}\,\kappa_{2}\,\kappa_{3}]^{\rm T}]$ is a unit normal vector that describes the sampled strain direction and L_j is the ray intersection path length measured over the compact support of $[{\boldvarepsilon}]$ . For scanning-3DXRD the formation of $[\gamma_{j}]$ from the raw diffraction image data has been described elsewhere (Henningsson & Hendriks, 2021). Following a flattened vector format similar to that of Henningsson & Hendriks (2021) we find the alternative measurement model

$[\gamma_{j} = \textstyle\int\limits_{{\cal R}_{j}}\bar{{\boldkappa}}_{j}^{\rm T}\bar{ {\boldvarepsilon}}({\bf x})\,{\rm d}{\bf x},\eqno(5)]$

where

$[\bar{{\boldvarepsilon}}({\bf x}) = \left[\matrix{\varepsilon_{1}({\bf x})\cr \varepsilon_{3}({\bf x})\cr \varepsilon_{2}({\bf x})\cr\varepsilon_{4}({\bf x})\cr \varepsilon_{5}({\bf x})\cr \varepsilon_{6}({\bf x})\cr }\right],\quad\bar{{\boldkappa}} = {{1}\over{L_{j}}}\left[\matrix{\kappa^{2}_{1j}\cr \kappa^{2}_{2j}\cr \kappa^{2}_{3j}\cr2\kappa_{1j}\kappa_{2j}\cr 2\kappa_{1j}\kappa_{3j}\cr 2\kappa_{2j}\kappa_{3j}\cr}\right].\eqno(6)]$

Let $[{\bar{\boldvarepsilon}}]$ be decomposed on $[{\bf x}]$ with n basis functions $[\varphi_{l}({\bf x})]$ as

$[{\bar{\boldvarepsilon}}({\bf x}) = \textstyle\sum\limits^{n}_{l = 1}\varphi_{l}({\bf x}){\boldalpha}_{l},\eqno(7)]$

where the basis coefficients $[{\boldalpha}_{l}]$ are defined as

$[{\boldalpha}_{l} = \left[\matrix{\alpha_{1l}&\alpha_{2l}&\alpha_{3l}&\alpha_{4l}&\alpha_{5l}&\alpha_{6l}}\right]^{\rm T}.\eqno(8)]$

In the following we select $[\varphi_{l}]$ to represent an equidistant grid of pixels such that $[\varphi_{l}({\bf x}) = 1]$ when $[{\bf x}]$ is in pixel number l and $[\varphi_{l} = 0]$ otherwise. By insertion of (7) into (5) we have

$[\gamma_{j} =\textstyle\int\limits_{{\cal R}_{j}}\bar{{\boldkappa}}_{j}^{\rm T}\sum\limits^{n}_{l = 1}\varphi_{l}({\bf x}){\boldalpha}_{l}\,{\rm d}{\bf x}.\eqno(9)]$

Reordering the integral and sum we can write

$[\gamma_{j} = \bar{{\boldkappa}}_{j}^{\rm T}\textstyle\sum\limits^{n}_{l = 1}{\boldalpha}_ {l}\int\limits_{{\cal R}_{j}}\varphi_{l}({\bf x})\,{\rm d}{\bf x}.\eqno(10)]$

We now introduce the vector $[{\bf w}_{j}]$ which contains the scalar weights of the ray integral with respect to basis functions,

$[{\bf w}_{j} = \left[\matrix{\int_{{\cal R}_{j}}\varphi_{1}(\bf{x})\,{\rm d}{\bf x}\cr\int_{{\cal R}_{j}}\varphi_{2}(\bf{x})\,{\rm d}{\bf x}\cr \vdots\cr \int_{{\cal R}_{j}}\varphi_{n}(\bf{x})\,{\rm d}{\bf x}}\right].\eqno(11)]$

Using the weights, $[{\bf w}_{j}]$ , we may form a matrix projection operator that projects the six components of the strain field along a single ray path as

$[{\bf R}_{j} = \underbrace{\left[\matrix{{\bf w}_{j}^{\rm T}& {\bf 0}&{\bf 0}&\ldots&{\bf 0}\cr {\bf 0}&{\bf w}_{j}^{\rm T}&{\bf 0}&\ldots&{\bf 0}\cr \vdots&\ddots&\ddots&\ddots&\vdots\cr {\bf 0}&\ldots&{\bf 0}&{\bf 0}&{\bf w}_{j}^{\rm T}}\right]}_{6\times 6n}.\eqno(12)]$

Additionally we introduce the vectors

$[{\bf a}_{1} = \left[\matrix{\alpha_{11}\cr \alpha_{12}\cr \vdots\cr\alpha_{1n}}\right],{\bf a}_{2} = \left[\matrix{\alpha_{21}\cr \alpha_{22}\cr\vdots\cr \alpha_{2n}\cr }\right]\ldots{\bf a}_{6} =\left[\matrix{\alpha_{61}\cr \alpha_{62}\cr \vdots\cr \alpha_{6n}\cr}\right],\eqno(13)]$

and stack the basis coefficients of the unknown strain tensor field in a single column vector as

$[{\bf s} =\underbrace{\left[\matrix{{\bf a}_{1}\cr {\bf a}_{2}\cr {\bf a}_{3}\cr {\bf a}_{4}\cr {\bf a}_{5}\cr {\bf a}_{6}\cr }\right]}_{6n\times 1}.\eqno(14)]$

We can now facilitate a fully vectorized and discretized format of the measurement model, equation (4), as

$[\gamma_{j} =\bar{{\boldkappa}}_{j}^{\rm T}{\bf R}_{j}{\bf s}.\eqno(15)]$

To arrive at a global format, in which several measurements, $[\gamma_{j}]$ , are considered simultaneously, we introduce the vector

$[\boldgamma =\left[\matrix{\gamma_{1}&\gamma_{2}&\ldots&\gamma_{m}}\right]^{\rm T}.\eqno(16)]$

Stacking the matrices $[\bar{{\boldkappa}}_{j}^{\rm T}]$ and $[{\bf R}_{j}]$ in the same fashion,

$[{\bf K} =\underbrace{\left[\matrix{\bar{{\boldkappa}}_{1}^{\rm T}& {\bf 0}&{\bf 0}&\ldots&{\bf 0}\cr {\bf 0}&\bar{{\boldkappa}}_{2}^{\rm T}&{\bf 0}&\ldots& {\bf 0}\cr \vdots&\ddots&\ddots&\ddots&\vdots\cr {\bf 0}&\ldots&{\bf 0}&{\bf 0}&\bar{{\boldkappa}}_ {m}^{\rm T}}\right]}_{m\times 6m},\quad{\bf V} =\underbrace{\left[\matrix{{\bf R}_{1}\cr {\bf R}_{2}\cr \vdots\cr {\bf R}_{m}}\right]}_{6m\times 6n},\eqno(17)]$

we find the global matrix formulation as

$[\boldgamma = {\bf K}{\bf V}{\bf s}.\eqno(18)]$

We note that in equation (18) the matrix $[{\bf A} = {\bf K}{\bf V}]$ is factorized in two terms: $[{\bf K}]$ , which contains information on the directional sampling of the strain field, and $[{\bf V}]$ , which holds information on the projections of the sampled fields.

3. Hexa-block-diagonal form

In scalar tomography the forward projection operator, $[{\cal P}]$ , is commonly block-partitioned over a series of projection views, $[{\cal P}_{i}]$ , as

$[{\cal P} =\left[\matrix{{\cal P}_{1}\cr {\cal P}_{2}\cr \vdots\cr {\cal P}_{k}}\right],\eqno(19)]$

where each projection view, $[{\cal P}_{i}]$ , represents an ordered set of parallel line integrals defined over a single scalar field. In contrast, we note that the ray integrals contained in $[{\bf V}]$ , as defined by equations (12) and (17), are neither ordered in complete views nor defined over a single scalar field. We therefore seek to reorder and partition the rays in $[{\bf V}]$ in a way that will allow our projection operator to be easily implemented using standard tomographic libraries. To this end, we note that the set of rows in $[{\bf V}]$ separated by a fixed multiple 6, with start at row number 1, forms the block-partitioned matrix

$[\underbrace{\left[\matrix{{\bf L}&{\bf 0}&{\bf 0}& {\bf 0}&{\bf 0}&{\bf 0}}\right]}_{m\times 6n},\eqno(20)]$

where $[{\bf L}]$ is now acting on the single scalar field $[{\bf a}_{1}]$ . If the measurements in $[\boldgamma]$ are selected to be stacked in complete projection views, we find that $[{\bf L} = {\cal P}]$ . Since the initial selected ordering of measurements in $[\boldgamma]$ is arbitrary, we shall assume that this ordering has been selected. Now, by simply repeating the row shifting operation with increasing row starting index, $[1,2,\ldots,6]$ , it is possible to mutate $[{\bf V}]$ into the block-diagonal matrix form,

$[{\bf P} = \underbrace{\left[\matrix{{\cal P}&{\bf 0}&{\bf 0}&{\bf 0}&{\bf 0}&{\bf 0}\cr {\bf 0}&{\cal P}&{\bf 0}&{\bf 0}& {\bf 0}&{\bf 0}\cr {\bf 0}&{\bf 0}&{\cal P}&{\bf 0}& {\bf 0}&{\bf 0}\cr {\bf 0}&{\bf 0}&{\bf 0}&{\cal P}& {\bf 0}&{\bf 0}\cr {\bf 0}&{\bf 0}&{\bf 0}&{\bf 0}&{\cal P}&{\bf 0}\cr {\bf 0}&{\bf 0}&{\bf 0}&{\bf 0}&{\bf 0}& {\cal P}\cr }\right]}_{6m\times 6n},\eqno(21)]$

which contains the reordered rows of $[{\bf V}]$ . Naturally, to maintain the global formulation in equation (17), we are required to now also modify $[{\bf K}]$ . The shifting of the rows of $[{\bf V}]$ , therefore, requires a corresponding shifting of columns in $[{\bf K}]$ , leading to the block-partitioned matrix

$[{\bf S} = \underbrace{\left[\matrix{{\bf S}_{1}& {\bf S}_{2}&{\bf S}_{3}&{\bf S}_{4}&{\bf S}_{5}& {\bf S}_{6}}\right]}_{m\times 6m},\eqno(22)]$

with diagonal blocks

$[{\bf S}_{1} = {\rm Diag}\left({{\kappa^{2}_{1j}}\over{L _{j}}}\right),\quad {\bf S}_{2} = {\rm Diag}\left({{ \kappa^{2}_{2j}}\over{L_{j}}}\right),\eqno(23)]$

$[{\bf S}_{3} = {\rm Diag}\left({{\kappa^{2}_{3j}}\over{L _{j}}}\right),\quad {\bf S}_{4} = {\rm Diag}\left({{ 2\kappa_{1j}\kappa_{2j}}\over{L_{j}}}\right),\eqno(24)]$

$[{\bf S}_{5} = {\rm Diag}\left({{2\kappa_{1j}\kappa _{3j}}\over{L_{j}}}\right),\quad {\bf S}_{6} = {\rm Diag}\left({{2 \kappa_{2j}\kappa_{3j}}\over{L_{j}}}\right).\eqno(25)]$

It is now possible to write

$[\boldgamma = {\bf S P s} = \textstyle\sum\limits^{i = 6}_{i = 1 }{\bf S}_{i}{\cal P}{\bf a}_{i}.\eqno(26)]$

In this factorization the execution of the forward operator, $[{\bf A} = {\bf S P}]$ , corresponds to six scalar forward projections followed by the application of $[{\bf S}]$ , which, due to its diagonal form, presents a modest 6m multiplications and additions. The implementation of $[{\cal P}]$ can be directly achieved by any ray-tracing library, e.g. the ASTRA-toolbox (Palenstijn et al., 2011; van Aarle et al., 2015, 2016). The implementation of $[{\bf S}]$ is trivial and, owing to its diagonal form, there is no need to assemble the matrix, as it suffices to store the six vectors of diagonal weights. Since the six projections being executed in $[{\bf P}]$ are independent, we note that the resulting arrays, $[{\bf a}_{1},\ldots,{\bf a}_{6}]$ , may be stacked and projected in parallel on a GPU. To service the diffraction imaging community, and to illustrate how equation (26) can be put to use to achieve an easily implementable GPU-accelerated diffraction model, we provide an open-source demo Python code at https://github.com/AxelHenningsson/flyxdm.

4. Generalizations

For the sake of clarity, we derived equation (26) in the setting of strain reconstruction. This setting features scalar measurements, $[\gamma_{j}]$ , which simplifies the exposition and allows us to focus on the core rearrangement of equations necessary to arrive at our block-partitioned factorized format. The same algebraic manipulations can be used to factorize a wider class of linear far-field diffraction models. We demonstrate the generality of our matrix factorization method in Appendix C where we have pursued an extended diffraction model originally suggested by Henningsson et al. (2020). In this alternative setting the intragranular orientation field is jointly reconstructed with the strain tensor field and the measurement associated to the ray integral is vector valued rather than scalar.

To reconstruct a target field, $[{\bf s}]$ , in practice, it is often desirable to introduce a measurement weight matrix, $[{\bf W}]$ , that describes the measurement precision. For instance, in the work of Henningsson et al. (2020) a diagonal weight matrix was used to reconstruct strain in a weighted least-squares sense. We note that our factorization is indifferent to the introduction of $[{\bf W}]$ and the global equation system would in practical application be extended as

$[{\bf W}\boldgamma = {\bf W S P s} ,\eqno(27)]$

where $[{\bf W}^{2} = {\boldSigma}_{\boldgamma}^{-1}]$ is the covariance of the measurements, $[\boldgamma]$ .

Another practical concern is the incorporation of constraints on the solution vector, $[{\bf s}]$ . One popular approach is to modify the basis of the unknown target field to encode the prior knowledge. To exploit our factorized format in these settings, we suggest introducing a rendering matrix, $[{\bf N}]$ , that maps the basis coefficients, $[{\bf q}]$ , from the constrained basis set back to the pixel basis coefficients, $[{\bf s}]$ , as

$[{\bf s} = {\bf Nq}.\eqno(28)]$

The resulting global equations now become

$[{\bf W}\boldgamma = {\bf W S P N q},\eqno(29)]$

where the columns of $[{\bf N}]$ can be interpreted as pixel images of the selected basis set. As the forward operator, $[{\bf W S P N}]$ , and the adjoint operator, $[{\bf N}^{\rm T}{\bf P}^{\rm T}{\bf S}^{\rm T}{\bf W}^{\rm T}]$ , still feature the desired multiplicative block-partitioned split between $[{\bf P}]$ and $[{\bf S}]$ , we conclude that the results in equation (26) can be exploited in a wide range of applications.

As a final note, on the topic of generalizations, we would like to mention that, just as the tensor components of the target field can be stacked into a 3D volume and projected in parallel on a GPU card, one may instead consider stacking grain slices into a volume and projecting each tensor component separately. This modification, reconstructing a full grain volume rather than a grain slice, has no impact on the algebraic format of equation (26). The rendering matrix, $[{\bf N}]$ , is then computing the coefficients of a set of voxels that are projected as a 3D volume by $[{\bf P}]$ .

5. Demonstration

To demonstrate the memory benefits that can be achieved using equation (26) compared with assembling and storing the sparse matrix, $[{\bf A}]$ , we consider a single-crystal diffraction simulation case study. The simulation is described in detail in Appendix A, together with illustrations of the reconstructions achieved when exploiting the format of equation (26) during regression (Appendices B and C). The supplementary code used to generate the simulation data as well as the reconstructions is openly available at https://github.com/AxelHenningsson/flyxdm.

In Fig. 1 we present the number of megabytes of computer RAM necessary to compute $[\boldgamma = {\bf A s}]$ using either a fully assembled sparse matrix $[{\bf A}]$ or, alternatively, the factorization $[{\bf A} = {\bf S P}]$ , where $[{\bf P}]$ is represented using pre-existing, on-the-fly, projection operators, available in the ASTRA-toolbox (van Aarle et al., 2015). Considering that the results presented in Fig. 1 represent the reconstruction of a single grain slice using 500 projection views (each corresponding to a diffraction event), it is evident that parallel, high-resolution, full volume/sample reconstructions are unfeasible using an assembled format of $[{\bf A}]$ . For instance, reconstructing, in parallel, a single, cubic-shaped grain volume, with a cross-sectional resolution of 256×256 pixels from ∼300 unique (with respect to Miller index) diffraction peaks would require ∼1 TB of computer RAM storage.

Figure 1
Number of megabytes of computer RAM necessary to compute $[\boldgamma = {\bf A s}]$ using either an assembled sparse matrix (dashed line) or, alternatively, the discussed factorization in conjunction with an on-the-fly projection operator, $[{\bf P}]$ (solid line). Note that the benchmarks correspond to a single-crystal grain slice and 500 diffraction peaks (projection views) as described in Appendix A

6. Conclusion

We have presented a system matrix factorization for strain tensor tomography in which the directional sampling of the strain tensor field is separated from the tomographic projection operator. The proposed format allows for the exploitation of standard tomographic ray-tracing libraries in the implementation of the forward operator. We have also shown how our factorization method can be generalized for other diffraction models, for example one in which strain and orientation are coupled. We have provided an openly available GPU implementation of the approach and demonstrated the computational efficiency of our factorization method through application to a model example. By enabling RAM-efficient, GPU-accelerated, on-the-fly strain/orientation tensor reconstruction, our results facilitate higher spatial resolution studies of intragranular deformation.

APPENDIX A

Demonstration example

To demonstrate the discussed matrix factorization in a practical application we have included a single-crystal X-ray diffraction simulation case study. Diffraction data were forward modelled from a 2D grain slice of α-quartz (SiO₂) subject to a spatially varying strain tensor field, $[{\boldvarepsilon}({\bf x})]$ , as well as a misorientation field, $[{\bf U}({\bf x})]$ , where $[{\bf U}({\bf x})]$ is the local crystal orientation matrix. The synthetic strain tensor field can be viewed in Fig. 2 together with the Bunge Euler angle variation, which was defined as

$[\eqalignno{\Delta\varphi_{1} &= \varphi_{1}-\langle\varphi_{1} \rangle,&\cr \Delta\phi & = \phi-\langle\phi\rangle,&\cr \Delta\varphi_{2}& = \varphi_{2}-\langle\varphi_{2}\rangle,&(30)}]$

where $[\langle\cdot\rangle]$ denotes volume average.

Figure 2
Strain (bottom) and mosaicity (top) in a simulated 2D grain slice of α-quartz (SiO₂). The Bunge Euler angles are displayed as variations around their respective mean values.

Using the space group of quartz (P3221) together with an X-ray wavelength of λ = 0.2846 Å, a total of 500 reflections were randomly (uniformly) selected in the Bragg angle interval θ = [4°, 13°] to be included in the simulation. Using the Laue equations,

$[{\bf G} = {\bf U B G}_{hkl},\eqno(31)]$

where the matrix $[{\bf B}]$ maps from the integer reciprocal-space Miller indices, $[{\bf G}_{hkl} = [h\,k\,l]^{\rm T}\in{\bb Z}^{3 \times 1}]$ , to crystal coordinates, diffraction vectors $[{\bf G} = [G_{1}\, G_{2}\, G_{3}]^{\rm T}]$ were formed. The theoretical diffraction vectors were then corrupted with zero mean Gaussian noise as

$[\eqalignno{& G_{1j}\leftarrow G_{1j}+n_{1j},\quad n_{1j}\sim {\cal N}(0,\sigma^{2}_{1j})&\cr & G_{2j}\leftarrow G_{2j}+n_{2j},\quad n_{2j}\sim{\cal N}(0, \sigma^{2}_{2j})&\cr & G_{3j}\leftarrow G_{3j}+n_{3j},\quad n_{3j}\sim{\cal N}(0, \sigma^{2}_{3j}),&(32)}]$

where $[\sigma_{1j} = \sigma_{2j} = \sigma_{3j} = 10^{-3}]$ for 90% of the measurements while $[\sigma_{1j} = \sigma_{2j} = \sigma_{3j} = 10^{-2}]$ for the remaining 10%, emulating the presence of outliers. The 10% selected for outlier noising was selected randomly uniformly from the full diffraction data set.

Using the ASTRA-toolbox we have implemented the matrix factorization derived in this paper in operator format. This means that we never need to form the explicit sparse matrices involved in the forward model, but, instead, implement an on-the-fly operator that operates on an input vector to produce the system matrix vector dot product (and the corresponding adjoint product). This operator implementation is only possible thanks to the algebraic results that constitute the contribution of this paper. The code used to implement our matrix factorization (on an NVIDIA GPU architecture), produce the simulated data and reconstruct the strain (and later orientation) field is openly available as a demo Python library https://github.com/AxelHenningsson/flyxdm. In the same supplementary demo code, we provide detailed comments on the simulation and reconstruction setup and provide the code used to generate all figures found in these Appendices (including the generalizations to orientation fields made in Appendix C).

APPENDIX B

Strain reconstruction

In this Appendix, we use the simulation data in conjunction with the Taylor expansion described by Henningsson & Hendriks (2021) and convert the diffraction vectors, $[{\bf G}]$ , into measurements of directional strain with the aim of reconstructing the intragranular strain field. In Appendix C we describe how the diffraction vectors can be used without transform to reconstruct strain and intragranular orientation variations jointly, again exploiting our matrix factorization method. Note that while the aim in Appendix B is to demonstrate our matrix factorization method for the case of strain reconstruction, the simulated data still originate from a grain that features both intragranular strain and orientation variation. As discussed elsewhere (Henningsson et al., 2020; Henningsson & Hendriks, 2021), this is not a problem as long as the mosaicity of the grain is moderate.

In Fig. 3, the result of strain tensor reconstruction can be viewed. Here a radial basis expansion was used to construct $[{\bf N}]$ [in equation (29)]. The true noise covariance of the diffraction vectors was propagated through the Taylor expansion to construct the directional strain covariance yielding $[{\bf W}]$ . The small residuals and root-mean-squared errors (RMSEs) found in the bottom row of Fig. 3 are expected as a consequence of noise and model mismatch. For further details we refer the reader to the supplementary code.

Figure 3
Strain tensor reconstruction (REC) of a single-crystal α-quartz grain slice. The matrix factorization in equation (29)

has been used to reconstruct the strain field without the need to assemble the global system matrix. The top row ground-truth simulation input (GT) is to be compared with the middle row reconstructed strain field (REC). The bottom row shows the residual between reconstructed and true strain fields (RES).

APPENDIX C

Generalization to coupled orientation–strain models

We shall now consider generalizing our matrix factorization method to a diffraction model discussed by Henningsson et al. (2020, Section 6 equation 10-16). Considering the Laue equations (31) in integral form, we find

$[\langle{\bf G}_{j}\rangle = {{1}\over{L_{j}}}\int\limits_{{\cal R}_{j}} {\bf U B G}^{(j)}_{hkl}\,{\rm d}{\bf x},\eqno(33)]$

where $[\langle\cdot\rangle]$ denotes volume average. The target intragranular field of deformation is now generalized to include the local rigid-body rotations as

$[{\bf U B} = {\bf U}({\bf x}){\bf B}({\bf x}) = {\bf f}({\bf x})\in{\bb R}^{3\times 3}.\eqno(34)]$

In contrast to the small strain tensor, $[{\boldvarepsilon}]$ , $[{\bf f}]$ is not symmetric and the number of unknowns per point, $[{\bf x}]$ , in the grain is now increased to 9 (compared with six strain tensor components). To reach a similar factorization as that of (26), we start by introducing a flattened format of (33) as

$[\langle{\bf G}_{j}\rangle = {\bf H}_{j}\textstyle\int\limits_{{\cal R}_{j}} {\bar{\bf f}}({\bf x})\,{\rm d}{\bf x},\eqno(35)]$

where

$[{\bar{\bf f}} = \left[\matrix{f_{11}({\bf x})\cr f_{21}({\bf x})\cr f_{31}({\bf x})\cr f_{12}({\bf x})\cr f_{22}({\bf x})\cr f_{32}({\bf x})\cr f_{13}({\bf x})\cr f_{23}({\bf x})\cr f_{33}({\bf x})}\right]\eqno(36)]$

and

$[{\bf H}_{j} = {{1}\over{L_{j}}}\left[\matrix{{\bf I}h_{j}& {\bf I}k_{j}&{\bf I}l_{j}}\right],\eqno(37)]$

and $[{\bf I}]$ is the 3×3 identity matrix. Let us now decompose the target field, $[{\bar{\bf f}}({\bf x})]$ , on a finite basis as

$[{\bar{\bf f}}({\bf x}) = \textstyle\sum\limits^{n}_{l = 1}{\boldbeta}_{l} \varphi_{l}({\bf x}),\eqno(38)]$

where $[\varphi_{l}\in{\bb R}]$ are the scalar (pixel/voxel) basis functions and

$[{\boldbeta}_{l} = \left[\matrix{\beta_{1l}\cr \beta_{2l}\cr \beta_{3l}\cr \beta_{4l}\cr \beta_{5l}\cr \beta_{6l}\cr \beta_{7l}\cr \beta_{8l}\cr \beta_{9l}\cr }\right].\eqno(39)]$

Insertion of equation (38) into equation (35) yields after rearrangement that

$[\langle{\bf G}_{j}\rangle = {\bf H}_{j}\textstyle\sum\limits^{n}_{l = 1}{\boldbeta}_{l}\int\limits_{{\cal R}_{j}}\varphi_{l}({\bf x})\,{\rm d}{\bf x},\eqno(40)]$

where the linearity of the involved operators was used. We now note that equation (40) is a higher-dimensional copy of equation (10). Defining $[{\bf w}_{j}]$ according to equation (11), we find in analogy with equation (12) that

$[{\bf Q}_{j} = \underbrace{\left[\matrix{{\bf w}_{j}^{\rm T}& {\bf 0}&{\bf 0}&\ldots&{\bf 0}\cr {\bf 0}&{\bf w}_{j}^{\rm T}&{\bf 0}&\ldots&{\bf 0}\cr \vdots&\ddots&\ddots&\ddots&\vdots\cr {\bf 0}&\ldots&{\bf 0}&{\bf 0}&{\bf w}_{j}^{\rm T}}\right]}_{9\times 9n}.\eqno(41)]$

We now introduce a set of coefficient vectors,

$[{\bf b}_{1} = \left[\matrix{\beta_{11}\cr \beta_{12}\cr \vdots\cr \beta_{1n}}\right],{\bf b}_{2} = \left[\matrix{\beta_{21}\cr \beta_{22}\cr \vdots\cr \beta_{2n}}\right],\ldots,{\bf b}_{9} = \left[\matrix{\beta_{91}\cr \beta_{92}\cr \vdots\cr \beta_{9n}}\right],\eqno(42)]$

and define a partitioned global coefficient vector as

$[{\boldtau} = \left[\matrix{{\bf b}_{1}\cr {\bf b}_{2}\cr \vdots\cr {\bf b}_{9}}\right].\eqno(43)]$

The ray integral equation (40) can now be cast as

$[\langle{\bf G}_{j}\rangle = {\bf H}_{j}{\bf Q}_{j} {\boldtau}.\eqno(44)]$

We now note that each row in equation (44) has an identical algebraic format compared with equation (15). Splitting the diffraction vector measurements into three separate vectors as

$[{\bf d}_{1} = \left[\matrix{\langle{\bf G}_{11}\rangle\cr \langle{\bf G}_{12}\rangle\cr \vdots\cr \langle{\bf G}_{1m}\rangle}\right],\quad{\bf d}_{2} = \left[\matrix{\langle{\bf G}_{21 }\rangle\cr \langle{\bf G}_{22}\rangle\cr \vdots\cr \langle{\bf G}_{2m}\rangle}\right],\quad{\bf d}_{3} = \left[\matrix{\langle{\bf G}_{31 }\rangle\cr \langle{\bf G}_{32}\rangle\cr \vdots\cr \langle{\bf G}_{3m}\rangle}\right],\eqno(45)]$

we are free to use the same arguments of row and column permutation as described in Section 3 to arrive at

$[{\bf d}_{1} = {\bf M}_{1}{\bf P}^{\prime}{\boldtau},\eqno(46)]$

where $[{\bf P}^{\prime}]$ now is a 9m×9n block-diagonal projection matrix in direct analogy to equation (21) and the block-diagonal matrix $[{\bf M}_{1}]$ holds the Miller indices weighted by path length as

$[{\bf M}_{1} = \underbrace{\left[\matrix{{\bf h}&{\bf 0}& {\bf 0}&{\bf k}&{\bf 0}&{\bf 0}&{\bf l}& {\bf 0}&{\bf 0}}\right]}_{m\times 9m},\eqno(47)]$

with

$[\eqalignno{\bf{h} &= {\rm Diag}\left({{h_{j}}\over{L_ {j}}}\right),&\cr {\bf k}& = {\rm Diag}\left({{k_{j}}\over{L_{j}}}\right), &\cr {\bf l}& = {\rm Diag}\left({{l_{j}}\over{L_{j}}}\right). &(48)}]$

The global system of equations can now be written as

$[\left[\matrix{{\bf d}_{1}\cr {\bf d}_{2}\cr {\bf d}_{3} }\right] = \left[\matrix{{\bf M}_{1}\cr {\bf M}_{2}\cr {\bf M}_{3} }\right]{\bf P}^{\prime}{\boldtau},\eqno(49)]$

with

$[\eqalignno{{\bf M}_{2} &= \left[\matrix{{\bf 0}& {\bf h}&{\bf 0}&{\bf 0}&{\bf k}&{\bf 0}& {\bf 0}&{\bf l}&{\bf 0}}\right],&\cr {\bf M}_{3} &= \left[\matrix{{\bf 0}&{\bf 0}& {\bf h}&{\bf 0}&{\bf 0}&{\bf k}&{\bf 0}& {\bf 0}&{\bf l}}\right].&(50)}]$

Alternatively, denoting the measurement vector as

$[{\bf d} = \underbrace{\left[\matrix{{\bf d}_{1}\cr {\bf d}_{2}\cr {\bf d}_{3}}\right]}_{3m\times 1}\eqno(51)]$

and the Miller sampling matrix as

$[{\bf M} = \underbrace{\left[\matrix{{\bf M}_{1}\cr {\bf M}_{2}\cr {\bf M}_{3}}\right]}_{3m\times 9m},\eqno(52)]$

we arrive at our final factorized diffraction model,

$[{\bf d} = {\bf M}{\bf P}^{\prime}{\boldtau}.\eqno(53)]$

The forward pass in equation (53) is defined by nine separate (scalar) projection operations followed by nine multiplications with the diagonal blocks ( $[{\bf h},{\bf k},{\bf l}]$ ) of $[{\bf M}]$ . This factorization therefore admits the same computational benefits that are discussed for the decoupled strain model in the main paper. The discussion on generalizations held in Section 4, introducing a weight matrix $[{\bf W}]$ and a change of basis matrix $[{\bf N}]$ , is likewise applicable to equation (53). To verify that our derivations are correct we dedicate the following section to applying equation (53) to our demonstration example presented in Appendix A.

C1. Application to demonstration example

For completeness we have implemented and solved equation (53) in our demo supplementary code for the same demonstration example that was considered in Appendices A and B. The same radial basis expansion as described in Appendix B was used during regression. Likewise, the true noise covariance matrix was used to solve for the unknown radial basis coefficients in a weighted least-squares sense. The resulting maximum likelihood reconstruction (REC) can be viewed in the middle row of Fig. 4 together with the ground truth (GT) and residual (RES).

Figure 4
Coupled strain–orientation reconstruction (middle row) in a single slice of α-quartz. The simulated ground-truth (GT) field and corresponding data are described in Appendix A

. The residual field (RES) can be viewed in the bottom row. Note that the Bunge Euler angles (left) are displayed as a deviation from their respective mean values, allowing for a shared colorbar.

Funding information

This work was funded by Vetenskapsrådet – Röntgen Ångström Cluster, project No. 2017-06719.

References

Aarle, W. van, Palenstijn, W. J., Cant, J., Janssens, E., Bleichrodt, F., Dabravolski, A., De Beenhouwer, J., Joost Batenburg, K. & Sijbers, J. (2016). Opt. Express, 24, 25129–25147. Web of Science PubMed Google Scholar
Aarle, W. van, Palenstijn, W. J., De Beenhouwer, J., Altantzis, T., Bals, S., Batenburg, K. J. & Sijbers, J. (2015). Ultramicroscopy, 157, 35–47. Web of Science PubMed Google Scholar
Hayashi, Y., Hirose, Y. & Seno, Y. (2015). J. Appl. Cryst. 48, 1094–1101. Web of Science CrossRef CAS IUCr Journals 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 PubMed Google Scholar
Hendriks, J. N., Wensrich, C. M. & Wills, A. (2020). Strain, 56, e12341. Google Scholar
Henningsson, A. & Hendriks, J. (2021). J. Appl. Cryst. 54, 1057–1070. Web of Science CrossRef CAS IUCr Journals 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
Korsunsky, A. M., Vorster, W. J. J., Zhang, S. Y., Dini, D., Latham, D., Golshan, M., Liu, J., Kyriakoglou, Y. & Walsh, M. J. (2006). Acta Mater. 54, 2101–2108. Web of Science CrossRef CAS Google Scholar
Lionheart, W. R. B. & Withers, P. J. (2015). Inverse Probl. 31, 045005. Web of Science CrossRef Google Scholar
Palenstijn, W. J., Batenburg, K. J. & Sijbers, J. (2011). J. Struct. Biol. 176, 250–253. Web of Science CrossRef CAS PubMed Google Scholar
Reischig, P. & Ludwig, W. (2020). Curr. Opin. Solid State Mater. Sci. 24, 100851. Web of Science CrossRef Google Scholar