Diffraction tomography of strain

This paper examines whether it is possible to recover the three dimensional strain tensor field tomographically from neutron and x-ray diffraction data for polycrystalline materials.

The tomography of symmetric rank two tensor fields defined on domains in ${{\mathbb{R}}^{3}}$ is well studied, in particular comprehensive results are presented in [18]. So far these results have not been extensively used in practice. Of course any attempt to recover the six independent spatially varying components of a tensor field will inevitably demand the collection of a large quantity of data, and the necessary computational resources are required for its inversion. As the computational challenges fall to the advance of Moore's law experimentalists are revisiting the requirements for data collection. In the case of polarized light tomography for the study of strain in transparent photoelastic objects systems already exist [3, 21, 24] and indeed there are explicit analytical [10] as well as numerical [20], reconstruction algorithms for the Truncated Transverse Ray transform, which describes that system for sufficiently small strains.

Recent interest has been shown in accurate non-destructive high resolution measurement of strain in polycrystalline materials. Traditionally strain fields are mapped in one, two or three dimensions by defining a gauge volume by slits or collimators that is then translated step by step in one two or three dimensions to recover a map of a strain component (figures 1, a, b). One disadvantage of this method is that it is time consuming while for synchrotron diffraction, at least, the gauge is often much longer along the transmitted raypath than lateral to it (figures 1, a, b) thereby limiting the spatial resolution (see [23]). Korsunsky et al point out [8, 9] that if strain tomography were possible then this would enable the reconstruction of the strain field across slices or volumes in three dimensions potentially alleviating both of these issues.

Zoom In Zoom Out Reset image size

Here our aim is to take a fundamental look at the potential for strain tomography from a mathematical viewpoint, to examine what information can be recovered and to identify how best to extract information about strain from the information accumulated by a ray. In particular we will show that the average of the transversal component of the strain tensor along a ray is related to a moment of superimposed, distorted Debye–Scherrer rings that can be measured in the experiment. We propose a reconstruction of this component from moment data available for all rays orthogonal to a chosen direction. For six suitably chosen directions the strain tensor can be uniquely recovered and we give a simple criterion to test if the directions are suitable. We also give a reconstruction algorithm assuming complete data. We see that the method of [8] is essentially correct except they use a different average over the Debye–Scherrer ring.

By contrast, we show that at least in its simplest interpretation, Bragg edge neutron diffraction tomography for polychromatic beams [1, 2, 4, 15–17, 22] can only measure the change in path length of the beam through the material (i.e. the average through thickness strain) and is unaffected by the distribution of strain and so strain cannot be simply tomographically reconstructed using Bragg edges.

Let $C_{0}^{\infty }({{\mathbb{R}}^{3}};{{S}^{2}}{{\mathbb{R}}^{3}})$ be the space of smooth compactly supported functions on ${{\mathbb{R}}^{3}}$ with values in the six dimensional vector space ${{S}^{2}}{{\mathbb{R}}^{3}}$ of symmetric two-tensors on ${{\mathbb{R}}^{3}}$. The simplest ray transform we can take is to apply the scalar ray transform to each component $X:C_{0}^{\infty }({{\mathbb{R}}^{3}};{{S}^{2}}{{\mathbb{R}}^{3}})\to {{C}^{\infty }}({{\mathbb{R}}^{3}}\times \mathbb{R}_{0}^{3};{{S}^{2}}{{\mathbb{R}}^{3}})$

$$\begin{eqnarray}&&Xf(x,\xi )=\int _{-\infty }^{\infty }f(x+s\xi ){\rm d}s,\end{eqnarray} \tag{ 1 }$$

where $\mathbb{R}_{0}^{3}={{\mathbb{R}}^{3}}\setminus \{0\}$. Without loss of generality, when describing the ray $x+t\xi$, we can take ξ to be a unit vector and $x\in {{\xi }^{\bot }}$. Following Sharafutdinov [18] we also define the transverse ray transform $J:C_{0}^{\infty }({{\mathbb{R}}^{3}};{{S}^{2}}{{\mathbb{R}}^{3}})\to {{C}^{\infty }}({{\mathbb{R}}^{3}}\times \mathbb{R}_{0}^{3};{{S}^{2}})$

$$\begin{eqnarray}&&Jf(x,\xi )=\int _{-\infty }^{\infty }{{P}_{\xi }}f(x+s\xi ){\rm d}s,\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&{{P}_{\xi }}(f)={{\Pi }_{\xi }}f{{\Pi }_{\xi }},{\rm with}\;{{\Pi }_{\xi }}={\bf I}-\frac{1}{|\xi {{|}^{2}}}\xi {{\xi }^{T}},\end{eqnarray} \tag{ 3 }$$

is the projection on to the plane ${{\xi }^{\bot }}$, we identify rank two tensors with matrices for the purposes of multiplication and ${\bf I}$ is the identity matrix.

We also define the longitudinal ray transform $I:C_{0}^{\infty }({{\mathbb{R}}^{3}};{{S}^{2}}{{\mathbb{R}}^{3}})\to {{C}^{\infty }}({{\mathbb{R}}^{3}})$

$$\begin{eqnarray}&&If(x,\xi )=\int _{-\infty }^{\infty }{{\xi }^{T}}f(x+s\xi )\xi {\rm d}s.\end{eqnarray} \tag{ 4 }$$

While we have formulated these transforms on compactly supported tensor fields it is equally possible to extend them to the Schwartz space of test functions, and to distributions for example using the formal adjoint as is done for X in [14]. It is also simple to show that, like X, J and I are bounded maps between suitable Sobolev spaces, see for example [20].

Our integral transforms of tensor fields will often enable us to recover a diagonal component in some frame at each point. With this in mind we notice that a symmetric matrix $A\in {{S}^{2}}{{\mathbb{R}}^{3}}$ determines a quadratic form $Q(x)={{x}^{T}}Ax$ and it is well known that the off-diagonal elements of A can be determined from Q using the polarization identity

$$\begin{eqnarray}&&{{a}_{ij}}=(Q({{e}_{i}}+{{e}_{j}})-Q({{e}_{i}})-Q({{e}_{j}}))/2.\end{eqnarray} \tag{ 5 }$$

