Diffraction tomography of strain

We consider whether it is possible to recover the three dimensional strain field tomographically from neutron and x-ray diffraction data for polycrystalline materials. We show that the distribution of strain transverse to a ray cannot be deduced from one diffraction pattern accumulated along that path, but that a certain moment of that data corresponds to the transverse ray transform of the strain tensor and so may be recovered by inverting that transform given sufficient data. We show that the whole strain tensor can be reconstructed from diffraction data measured using rotations about six directions that do not lie on a projective conic. In addition we give an inversion formula for complete data for the transverse ray transform. We also show that Bragg edge transmission data, which has been suggested for strain tomography with polychromatic data, cannot provide the strain distribution within the material but only the average along the ray path.


Introduction
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 ℝ 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 Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
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.
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][16][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.

Ray transforms of rank 2 tensor fields
2 3 be the space of smooth compactly supported functions on ℝ 3 with values in the six dimensional vector space ℝ S 2 3 of symmetric two-tensors on ℝ 3 . The simplest ray transform we can take is to apply the scalar ray transform to each component where ℝ = ℝ ⧹{0} . Without loss of generality, when describing the ray ξ + x t , we can take ξ to be a unit vector and ξ ∈ ⊥ x . Following Sharafutdinov [18] we also define the transverse ray transform T 2 is the projection on to the plane ξ ⊥ , we identify rank two tensors with matrices for the purposes of multiplication and I is the identity matrix.
We also define the longitudinal ray transform T 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].

Determination of a symmetric matrix from its quadratic form
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 S 2 3 determines a quadratic form = Q x x Ax ( ) T and it is well known that the off-diagonal elements i T i are independent as linear forms in the a ij , or equivalently that the outer products x x i i T span ℝ S 2 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) 6. This is precisely the same as saying that the six points on the projective plane P 2 lie on a projective conic defined by (non-zero scalar multiples of) A. Indeed any five distinct points in P 2 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, hence points on P 2 are simply directions in ℝ 3 . Projective lines in P 2 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]. One important special case is where the six directions x i lie on some antipodal pair of small circles defined by | | = x k a · for some fixed unit vector k and < < a 0 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

x-ray diffraction and strain
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 1 normal to a set of one of more unit vectors 2 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 θ ξ i for some ∈ ℕ 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 to ξ. The integer n is called the order of diffraction and 1 In crystallography the planes are indexed by triples of integers called Miller indices but we shall take them to have been enumerated by a single index. 2 There are many clashes of notation between crystallography and tensor tomography. In crystallography it is common for our k to be called g and our ξ to be called k Inverse Problems 31 (2015) 045005 W R B Lionheart and P J Withers In vectors Note that k i is defined only up to sign but the choice of sign does not effect η.
In practice one might expect to observe diffraction from between one and three orders, and ⩽ n 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.

R B Lionheart and P J Withers
A polycrystalline material is understood to be a bounded region ⊂ ℝ D 3 and a finite set of subsets α D (crystals) each with a positive measure and a Lipschitz boundary such ⋃ = α α D D and their pair-wise intersection is measure zero. Each of the α D 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 SO (3) relative to some reference configuration. Metals form typical examples of such polycrystalline materials and in this case the α D 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 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 sin 2 i n , , for the i the crystallographic plane and the n the order diffraction. The intensity of the observed ray in a circle is the measure of the intersection of the beam and the crystals α D such that α R k i results in a diffracted beam at y on the plane, that is 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 F D : ( ) 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 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.
is the derivative of F, be the deformed separation vector we have is the infinitesimal strain tensor. Hence we have the relative change in d T 2 Often for high energy, small λ, x-rays (called hard x-rays) and low orders n the Bragg angles θ 2 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 such that

By contrast
T 2 or with the small strain approximation.
and using the small θ approximation in k we have approximatly 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.