Here Q evaluated on the unit basis vectors ${{e}_{1}},{{e}_{2}},{{e}_{3}}$ as well as the vectors ${{e}_{1}}+{{e}_{2}},{{e}_{2}}+{{e}_{3}}$ and ${{e}_{3}}+{{e}_{1}}$, but more generally it is convenient to have geometric criteria that six non-zero vector x_i, $i=1,\ldots ,6$ are sufficient for $Q({{x}_{i}})$ to uniquely determine A. Of course one can simply say that the expressions

$$\begin{eqnarray}&&x_{i}^{T}A{{x}_{i}},\quad i=1,\ldots ,6\end{eqnarray} \tag{ 6 }$$

are independent as linear forms in the a_ij, or equivalently that the outer products ${{x}_{i}}x_{i}^{T}$ span ${{S}^{2}}{{\mathbb{R}}^{3}}$, but it is useful to have a geometric rather than algebraic criterion. Notice that if the equations are not independent there must be some non-zero symmetric matrix A in the null space of the linear operator defined by (6), with $x_{i}^{T}A{{x}_{i}}=0$ for all $i=1,\ldots ,6$. This is precisely the same as saying that the six points on the projective plane P² lie on a projective conic defined by (non-zero scalar multiples of) A. Indeed any five distinct points in P² lie on a unique projective conic so what is needed is a criterion to tell if the sixth point lies on the conic defined by the other five. We have exactly what we need in the classical Pascal's theorem (see figure 2) and its converse, from projective geometry. Recall that the projective plane can be considered as simply the unit sphere $\{x\in {{\mathbb{R}}^{3}}\left| \;|x|=1 \right.\}$ under the identification $x\sim -x$, hence points on P² are simply directions in ${{\mathbb{R}}^{3}}$. Projective lines in P² are the pairs of antipodal points lying on the intersection of a plane through the origin with the unit sphere. We now have Pascal's theorem and its converse [5].

Theorem 3.1. Six distinct points ABCDEF in the projective plane lie on a conic if and only if points of intersection of the lines AB with DE, CD with AF, and BC with FE are colinear.

One important special case is where the six directions x_i lie on some antipodal pair of small circles defined by $|x\cdot k|=a$ for some fixed unit vector k and $0\lt a\lt 1$. Such points do lie on a projective conic and so A is not determined by the $Q({{x}_{i}})$. On the other hand given five distinct points on such a pair or circles and a sixth that does not, A is determined by the $Q({{x}_{i}})$

A crystal consists of a periodic spatial arrangement of atoms (which we take to be located at points). For our purpose the important feature is that there is a family of planes¹ normal to a set of one of more unit vectors² k_i, separated by a distances d_i. It is usual to call each family of parallel planes a crystallographic plane, thought of as representative of the equivalence class. When a parallel beam of monochromatic x-rays, with wavelength λ, in direction ξ are incident on the crystal such that $\theta :=-{{{\rm sin} }^{-1}}\xi \cdot {{k}_{i}}$ satisfies Bragg's law (see figure 3),

$$\begin{eqnarray}&&2{{d}_{i}}{\rm sin} \theta =n\lambda ,\end{eqnarray} \tag{ 7 }$$

for some $n\in \mathbb{N}$, part of the beam is diffracted and continues in the direction η in the plane defined by k_i and ξ and at an angle $2\theta$ to ξ. The integer n is called the order of diffraction and we will denote $-{{{\rm sin} }^{-1}}(n\lambda /2{{d}_{i}})$ by ${{\theta }_{i,d}}$. In vectors

$$\begin{eqnarray}&&\eta =-\left( {{k}_{i}}\cdot \xi \right){{k}_{i}}+{{k}_{i}}\times \left( \xi \times {{k}_{i}} \right).\end{eqnarray} \tag{ 8 }$$

Note that k_i is defined only up to sign but the choice of sign does not effect η.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In practice one might expect to observe diffraction from between one and three orders, and $n\leqslant 3$. Bragg diffraction is thought of as elastic scattering in that the diffracted x-rays have the same wavelength (and hence energy) as the incident x-rays. It is also described as coherent scattering as phase relationships are maintained on the length scale of the distance between planes. We will ignore any attenuation of the diffracted x-rays by the material.

A polycrystalline material is understood to be a bounded region $D\subset {{\mathbb{R}}^{3}}$ and a finite set of subsets ${{D}_{\alpha }}$ (crystals) each with a positive measure and a Lipschitz boundary such ${{\bigcup }_{\alpha }}{{D}_{\alpha }}=D$ and their pair-wise intersection is measure zero. Each of the ${{D}_{\alpha }}$ is a crystal with the same set of planes and distances between planes except that the normals to the planes are rotated by some element of the rotation group ${{R}_{\alpha }}\in SO(3)$ relative to some reference configuration. Metals form typical examples of such polycrystalline materials and in this case the ${{D}_{\alpha }}$ are called grains and we are assuming the grains all have the same crystal structure and there are no gaps between them. We note also that is an idealisation. For example real crystals can have defects and impurities where the periodic structure breaks down, they might have residual stress that distorts the crystal structure, and except at absolute zero temperature the atoms in the crystal vibrate about their equilibrium positions.

We now describe the effect of x-ray diffraction on our sample of polycrystalline material. We will assume that the x-rays incident on D form a narrow beam (for example cylindrical) centred on some ray. This is possible for example using a synchrotron source. We will assume that the width of the beam is sufficient that its cross section includes a large number of crystals across the width of the beam. On the other hand the crystals themselves consist of sufficiently many planes for well defined Bragg diffraction from one or more planes. We are of course considering the beam to consist of plane waves propagating in a given direction and yet confined to a beam (collimated). This standard ray optics approximation is justified as the wavelength of x-rays is similar to the spacing of the atoms, and we have many crystals, let alone atoms, in the width of the beam.

Let us now assume that we are observing the transmitted and diffracted x-rays where they are incident on a screen normal to the incident direction ξ, at a distance large compared to the size of D. Although clearly it is important for constructive interference between x-rays on the scale of the distance between crystallographic planes, x-rays from our synchrotron will not be coherent over the path length between the sample and detector. Consequently the detector will measure the sum of the squared magnitude of all the x-ray waves incident at that point, or more realistically over a small square pixel containing that point.

For a beam along a given ray $x+t\xi$ the diffracted rays when observed at a distance L on a plane normal to the ray lie on circle centred on the projection of x on to this plane with a radius $L{\rm sin} 2{{\theta }_{i,n}}$, for the i the crystallographic plane and the n the order diffraction. The intensity of the observed ray in a circle $\{y\in {{\xi }^{\bot }}\left| \;|y|=L{\rm sin} 2{{\theta }_{i,n}} \right.\}$ is the measure of the intersection of the beam and the crystals ${{D}_{\alpha }}$ such that ${{R}_{\alpha }}{{k}_{i}}$ results in a diffracted beam at y on the plane, that is

$$\begin{eqnarray}&&y=L\left( {\rm sec} 2{{\theta }_{i,n}}\eta -\xi \right),\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&{{\eta }_{\alpha }}=-\left( {{R}_{\alpha }}{{k}_{i}}\cdot \xi \right){{R}_{\alpha }}{{k}_{i}}+{{R}_{\alpha }}{{k}_{i}}\times \left( \xi \times {{R}_{\alpha }}{{k}_{i}} \right).\end{eqnarray} \tag{ 10 }$$

Of course as we have described the situation this will result in a delta function measure for each of the finite number of crystals, but in practice for sufficiently fine crystals with orientations uniformly distributed in $SO(3)$ what is observed are circles of uniform intensity around the ring, and we must regard this as a limiting case of infinitely many crystals. Also in practice one does not observe a delta measure on the circle but some distribution of intensities with the maximum on the circles. This is due to the various approximations we have made including the distance L being finite.

We now consider the effect of elastic deformation on the material. Mathematically this means that we have a diffeomorphism $F:D\to F(D)\subset {{\mathbb{R}}^{3}}$ that takes the unstrained to the strained material. Of course the intersections of the crystallographyic planes are also preserved by this deformation. Under the usual small strain assumption we write $F(x)=x+U(x)$ for some small vector field U. We will assume that the crystals are sufficiently small that $\nabla U$ can be taken to be constant over a crystal for the purposes of studying the diffraction.

We now consider the effect of the deformation on the separation d of the planes normal to k in one crystal. Let $d^{\prime} k^{\prime} =DFdk$, where $DF={\bf I}+\nabla U$ is the derivative of F, be the deformed separation vector we have

$$\begin{eqnarray}&&d{{^{\prime} }^{2}}={{d}^{2}}(1+{{k}^{T}}(\nabla {{U}^{T}}+\nabla U)k+O(|\nabla U{{|}^{2}})).\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&d^{\prime} =d\sqrt{1+{{k}^{T}}2\epsilon k+O(|\nabla U{{|}^{2}})}\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&=\;d(1+{{k}^{T}}\epsilon k)+O(|\nabla U{{|}^{2}}),\end{eqnarray} \tag{ 13 }$$

where $\epsilon =(\nabla {{U}^{T}}+\nabla U)/2$ is the infinitesimal strain tensor. Hence we have the relative change in d

$$\begin{eqnarray}&&\frac{d^{\prime} -d}{d}={{k}^{T}}\epsilon k+O(|\nabla U{{|}^{2}}).\end{eqnarray} \tag{ 14 }$$

By contrast $k^{\prime} \cdot k=1+{{k}^{T}}\nabla {{U}^{T}}k+O(|\nabla U{{|}^{2}})$ so to the leading order term the direction is unchanged. We see that the circles in the diffraction pattern of the unstrained case, Debye–Scherrer rings, are distorted by the strain. Let $y=|y|\hat{y}$ with $\hat{y}$ the unit vector in ${{\xi }^{\bot }}$ in the plane of k and ξ then

$$\begin{eqnarray}&&|y|=L{\rm tan} 2\theta ^{\prime} ,\;{\rm where}\;\theta ^{\prime} ={{{\rm csc} }^{-1}}\frac{d(1+{{k}^{T}}\epsilon k)}{n\lambda }.\end{eqnarray} \tag{ 15 }$$

Often for high energy, small λ, x-rays (called hard x-rays) and low orders n the Bragg angles $2\theta$ are small, perhaps two or three degrees. On the other hand it is possible in a synchrotron beam to have the the detector screen many meters away. This means that while the Debye–Scherrer rings can still be measured the normal to the plane k_i, giving rise to the diffraction can be considered approximately parallel to the screen. This gives rise to a the deformed Debye–Scherrer ring being simply an ellipse defined by the points $y\in {{\xi }^{\bot }}$ such that

$$\begin{eqnarray}&&|y|=\frac{Ln\lambda }{(1+{{k}^{T}}\epsilon k)}+O({{\theta }^{2}}),\end{eqnarray} \tag{ 16 }$$

or with the small strain approximation.

$$\begin{eqnarray}&&|y|=Ln\lambda (1-{{k}^{T}}\epsilon k)+O({{\theta }^{2}},|\nabla U{{|}^{2}})\end{eqnarray} \tag{ 17 }$$

and using the small θ approximation in k we have approximatly

$$\begin{eqnarray}&&{{y}^{T}}({\bf I}-{{P}_{\xi }}\epsilon )y=R_{i,n}^{2},\end{eqnarray} \tag{ 18 }$$

where ${{R}_{i,n}}=nL\lambda /{{d}_{i}}$ is the radius of the unstrained Debye–Scherrer ring for the crystallographic plane k_i and order n. This describes the situation for uniform strain. The diffraction pattern observed will be the superposition of (slightly blurred) ellipses, and any given point on the detector plane will include x-rays diffracted from points along the line with different strains that meet the Bragg condition.

We will consider the case for small θ, in which the diffraction pattern contains only information about the strain transverse to ξ, ${{P}_{\xi }}\epsilon$. We might consider this a case of rich tomography in that for each ray we have not a scalar or a vector value but a function of two variables. One approach would be simply to solve the linear equations relating the strains in a grid of voxels to the intensities measured on the screens. Even for only one rotation axis, for similar pixel and voxel dimensions, we would have far more equations than variables and solving such a system would be inefficient. Our aim therefore is reducing the redundancy of the data while still retaining sufficient information for a unique reconstruction.