Extracting information from the diffraction pattern
We will consider the case for small θ, in which the diffraction pattern contains only information about the strain transverse to ξ, ϵ ξ P . 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 Ay 1 0 T . In our case A(s) is the matrix in a suitable basis on ξ ⊥ The unit density on the ellipse is therefore the delta measure y y Ay T : 1 0 T we use the notation δ 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 δ {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 y y A s y : ( ) 1 0 T 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 a , 11  which is the equation of a plane in a a a ( , , ) 11 12 22 space normal to the vector ψ ψ ψ ψ ψ = K ( ) ( cos , 2 sin cos , sin ) 2 2 and a distance ψ | | r K 1 ( ) 2 from the origin. As r varies for fixed ψ this gives us a family of parallel planes, while as ψ varies the ψ K ( ) describes a circle in the plane + = a a 1 11 22 . 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 a a ( , , ) 11 12 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.
The case in which + = a a t 2 11 22 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 ( , )2 11 22 . 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 a a ( , , ) 11 12 22 space. This we can do from only three choices of direction ψ K ( ) as we will show.  11 12 22 illustrating that the distribution of A values along one beam will cluster around a curve.

Average strain from moments
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 A  noting that as A is positive definite > a 0 11 . Applying a rotation of coordinates we see that for a general unit vector y T in particular (27) applies with 2 in place of 1. The defining property of the density g guarantees ij S ij 2 2 so we have from these three moment calculations along radial directions of the diffraction pattern the transverse ray transform ε ξ J x ( , ).

Measurement and reconstruction
For the transverse ray transform J a simple explicit reconstruction method is suggested by [18]. Let ζ be any unit vector and consider any ξ ζ ∈ ⊥ . We then have T T as the projection leaves ζ unchanged. We note this is exactly the two dimensional Radon transform of the scalar ζ εζ T 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 ε ε = e e ii i T 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 ζ = + e e ( ) 2 1 2 gives ζ εζ ε ε ε = + + 1 2 hence we obtain the off diagonal term ε 12 by a further rotation about this tilted axis, and the other off diagonal terms are obtained similarly. In general any six vectors ζ i chosen according to the criterion in section 3 will yield sufficient data. We now have 1. For any plane π normal to a unit vector ζ the components of ζ εζ | π T can be obtained as the inverse Radon transform of the moment data  ζ ( ) for the diffraction patterns for the rays in π. 2. Given six unit vectors ξ i that do not all lie on the same projective conic, the moment data  ζ ( ) for all the rays ξ + x t with ξ ζ ∈ ⊥ and with non empty intersection with D determine ζ εζ i T i and hence six linear equations for known ζ εζ i T 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.

Inversion of the transverse ray transform with complete data
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 where Ω is the unit sphere and ω its standard measure. On 3 , the space of smooth symmetric rank k tensors fields we denote by δ (we use the convention of summation over repeated indices) and the negative of its formal adjoint, the symmetric derivative where σ is the symmetrization operator on tensors. We will denote by Δ the Laplacian taken on components. The inverse Δ −1 is defined by convolution with the fundamental solution. The trace of a rank 2 tensor field f is denoted by = f f Tr ii . We are now ready to state explicit inversion formulae using complete data for the transverse ray transform Using symmetry in t and the substition ξ For m = 0 this is just the classical result giving the inverse of the scalar ray transform in three dimensions: 1 2 . For example in [11] equation (2.2) on p19 with α = 1, indeed we have immediately 0 and using Eq (2.11.6) of [18] we have also for f a vector field For m = 2 we need the linear operator (3) (4) 1 2 over all permutations π of four symbols. Of the 24 permutations 16 give indecies i 1 and i 2 on separate factors Π whereas 8 result in i 1 and i 2 on the same factor. In [18] where f is zero outside D so has a discontinuity on ∂D so we will give an explicit argument. In particular we note that in the interior ε = U (1 2)d . For simplicity take ξ = e 3 and assume D is convex where + x 3 and − x 3 are the x 3 coordinates of the intersection of the ray with ∂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.

Conclusions
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.
funding the visiting professorship at the Danish Technical University during which the paper was completed.