Given the wealth of data in a complete diffraction pattern it might be supposed that we could at least reconstruct the distribution of strains along the x-ray beam. A typical diffraction pattern from a small area with uniform stress consists of a series concentric ellipses for different orders. For a complete beam path we have the sum of such ellipses weighted by the prevalence of each transverse strain along the beam path. We might suspect that we could fit a weighted sum of ellipses to this pattern and then interpret the weights of the proportion of each transverse strain tensor present along the ray. We will show however that this is not the case as different distributions of ellipses can produce the same image.

An ellipse centred at the origin in the plane can in general be represented as a symmetric 2 × 2 matrix A with the points on the ellipse satisfying ${{y}^{T}}Ay-1=0$. In our case A(s) is the matrix

$$\begin{eqnarray}&&\frac{{{d}^{2}}}{n{{\lambda }^{2}}{{L}^{2}}}{{P}_{\xi }}\left( {\bf I}-\varepsilon (x+s\xi ) \right),\end{eqnarray} \tag{ 19 }$$

in a suitable basis on ${{\xi }^{\bot }}$

The unit density on the ellipse is therefore the delta measure

$$\begin{eqnarray}&&{{\delta }_{\{{{y}^{\prime }}:{{y}^{\prime T}}A{{y}^{\prime }}-1=0\}}}(y)=\frac{1}{|2Ay|}\delta ({{y}^{T}}Ay-1),\end{eqnarray} \tag{ 20 }$$

we use the notation ${{\delta }_{S}}$ for the unit measure on a set S and just δ for the standard Dirac delta measure at the origin in one dimensional space, that is an abbreviation for ${{\delta }_{\{0\}}}$. Ignoring the spread of the of the Debye–Scherrer ring our diffraction pattern is the sum of such densities as A runs over the A(s) be represented (in this idealized case) by

$$\begin{eqnarray}&&f(y)=\int _{-\infty }^{\infty }{{\delta }_{\{{{y}^{\prime }}:{{y}^{\prime T}}A(s){{y}^{\prime }}-1=0\}}}(y){\rm d}s,\end{eqnarray} \tag{ 21 }$$

where we have suppressed x and ξ as for the present discussion we regard them as fixed. We seek then information about the distribution g(A) of A that gave rise to this. In a very general sense (see for example [6, 12]) a Funk–Radon transform is the integral of a function on the plane over a two parameter family of curves. In this respect the two dimensional Radon transform (the integral over all lines) and its adjoint, also known as the backprojection operator, (the integral over translated sine waves of the same frequency over half a period) are both Funk–Radon transforms.

When straight lines occur in a photograph one way to detect them automatically is to take a filtered Radon transform, known as the Hough transform in computer vision, as lines densities in the image are mapped to delta functions by the radon transform. The Hough transform has been generalized to the detection of ellipses in a photograph [25] and one might expect that this would work for x-ray diffraction of polycrystalline materials. However the family of concentric ellipses is parameterized by three variables ${{a}_{11}},{{a}_{12}}$ and a₂₂, which satisfy ${{a}_{11}}\gt 0,{{a}_{22}}\gt 0$ and $a_{12}^{2}\lt {{a}_{11}}{{a}_{22}}$ to ensure the matrix is positive definite. Alternatively the semi-major axes and the rotation of the principal axis can be used as three parameters.

Consider polar coordinates ${\bf y}={{(r{\rm cos} \psi ,r{\rm sin} \psi )}^{T}}$ which transforms ${{y}^{T}}Ay-1=0$ to

$$\begin{eqnarray}&&{{a}_{11}}{{{\rm cos} }^{2}}\psi +{{a}_{12}}2{\rm sin} \psi {\rm cos} \psi +{{a}_{11}}{{{\rm sin} }^{2}}\psi =1/{{r}^{2}},\end{eqnarray} \tag{ 22 }$$

which is the equation of a plane in $({{a}_{11}},{{a}_{12}},{{a}_{22}})$ space normal to the vector $K(\psi )=({{{\rm cos} }^{2}}\psi ,2{\rm sin} \psi {\rm cos} \psi ,{{{\rm sin} }^{2}}\psi )$ and a distance $1/{{r}^{2}}|K(\psi )|$ from the origin. As r varies for fixed ψ this gives us a family of parallel planes, while as ψ varies the $K(\psi )$ describes a circle in the plane ${{a}_{11}}+{{a}_{22}}=1$. Our diffraction data is thus a two parameter family of plane integrals and is insufficient to determine a general density in the three dimensional $({{a}_{11}},{{a}_{12}},{{a}_{22}})$ space. If the set of possible densities were sufficiently restricted it may well be possible to recover it from integrals over a two parameter family of curves. We could assume that it is concentrated on a (possibly self intersecting) curve as in our case it is parameterized by s, but as this curve is unknown this is of no obvious help. See figure 4.

Zoom In Zoom Out Reset image size

The case in which ${{a}_{11}}+{{a}_{22}}=2t$ for t a known constant has a unique solution as our planes intersect this plane in all possible lines, in fact this is exactly the two dimensional Radon transform in the variables $(s,{{a}_{12}})$ where $s=({{a}_{11}}-{{a}_{22}})/2$. However the determination of the distribution of the perpendicular strain along a ray from one diffraction pattern (for one energy of x-ray) is clearly too ambitious in general. Indeed, even in special cases where it is possible it only gives the information on which plane strains tensors arise and how frequently, not where they are on the line.

Instead we set a more modest aim: to determine the centre of mass in of the distribution in $({{a}_{11}},{{a}_{12}},{{a}_{22}})$ space. This we can do from only three choices of direction $K(\psi )$ as we will show.

Korsunsky [8] assumes that the average strain along a ray in a given direction normal to the ray is proportional the deflection of the diffraction pattern and the intuition behind this is sound, but as the diffraction pattern is an average of elliptical Debye–Scherrer rings and each of those rings is itself blurred it is not clear what average value for deflection should be taken. In this section we show that a certain moment of the diffraction pattern in the required direction is the correct choice.

We have in general

$$\begin{eqnarray}&&\mathcal{I}^{\prime} (y)={{\int }_{\{A:{{y}^{T}}Ay=1\}}}g(A)\;\sqrt{2}{\rm d}{{a}_{11}}{\rm d}{{a}_{22}}{\rm d}{{a}_{12}}.\end{eqnarray} \tag{ 23 }$$

We have used the measure ${\rm d}A=\sqrt{2}{\rm d}{{a}_{11}}{\rm d}{{a}_{22}}{\rm d}{{a}_{12}}$ as it is the measure associated with the Fröbenius norm on ${{S}^{2}}{{\mathbb{R}}^{2}}$ and hence invariant under the action of $SO(2)$.

Without loss of generality we now take $y={{y}_{1}}{{e}_{1}}$ so that ${{a}_{11}}=1/y_{1}^{2}$ and (23) becomes

$$\begin{eqnarray}&&\mathcal{I}^{\prime} ({{y}_{1}}{{e}_{1}})={{\int }_{\{A:{{a}_{11}}=1/y_{1}^{2}\}}}g(A)\;\sqrt{2}{\rm d}{{a}_{11}}{\rm d}{{a}_{22}}{\rm d}{{a}_{12}}\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}=\sqrt{2}\int g(\left( \begin{array}{ccccccccccccccc} 1/y_{1}^{2} & {{a}_{12}} \\ {{a}_{12}} & {{a}_{22}} \\ \end{array} \right)){\rm d}{{a}_{22}}{\rm d}{{a}_{12}}.\end{eqnarray} \tag{ 25 }$$

We now define the moment of $\mathcal{I}^{\prime}$ in the direction of any y as

$$\begin{eqnarray}&&\mathcal{M}(y)=\mathop{\int }\limits_{r=0}^{\infty }r\mathcal{I}^{\prime} ({{r}^{-1/2}}{\mkern 1mu} y)\;{\rm d}r.\end{eqnarray} \tag{ 26 }$$

As the substitution ${{y}_{1}}={{r}^{-1/2}}$ in (24) gives ${{a}_{11}}=r$.

$$\begin{eqnarray}\mathcal{M}({{e}_{1}})=\int _{r=0}^{\infty }r\int g(\left( \begin{array}{ccccccccccccccc} r & {{a}_{12}} \\ {{a}_{12}} & {{a}_{22}} \\ \end{array} \right)){\rm d}{{a}_{22}}{\rm d}{{a}_{12}}{\rm d}r=\int \int \int {{a}_{11}}g(A){\rm d}{{a}_{22}}{\rm d}{{a}_{12}}{\rm d}{{a}_{11}},\end{eqnarray} \tag{ 27 }$$

noting that as A is positive definite ${{a}_{11}}\gt 0$. Applying a rotation of coordinates we see that for a general unit vector y

$$\begin{eqnarray}&&\mathcal{M}(y)=\int {{y}^{T}}Ay\;g(A){\rm d}A,\end{eqnarray} \tag{ 28 }$$

in particular (27) applies with 2 in place of 1.

The defining property of the density g guarantees

$$\begin{eqnarray}&&\int _{-\infty }^{\infty }{{a}_{ij}}(s){\rm d}s={{\int }_{{{S}^{2}}{{\mathbb{R}}^{2}}}}{{a}_{ij}}g(A){\rm d}A,\end{eqnarray} \tag{ 29 }$$

so we have from these three moment calculations along radial directions of the diffraction pattern the transverse ray transform $J\varepsilon (x,\xi )$.

For the transverse ray transform J a simple explicit reconstruction method is suggested by [18]. Let ζ be any unit vector and consider any $\xi \in {{\zeta }^{\bot }}$. We then have

$$\begin{eqnarray}&&\zeta \cdot J\varepsilon (x,\xi )\zeta =\int _{-\infty }^{\infty }{{\zeta }^{T}}{{P}_{\xi }}\left( \varepsilon (x+s\xi ) \right)\zeta \;{\rm d}s=\int _{-\infty }^{\infty }{{\zeta }^{T}}\varepsilon (x,\xi )\zeta \;{\rm d}s,\end{eqnarray} \tag{ 30 }$$

as the projection leaves ζ unchanged. We note this is exactly the two dimensional Radon transform of the scalar ${{\zeta }^{T}}\varepsilon \zeta$ in the plane through x normal to ζ. Taking ζ as each of the unit basis vectors e_i in turn we can recover the diagonal terms ${{\varepsilon }_{ii}}=e_{i}^{T}\varepsilon {{e}_{i}}$ by Radon transform inversion plane by plane. Experimentally we would achieve this by rotating the object half a turn about each the coordinate axes. To obtain the off diagonal terms (again a special case of the polarization identity(5) note that the choice $\zeta =({{e}_{1}}+{{e}_{2}})/\sqrt{2}$ gives

$$\begin{eqnarray}&&{{\zeta }^{T}}\varepsilon \zeta =(\frac{1}{2}){{\varepsilon }_{11}}+(\frac{1}{2}){{\varepsilon }_{22}}+{{\varepsilon }_{12}},\end{eqnarray} \tag{ 31 }$$

hence we obtain the off diagonal term ${{\varepsilon }_{12}}$ by a further rotation about this tilted axis, and the other off diagonal terms are obtained similarly. In general any six vectors ${{\zeta }_{i}}$ chosen according to the criterion in section 3 will yield sufficient data. We now have

Theorem 7.1. 

• 1.  
  For any plane π normal to a unit vector ζ the components of ${{\zeta }^{T}}\varepsilon \zeta {{|}_{\pi }}$ can be obtained as the inverse Radon transform of the moment data $\mathcal{M}(\zeta )$ for the diffraction patterns for the rays in π.
• 2.  
  Given six unit vectors ${{\xi }_{i}}$ that do not all lie on the same projective conic, the moment data $\mathcal{M}(\zeta )$ for all the rays $x+t\xi$ with $\xi \in {{\zeta }^{\bot }}$ and with non empty intersection with D determine $\zeta _{i}^{T}\varepsilon {{\zeta }_{i}}$ and hence six linear equations for known $\zeta _{i}^{T}\varepsilon {{\zeta }_{i}}$ can be solved to give everywhere in D.

While this method gives a proof sufficiency of data, one suspects that a reconstruction would be possible with less than six rotation axes. Indeed for each ray the above algorithm uses a moment of the diffraction pattern only in the direction of the rotation, a single measurement of a one-dimensional section of the diffraction pattern. It would be desirable for practical purposes to have an explicit reconstruction algorithm that uses rotation about fewer axes, but makes better use of the data for each ray. We do not currently have such an algorithm, but of course one could take a brute force approach by solving a very large sparse system by an iterative method from numerical linear algebra. At the other extreme, in the next section, we give a reconstruction formula using data from all rays.

For the scalar x-ray transform there is an explicit inversion formula [11] for data for all rays passing through the support of the function to be recovered. This formula is a generalization of the filtered back projection formula for the two dimensional case, and it has several uses even though it is not usually practical to measure over a set of rays forming a reasonably uniform sample of all the four dimensional space of possible rays. For example if one can measure an open set in the space of rays, by tilting a sample about two axes in a parallel beam with a planar detector, one can show that the inverse ray transform has a unique solution using the Paley–Weiner theorem and analytic continuation. Of course this reconstruction is unstable in any Sobolev spaces [7] but it can be analysed using microlocal methods to understand the artefacts that result from such a truncated data reconstruction. In anticipation that similar limited data may be available for x-ray diffraction we derive an explicit filtered back projection inversion formula for the transverse ray transform, following closely the method used in [18] for the Truncated Transverse Ray transform. This section and the one that follows uses several concepts from the general theory in that book and is not needed for the understanding of the simple reconstruction procedure detailed in the previous section.

First some essential notation. We will denote by $B:{{C}^{\infty }}({{\mathbb{R}}^{3}}\times \mathbb{R}_{0}^{3};{{S}^{2}}{{\mathbb{R}}^{3}})\to C_{0}^{\infty }({{\mathbb{R}}^{3}};{{S}^{2}}{{\mathbb{R}}^{3}})$ the backprojection operator

$$\begin{eqnarray}&&B\phi \;(x)=\frac{1}{4\pi }{{\int }_{\Omega }}\phi (x,\xi ){\rm d}\omega (\xi ),\end{eqnarray} \tag{ 32 }$$

where Ω is the unit sphere and ω its standard measure. On $C_{0}^{\infty }({{\mathbb{R}}^{3}};{{S}^{k}}{{\mathbb{R}}^{3}})$, the space of smooth symmetric rank k tensors fields we denote by $\delta :C_{0}^{\infty }({{\mathbb{R}}^{3}};{{S}^{k}}{{\mathbb{R}}^{3}})\to C_{0}^{\infty }({{\mathbb{R}}^{3}};{{S}^{k-1}}{{\mathbb{R}}^{3}})$ the generalization of divergence

$$\begin{eqnarray}&&{{(\delta f)}_{{{i}_{1}}...{{i}_{k-1}}}}=\frac{\partial }{\partial {{x}_{j}}}{{f}_{{{i}_{1}}...{{i}_{k-1}}j}}.\end{eqnarray} \tag{ 33 }$$

(we use the convention of summation over repeated indices) and the negative of its formal adjoint, the symmetric derivative $d:C_{0}^{\infty }({{\mathbb{R}}^{3}};{{S}^{k}}{{\mathbb{R}}^{3}})\to C_{0}^{\infty }({{\mathbb{R}}^{3}};{{S}^{k+1}}{{\mathbb{R}}^{3}})$

$$\begin{eqnarray}&&{{(df)}_{{{i}_{1}}...{{i}_{k+1}}}}=\sigma \left( \frac{\partial }{\partial {{x}_{{{i}_{k+1}}}}}{{f}_{{{i}_{1}}...{{i}_{k}}}} \right),\end{eqnarray} \tag{ 34 }$$

where σ is the symmetrization operator on tensors. We will denote by Δ the Laplacian taken on components. The inverse ${{\Delta }^{-1}}$ is defined by convolution with the fundamental solution. The trace of a rank 2 tensor field f is denoted by ${\rm Tr}f={{f}_{ii}}$.

We are now ready to state explicit inversion formulae using complete data for the transverse ray transform

Theorem 8.1. Let $f\in C_{0}^{\infty }({{S}^{2}})$ be a symmetric tensor field then

$$\begin{eqnarray}&&f=\frac{4}{\pi }\left( 4{\bf I}-{\bf I}{\rm Tr}+4{{\Delta }^{-1}}d\delta -{{\Delta }^{-1}}{{d}^{2}}{\rm Tr}-{{\Delta }^{-1}}{\bf I}{{\delta }^{2}} \right){{(-\Delta )}^{1/2}}BJf.\end{eqnarray} \tag{ 35 }$$

We begin with some formulae for the Fourier transforms of terms that arise in the operators BJ which are simplifications of those given in [Sec 2.11, 18] for the specific cases we need. For m = 0,1, or 2 define ${{\Phi }^{m}}:C_{0}^{\infty }({{\mathbb{R}}^{3}};{{S}^{m}}{{\mathbb{R}}^{3}})\to C_{0}^{\infty }({{\mathbb{R}}^{3}};{{S}^{m}}{{\mathbb{R}}^{3}})$ by

$$\begin{eqnarray}&&{{({{\Phi }^{m}}f)}_{{{i}_{1}},{{i}_{m}}}}=\frac{1}{4\pi }{{\int }_{\Omega }}{{\xi }_{{{i}_{1}}}}{{\xi }_{{{i}_{m}}}}\int _{-\infty }^{\infty }{{f}_{{{j}_{1}}{{j}_{m}}}}(x+\xi t){{\xi }_{{{j}_{1}}}}{{\xi }_{{{j}_{m}}}}\;{\rm d}t{\mkern 1mu} {\rm d}\omega (\xi ).\end{eqnarray} \tag{ 36 }$$

Using symmetry in t and the substition $x^{\prime} =x+t\xi$ we obtain

$$\begin{eqnarray}&&{{({{\Phi }^{m}}f)}_{{{i}_{1}},{{i}_{m}}}}=\frac{1}{2\pi }{{\int }_{\Omega }}{{\xi }_{{{i}_{1}}}}{{\xi }_{{{i}_{m}}}}\int _{0}^{\infty }{{f}_{{{j}_{1}}{{j}_{m}}}}(x+\xi t){{\xi }_{{{j}_{1}}}}{{\xi }_{{{j}_{m}}}}\;{\rm d}t{\mkern 1mu} {\rm d}\omega (\xi )\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&=\frac{1}{2\pi }{{\int }_{\Omega }}\int _{0}^{\infty }\frac{{{(x^{\prime} -x)}_{{{i}_{1}}}}{{(x^{\prime} -x)}_{{{i}_{m}}}}}{|x^{\prime} -x{{|}^{m}}}{{f}_{{{j}_{1}}{{j}_{m}}}}(x^{\prime} )\frac{{{(x^{\prime} -x)}_{{{j}_{1}}}}{{(x^{\prime} -x)}_{{{j}_{m}}}}}{|x^{\prime} -x{{|}^{m+2}}}\;{{t}^{2}}\;{\rm d}t{\mkern 1mu} {\rm d}\omega (\xi )\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&=\frac{1}{2\pi }{{f}_{{{j}_{1}}{{j}_{m}}}}*\frac{{{x}_{{{i}_{1}}}}{{x}_{{{i}_{m}}}}{{x}_{{{j}_{1}}}}{{x}_{{{j}_{m}}}}}{|x{{|}^{2m+2}}}.\end{eqnarray} \tag{ 39 }$$

For m = 0 this is just the classical result giving the inverse of the scalar ray transform in three dimensions: $f={{(-\Delta )}^{1/2}}BXf$. For example in [11] equation (2.2) on p19 with $\alpha =1$, indeed we have immediately

$$\begin{eqnarray}&&\widehat{{{\Phi }^{0}}f}(y)=\hat{f}(y)\frac{\pi }{|y|},\end{eqnarray} \tag{ 40 }$$

and using Eq (2.11.6) of [18] we have also for f a vector field

$$\begin{eqnarray}&&\widehat{{{\Phi }^{1}}f}(y)={{\Pi }_{y}}\hat{f}(y)\frac{\pi }{2|y|}.\end{eqnarray} \tag{ 41 }$$

For m = 2 we need the linear operator

$$\begin{eqnarray}&&\tilde{{{\Pi }_{y}}}(f)=\frac{2}{3}{{\Pi }_{y}}f(y){{\Pi }_{y}}+\frac{1}{3}{\rm Tr}({{\Pi }_{y}}f){{\Pi }_{y}},\end{eqnarray} \tag{ 42 }$$

which is $\sigma ({{\Pi }_{y}}\otimes {{\Pi }_{y}})$ contracted with f over two indices. To see this notice that in components this is a sum of terms of the form ${{\Pi }_{{{i}_{\pi (1)}}{{i}_{\pi (2)}}}}{{\Pi }_{{{i}_{\pi (3)}}{{i}_{\pi (4)}}}}{{f}_{{{i}_{1}}{{i}_{2}}}}$ over all permutations π of four symbols. Of the 24 permutations 16 give indecies i₁ and i₂ on separate factors Π whereas 8 result in i₁ and i₂ on the same factor. In [18] $\tilde{{{\Pi }_{y}}}(f)$ is denoted by ${{\varepsilon }^{2}}(y)/f$. We have then

$$\begin{eqnarray}&&\widehat{{{\Phi }^{2}}f}(y)={{\tilde{\Pi }}_{y}}(\hat{f}(y))\frac{3\pi }{8|y|}.\end{eqnarray} \tag{ 43 }$$

To apply these results to the transverse ray transform first we expand equation 3 as

$$\begin{eqnarray}&&{{P}_{\xi }}{\mkern 1mu} f=f-\frac{\xi {{\xi }^{T}}f+f\xi {{\xi }^{T}}}{|\xi {{|}^{2}}}+\frac{({{\xi }^{T}}f\xi )\xi {{\xi }^{T}}}{|\xi {{|}^{4}}}.\end{eqnarray} \tag{ 44 }$$

Now we write

$$\begin{eqnarray}&&BJf={{G}_{1}}{\mkern 1mu} f+{{G}_{2}}{\mkern 1mu} f+{{G}_{3}}{\mkern 1mu} f,\end{eqnarray} \tag{ 45 }$$

corresponding to the three terms in equation 44. We notice G₁ is simply BX and hence

$$\begin{eqnarray}&&{{G}_{1}}{\mkern 1mu} f=\frac{1}{2\pi }f*|x{{|}^{-2}},\end{eqnarray} \tag{ 46 }$$

and the third term is exactly ${{\Phi }^{2}}f$ so

$$\begin{eqnarray}&&{{({{G}_{3}}{\mkern 1mu} f)}_{ij}}=\frac{1}{2\pi }{{f}_{kl}}*\frac{{{x}_{i}}{{x}_{j}}{{x}_{k}}{{x}_{l}}}{|x{{|}^{6}}}.\end{eqnarray} \tag{ 47 }$$

For the second term one must look again at equation 37 for m = 1, but applied to ${{f}_{ip}}{{\xi }_{p}}$. We find

$$\begin{eqnarray}&&{{({{G}_{2}}{\mkern 1mu} f)}_{ij}}=\frac{1}{2\pi }\left( {{f}_{ik}}*\frac{{{x}_{j}}{{x}_{k}}}{|x{{|}^{4}}}+{{f}_{jk}}*\frac{{{x}_{i}}{{x}_{k}}}{|x{{|}^{4}}} \right).\end{eqnarray} \tag{ 48 }$$

Taking the Fourier transforms of these expressions gives us the following lemma:

Lemma 9.1. Let $g=(2/\pi )BJf$ then

$$\begin{eqnarray}&&|y|\hat{g}=\hat{f}-\frac{1}{2}\left( {{\Pi }_{y}}\hat{f}+\hat{f}{{\Pi }_{y}} \right)+\frac{3}{8}{{\tilde{\Pi }}_{y}}(\hat{f}).\end{eqnarray} \tag{ 49 }$$

$z=y/|y|$

${{\Pi }_{y}}={\bf I}-z{{z}^{T}}$

$$\begin{eqnarray}&&{{\tilde{\Pi }}_{y}}f=\frac{2}{3}f+\frac{1}{3}{\rm Tr}f\;{\bf I}-\frac{2}{3}(z{{z}^{T}}f+fz{{z}^{T}})+({{z}^{T}}fz)z{{z}^{T}}-\frac{1}{3}{\rm Tr}fz{{z}^{T}}-\frac{1}{3}({{z}^{T}}fz){\bf I}.\end{eqnarray} \tag{ 50 }$$

Substitution of (50) in (49) gives

$$\begin{eqnarray}&&|y|\hat{g}=\frac{1}{4}\hat{f}+\frac{1}{4}(z{{z}^{T}}\hat{f}+\hat{f}z{{z}^{T}})+\frac{3}{8}({{z}^{T}}\hat{f}z)z{{z}^{T}}+\frac{1}{8}{\bf I}{\rm Tr}\hat{f}-\frac{1}{8}z{{z}^{T}}{\rm Tr}\hat{f}-\frac{1}{8}{{z}^{T}}\hat{f}z{\bf I}.\end{eqnarray} \tag{ 51 }$$

We now verify by direct substitution that

$$\begin{eqnarray}&&\hat{f}=\frac{1}{|y|}\left( 4\hat{g}-2(z{{z}^{T}}\hat{g}+\hat{g}z{{z}^{T}})-{\bf I}{\rm Tr}\hat{g}-z{{z}^{T}}{\rm Tr}\hat{g}+({{z}^{T}}\hat{g}z){\bf I} \right).\end{eqnarray} \tag{ 52 }$$

theorem 8.1 now follows from the definitions of d and δ. Such calculations are simplified by the observation that the algebra generated by the five terms in the brackets in (35) is finite dimensional over the field generated by δ. This means for elements that are units one has to find only five coefficients to compute an inverse.

We note that the pseudo differential operator in theorem 8.1 is elliptic as it has an explicit (left and right) inverse for all $y\ne 0$. This means that micro-local methods, such as Lambda Tomography, may be used to recover some of the singular support of a non-smooth from limited data in a straightforward extension of the method used for the scalar x-ray transform [13]. Indeed the operator ${{(-\Delta )}^{1/2}}B$ is simply the inverse of the scalar x-ray transform applied to each component. This means that any limited data inversion method for the x-ray transform, for example methods for lines passing through a curve external to the support, yields a reconstruction algorithm for J.

In the approach we have described so far for each ray the data observed is the diffraction pattern. Another approach, suggested for neutron strain tomography, is apply polychromatic beam containing an interval of wavelengths λ at known intensities and use a single detector for that ray that can resolve the spectrum–the energy transmitted at each frequency. The decrease in the transmitted energy is due to the diffraction of the neutrons. For a fixed ray, wavelength, order and crystallographic plane some neutrons will be scattered when there is a θ in $(-\pi ,\pi )$ such that $2d^{\prime} {\rm sin} \theta =n\lambda$ given that there is a sufficient range of distorted distances d'. What is observed [17] is that the transmitted energy against λ curve has a discontinuity, a sudden increase, as theta approaches $\pi /2$, so the neutrons are 'back scattered', and no further diffraction happens from that plane and order. This discontinuity is termed the Bragg edge. The treatment in the literature [2, 15–17, 19] is quite heuristic and the various curve fitting procedures to quantify the mean of $(d^{\prime} -d)/d$ along the line $x+t\xi$. However for our purposes what matters is that the measurements give only

$$\begin{eqnarray}&&\int _{-\infty }^{\infty }{{\xi }^{T}}\varepsilon (x+s\xi )\xi {\rm d}s=I\varepsilon (x,\xi ).\end{eqnarray} \tag{ 53 }$$

An important distinction between J and I is that even for complete data I has a null space consisting of f such that f = dh for some smooth vector field h. We will need to treat the case where f is zero outside D so has a discontinuity on $\partial D$ so we will give an explicit argument. In particular we note that in the interior $\varepsilon =(1/2){\rm d}U$. For simplicity take $\xi ={{e}_{3}}$ and assume D is convex

$$\begin{eqnarray}&&\int _{-\infty }^{\infty }e_{3}^{T}\varepsilon (x+s{{e}_{3}}){{e}_{3}}{\rm d}s=\int _{-\infty }^{\infty }{{\varepsilon }_{3}}(x+s{{e}_{3}}){\rm d}s,\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&=\int _{-\infty }^{\infty }\partial {{U}_{3}}/\partial {{x}_{3}}{\rm d}{{x}_{3}},\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&={{U}_{3}}(x_{3}^{+})-{{U}_{3}}(x_{3}^{-}),\end{eqnarray} \tag{ 56 }$$

where $x_{3}^{+}$ and $x_{3}^{-}$ are the x₃ coordinates of the intersection of the ray with $\partial D$. We see that such a technique is only sensitive to the overall change in thickness of the sample along the ray and not the distribution of strains. The restriction to a convex body is of course unnecessary, we would simply get the sum of the thickness change in each component of the intersection of the ray with D. In the case of a void within an object with known boundary this may well be useful as the displacement of the void could be be measured more accurately using Bragg edge than conventional tomographic methods.

In conclusion we see that it is indeed possible to obtain the component of the strain normal to a plane from observing the diffraction due to a narrow beams in that plane. We recommend that rather than using a curve fitting procedure the appropriate radial moment of the diffraction profile is used. If six well chosen rotation axes are used and for each axis such diffraction measurements are made for a set of lines in the plane, then all six components of the strain can be reconstructed at points on the intersection six of these planes. To reconstruct the stress throughout the sample this way requires a time consuming measurement process. It remains to be seen if there are explicit reconstruction procedures using more of the diffraction data for fewer rotation axes. Certainly it would be possible to construct a system of sparse linear equations for a discretization of the problem and attempt a solution using standard iterative algorithms.

We have also shown, contrary to previous assertions, it is not possible to simply reconstruct the internal strain distribution within a body from Bragg edge measurements without some form of extra information because it just gives the average fractional change in dimension (strain) along the ray path.

This work was partly funded by the EPSRC building global engagements grant MIRAN (EP/K00428X/1), the Tomography CCPi (EP/J010456/1), the Next Generation Multi-Dimensional x-ray Imaging grant (EP/M010619/1), the Residual Stress & Damage Characterization Extending Across. Length and Time Scales grant (EP/F028431/1), the Structural evolution across multiple time and length scales grant (EP I02249X/1) and the authors are grateful for the helpful discussions on this topic with Vladimir Sharafutdinov while visiting under the MIRAN programme. The authors would like to thank Alexander Korsunsky and Brian Abbey for helpful discussions and for feedback on earlier drafts, and the anonymous referees for their suggested improvements. WL would also like to thank the Otto Mønstead foundation for funding the visiting professorship at the Danish Technical University during which the paper was completed.