An explicit reconstruction algorithm for the transverse ray transform of a second rank tensor field from three axis data

We give an explicit plane-by-plane filtered back-projection reconstruction algorithm for the transverse ray transform of symmetric second rank tensor fields on Euclidean 3-space, using data from rotation about three orthogonal axes. We show that in the general case two axis data is insufficient but give an explicit reconstruction procedure for the potential case with two axis data


Introduction
The transverse ray transform (TRT) of rank two symmetric tensor fields in three dimensional Euclidean space has recently been shown to be of importance in x-ray diffraction strain tomography [5], however currently reconstruction algorithms are known only for data from rotations about six axes [8,Sec 5.1.6] or complete data [5]. Given that, in the proposed application, each projection is acquired laboriously using a raster scan, it is advantageous to perform the reconstruction using data from a minimum number of axes and making the most of the data collected from each projection. In this paper we give an explicit reconstruction formula for three axis data using similar techniques to those employed by [4] for the truncated transverse ray transform (TTRT). The inversion formula uses plane-by-plane filter and back-projection operations familiar from inversion of the parallel x-ray transform of a scalar field. We go on to show that data from only two rotation axes are insufficient in the general case, but give an explicit reconstruction for the potential case with two axis data. We present some numerical results for our reconstruction algorithms using simulated data.

Definitions and notation
We denote the complex vector space of symmetric R-bilinear maps R 3ˆR3 Ñ C by S 2 C 3 . The elements of this space are (complex-valued ) symmetric tensors of second rank on R 3 . We will identify a complex symmetric tensor f P S 2 C 3 with the C-linear operator, f : C n Ñ C n , as usual xf ξ, ηy " xf η, ξy " f pξ, ηq where η, ξ P R 3 .

1
Defining the Schwartz space of rapidly decreasing functions SpR n q on n-dimensional Euclidean (n " 2 or 3) space in the usual way, we extend this definition to (complex) vector fields SpR n ; C n q and complex symmetric tensor fields as SpR n ; S 2 C n q. The choice of Schwartz spaces and the use of complex vectors and tensors is convenient as we will rely heavily on the Fourier transform. Extension to Sobolev spaces for compactly supported fields follows using the usual apparatus as applied to scalar ray transforms [7]. We define the Fourier transform F : SpR 3 q Ñ SpR 3 q, by e´i xy,xy f pxq dx.
Given an orthonormal basis pe 1 , e 2 , e 3 q of R 3 , a tensor f P S 2 C 3 can be represented by the symmetric 3ˆ3 matrix pf jk q, f jk " f pe j , e k q. The Hermitian scalar product on S 2 C 3 can be written as xf, gy " ř 3 j,k"1 f jkḡjk independently of the choice of an orthonormal basis. Since only orthonormal bases will be used in this paper, we do not distinguish between co-and contravariant tensors.
We will need the partial Fourier transform F V : SpR 3 q Ñ SpR 3 q, for any kdimensional vector subspace V Ă R 3 . Given Cartesian coordinates px 1 , x 2 , x 3 q in R 3 such that V " tx|x k`1 " .... " x 3 " 0u. the partial Fourier Transform can be written as This result is independent of the choice of orthonormal coordinates and satisfies the commutation law F " The oriented lines in R 3 can be parameterized by points of the tangent bundle of the unit sphere S 2 in R 3 T S 2 " tpξ, xq P R 3ˆR3 | |ξ| " 1, xξ, xy " 0u Ă R 3ˆR3 .
An oriented line l Ă R 3 is uniquely represented as l " tx`tξ | t P Ru with pξ, xq P T S 2 . The Schwartz space SpT S 2 q is defined as in [4].
Our main object of interest is the Transverse Ray Transform (TRT) of symmetric rank two tensor fields which is defined by Jf pξ, xq " where P ξ : S 2 C 3 Ñ S 2 C 3 is the orthogonal projection onto the subspace tf P S 2 C 3 | f ξ " 0u. For example, for an orthonormal basis of the form pe 1 , e 2 , e 3 " ξq, the projection is expressed by 2 By contrast the scalar x-ray transform (or simply the ray transform) of a function f P SpR 3 q is defined by Xf pξ, xq " We also have the Longitudinal Ray Transform (LRT), defined on vector and tensor fields respectively by If pξ, xq " The x-ray transform and LRT can also be defined on a plane in R 3 . For η P S 2 , let η K " tξ P R 3 | xξ, ηy " 0u, Rη " tsη | s P Ru, η K C be the complexification of η K and S 1 η " tξ P η K | |ξ| " 1u be the unit sphere in ξ K . Given s P R, let sη`η K be the plane through sη parallel to η K and ι s,η : sη`η K Ă R 3 be the identical embedding. The family of oriented lines in the plane sη`η K is parameterized by points of the manifold T S 1 η " tpξ, xq | ξ P S 1 η , x P η K , xξ, xy " 0u such that a point pξ, xq P T S 1 η corresponds to the line tsη`x`tξ | t P Ru. We define the X-ray transform on a plane X η,s f pξ, xq and the LRT on the plane sη`η K X η,s : by the following formulae I η,s f pξ, xq " ş 8 8 xf psη`x`tξq, ξy dt, I η,s f pξ, xq " ş 8 8 xf psη`x`tξqξ, ξy dt (9) respectively. One can see that operators (3) and (9) are related. If f P SpR 3 ; S 2 C 3 q and ιη ,s f is the slice of f by the plane sη`η K , then I η,s pιη ,s f qpξ, xq " If pξ, sη`xq, for pξ, xq P T S 1 η . A typical experimental situation would involve rotation of the specimen (or equivalently the source and detector) about some finite collection of axes. In the scalar case of the x-ray transform rotation about one axis is sufficient as the problem reduces to the Radon transform in the plane, for a family of planes normal to the rotation axis. Consider first the case n " 2 in which X is identical to the Radon transform. In this case the formal adjoint B " X˚: SpT S 1 q Ñ C 8 pR 2 q, the back-projection operator, is well defined and given by 3 for φ P SpT S 1 q. We then have the inversion formula [7] for data φpξ, xq " X η,s f pξ, xq in the range of X F rf pxqspyq " |y|F rBφpxqs, which means that inversion is performed by a ramp filter (also known as a Riesz potential) applied to the back-projected data. This operation can be performed slice by slice to invert the x-ray transform for n " 3, in which case data is needed only for ξ P η K for some fixed rotation axis η P S 2 . In this case the slice-by-slice back-projection operator B η : SpRˆT S 1 η q Ñ C 8 pR 3 q. The reconstruction formula (11) becomes where s " xx, ηy.
Notice that the component xη, Jf pξ, xqηy " Xxη, f ηypx, ξq, for ξ P η K is simply the scalar x-ray transform in the plane through x normal to η. As observed in [8, Sec 5.1.6] this component can be reconstructed using any inversion formula for the planar Radon transform inversion plane by plane, including the one given in (12). Choosing six rotation axes η i so that the outer products ηη T are linearly independent in S 2 R 3 recovers f everywhere. As mentioned in the introduction this procedure is likely to be time-consuming as rotations are performed about six axes and yet for each ray only one measurement (out of a possible three) is used. The aim of this paper is to show, via a constructive inversion procedure, that f can be determined uniquely from the data Jf px, ξq for ξ P e K i , i " 1, ..., 3. Of course the diagonal elements f ii are already determined as above so our main task is to provide a reconstruction procedure for the off-diagonal elements.
For a given rotation axis η and direction ξ P η K we have in addition to the 'axial' component xη, Jf pξ, xqηy the 'non axial' components xζ, Jf pξ, xqηy where ζ " ξˆη and xζ, Jf pξ, xqζy with which to reconstruct the off-diagonal (in a basis including η) elements of f .

Relations between transforms
The aim of this section is to write the non-axial components of Jf pξ, xq, ξ P η K in terms of longitudinal ray transforms on transaxial planes.
We have P ξ pf q " Π ξ f Π ξ where Π ξ is the orthogonal projection matrix onto ξ K , and the 'off diagonal' component can be expressed as xpJf qpξ, xqη, ξˆηy " Since the vector field ηˆf η is orthogonal to η, its restriction to every plane sη`η K can be considered as a vector field on the plane, i.e., pηˆf ηq| sη`η K P Spsη`η K ; η K C q.
As we have seen in (2), that the TRT depends only on the projection normal to the direction of the ray. Now working in pη, ζ " ξˆη, ξq coordinates, let us consider, xJf pξ, xqζ, ζy, which can be transformed into xJf pξ, xqζ, ζy " Let us parameterize ξ in the usual sense as " ξAdj e K 3 pf qξ, where Adj e K 3 pf q denotes the adjugate matrix of the slice of f restricted to the plane; of course this is nothing other than the conjugation of P e3 f with a right angle rotation about the e 3 axis. Hence using the above, we see that for ξ P η K xJf pξ, xqζ, ζy " which is the LRT transform of the two dimensional adjugate of the slice of f restricted to the plane. We notice that this is also exactly the transverse ray transform in the planar case. The results of this section can be summarized in the following Lemma: Lemma 1 Let f P SpR 3 ; S 2 C 3 q be a symmetric tensor field. The equations hold for every s P R and η P S 2 , where pJ 1 η,s f q " xpJf qpξ, sη`xqη, ξˆηy, pJ 2 η,s f q " xpJf qpξ, sη`xqζ, ζy.
In the next section, we will transform (16) and (17) to algebraic equations by applying the Fourier transform to back-projected data.

Main algebraic equations 4.1 Curl Components of Tensor and Vector Fields
We can now transform (16) and (17) to algebraic equations by applying the Fourier transform. We only require what [8] refers to as tangential component τ g P C 8 pR 2 q of a vector field g P C 8 pR 2 ; C 2 q, which is defined by Here the vector y K is the result of rotating y by π{2 in the positive direction. Of course, one can understand (18) as the 2D curl of a vector field in Fourier (frequency) space. The manifold T S 1 can be identified with RˆS 1 by the diffeomorphism pp, ξq Þ Ñ pξ, pξ K q for pp, ξq P RˆS 1 . Therefore the derivative B Bp : SpT S 1 q Ñ SpT S 1 q is well defined. For a vector field f P SpR 2 ; C 2 q, the tangential component of the Fourier Transform F rf s is recovered by the LRT, If , by the formula We see in [4] and [9], the tangential component, τ g P C 8 pR 2 q, of a tensor field g P C 8 pR 2 ; S 2 C 2 q is defined by pτ gqpyq " |y| 2 tr g´xgpyqy, yy. (20) This is exactly the Fourier transform of the single unique non-zero component of the compatibility tensor of Barré de Saint-Venant in the plane which is also sometimes described as the curl of a symmetric tensor field. For f P SpR 2 ; S 2 C 2 q, the tangential component of the Fourier transform F rf s is recovered from the LRT, If, as For φ P SpT S 1 q, the function Bφpxq is C 8 -smooth and bounded on R 2 but does not decay fast enough to be in the Schwartz class. Thus we understand the Fourier transform in the distribution sense in (19) and (22).

Derivation of System of Equations
Let f P SpR 3 ; S 2 C 3 q be a symmetric tensor field and denote by f 1 " F η K rf s P SpR 3 ; S 2 C 3 q the partial Fourier transform of f . For any s P R, the restriction of the vector field ηˆf 1 η to the plane sη`η K coincides with the 2D Fourier transform of pηˆf ηq| sη`η K . This is We then apply formula (19) to the vector field pηˆf ηq| sη`η K , giving Note that (18) gives τ ppηˆf 1 ηq| sη`η K qpsη`yq " xηˆf 1 psη`yqη, ηˆyy " xf 1 psη`yqη, yy.
Upon substitution of (25) into the LHS of (24) and applying the one-dimensional Fourier transform F Rη taking s to σ gives xf pση`y 1 qη, y 1 y " wheref is the three-dimensional Fourier transform F rf s. Since y 1 P η K and σ P R, we let y " ση`y 1 , where y 1 " Π η y. Hence the previous formula (26) can be written as Note that (27) is identical to the off-diagonals for the TTRT operator case in [4]. Moreover this will just reconstruct the solenoidal part of the off-diagonals since the Fourier transform interweaves with the solenoidal part. For any s P R, the slice ιη ,s f 1 coincides with the 2D Fourier transform of the slice ιη ,s f , i.e., ιη ,s f 1 " F η K rιη ,s f s, where the Fourier transform on the plane sη`η K . Henceforth, we refer to the adjugate of the slice of f 1 restricted to the plane as Adj η K pιη ,s f 1 q " h 1 . Upon application of formula (22) to h 1 , we see Using Lemma 1 we can rewrite the above as rτ ph 1 qspsη`yq " Now, we apply formula (20) to the field g " h 1 P Spsη`η K ; S 2 η K C q to give rτ ph 1 qspsη`yq " |y| 2 tr h 1 psη`yq´xh 1 psη`yqy, yy for y P η K . (30) Next, substitute (30) into (29) to give |y| 2 tr h 1 psη`yq´xh 1 psη`yqy, yy " By applying the one-dimensional Fourier transform on Rη to the above, we obtain |y 1 | 2 trĥpση`y 1 q´xĥpση`y 1 qy 1 , y 1 y " for y P η K . As before, employing a change of variables, y " ση`y 1 , transforms the above to In the following statement, the results are summarized.
Lemma 2 Let p f be the 3D Fourier transform of a symmetric tensor field f P SpR 3 ; S 2 C 3 q. For a unit vector η P S 2 , the following equations hold with the additional condition that p h P SpR 2 ; S 2 C 2 q, is defined to be the 2D adjugate of f restricted to the plane x p f pyqη, Π η yy " λ η pyq and (34) hold on R 3 , with the right hand sides defined by The partial derivative B Bp : SpRˆT S 1 η q Ñ SpRˆT S 1 η q is defined with the help of the diffeomorphism R 2ˆS1 η Ñ RˆRˆT S 1 η , ps, p, ξq Þ Ñ ps, ξ, pξˆηq. Given the data Jf | η K , right-hand sides λ η pyq and µ η pyq of equations (34) -(35) can be effectively recovered by formulas (36) -(37).

8
The system of equations (38) can be solved to give The result is summarized in the following theorem Theorem 1 A symmetric tensor field f P SpR 3 ; S 2 C 3 q is uniquely determined by the data Jf pξ, xq for ξ P η K i , i " 1, 2, 3 where pη 1 , η 2 , η 3 q forms an orthogonal basis.

Alternative formulae for diagonal components
While the diagonal components f ii are easily determined as we have seen, there is an alternative more complicated procedure to recover them. As this uses different data it can also be viewed as a compatibility condition on the three axis data. Consider (35) and (37). When η " e 1 , we have trĥ "f 22`f33 , and h "¨0 In the same manner as above (λ i ), we achieve a system of equations for µ i Let us relabel the RHS of the above as In this way the solution of (43) can be written aŝ Upon substitution of the off-diagonals and µ's into (45), we obtain 6 Insufficiency for two axes To reduce the data acquisition time, experimentalists would want to rotate the specimen on its axis as few times as possible. We show that in the general case two orthogonal axes are insufficient by considering components in the null space of the TRT. Thus J η f " 0. If we had two orthogonal axes, say η " e 1 , e 2 , then xη, J η f ηy " 0. This implies that f 11 " f 22 " 0. From the definition of λ η pyq and µ η pyq, λ 1 " λ 2 " 0 and µ 1 " µ 2 " 0. The system of equations for the off-diagonals (38) gives us Moreover the system of equations corresponding to the other non-axial components, (44), gives Due to the argument at the start we can rearrange (50) as From the above, sayf 33 is arbitrary and consequentlyf 13 andf 23 are determined aŝ f 23 "´y 3 2y 2f 33 andf 13 "´y 3 2y 1f 33 .
Using the values obtained in (53) and substituting into (48) we can writef 12 aŝ Thus all the off-diagonal components in the tensor field are determined throughf 33 which is arbitrary. Hence two axes are insufficient.

10
In the potential case f ij " Bu i {Bx j`B u j {Bx i where u P SpR 3 ; C 3 q. This is important for applications in that a linear strain tensor f has this form where u is twice the displacement field. Without loss of generality suppose that data is known only for rotations about η " e 1 , e 2 . We have immediately f 11 , f 22 and hence by direct integration twice u 1 and u 2 . We now also have f 12 from the partial derivatives of u 1 and u 2 . It remains only to find u 3 . Multiplyingf 12 by y 1 y 2 andf 13 by y 1 y 3 and adding both of them in (41) gives This gives us f 13 in terms of known data and as Bu 1 {Bx 3 is known we have Bu 3 {Bx 1 and hence u 3 . We summarise in the theorem Theorem 2 A potential f P SpR 3 ; S 2 pC 3 qq is determined uniquely from Jf pξ, xq restricted to ξ P η K 1 and ξ P η K 2 where η 1 and η 2 are orthogonal. This result is of considerable practical importance as it means that stain tensors, in a scheme such as that envisaged in [5], can be recovered from rotations about only two axes.
We now show that in general a potential f cannot be recovered uniquely from a one-axis rotation by constructing a general element of the null space. Suppose we rotate only about e 1 we have immediately f 11 " 0 and asf ij " y iûj`yjûi we see u 1 " 0. Now as λ 1 and u 1 are zero and as µ 1 " 0 y 2û2`y3û3 " µ 1 2py 2 2`y 2 giving no new information. So u must satisfy u 1 " 0 and Bu 2 {Bx 2`B u 3 {Bx 3 " 0. For example if u 2 is arbitrarily specified, then 8 Numerical results

Forward model
In order to generate data sets, we need to implement a discretized version of the operator J described in (1) as a matrix which will calculate integrals of projections and act upon generated strain fields represented by an array. Instead of calculating the whole matrix, we generate it one row at a time (on the fly) which corresponds to one individual source-detector pair for one of the three components in P ξ f described by (2).

Discrete representation of the tensor field
The discretized tensor field is stored as a NˆNˆNˆ6 vector, containing the 6 distinct values of the symmetric second rank tensor field for each voxel in a NˆNˆN voxel grid. We increment first by the tensor component number, then the position x 1 , x 2 and finally x 3 . Furthermore, the data is represented by a 3ˆn θˆhˆwˆ3 multidimensional array (5 dimensional), where we use 3 rotation axes pη " e 1 , e 2 and e 3 q, n θ angles steps for tomographic acquisition around each axis and hˆw represents how many rays in the horizontal and vertical direction. The factor of 3 is the number of independent values of Jf in (2) which we integrate along each ray.

Methodology
We simulate an experimental setup with parallel rays passing through a specimen. Sources and detectors consist of arrays in an equally spaced grid, either side of the object being scanned. The ratio of the number of rays in the horizontal to vertical direction is 4 : 3. The source-detector pair is kept fixed and the object is rotated. We follow the procedure of [10, Chapter 5.1.4] to calculate the approximate integral along a line through a voxel grid. This will give us the contribution of each voxel to the total integral for a given ray. For a given tensor we need to calculate the projection on to the plane perpendicular to the ray. To simplify this, we rotate the coordinate system. We extract the two appropriate components (relating to the axial and non-axial components). Since we have the contribution of each voxel to the integral, the length of intersection of the ray with the voxel, from ray tracing we can form the sum of these intersection lengths with the voxel values to form the approximate integral. For our numerical experiments, phantoms were generated inside a cubic grid measuring 405ˆ405ˆ405 voxels and measurements were simulated for a source/detector grid with 405ˆ540 pixels. The number of rays in the vertical direction was taken to be the image height (i.e. 405 pixels). The pixel grid was then down-sampled by a factor of 3 to 135ˆ180 pixels by the process of binning. Finally 1% Gaussian pseudo-random noise was added before reconstructing on a courser voxel grid measuring 135ˆ135ˆ135. The number of angles (projections) was 240 per rotation axis.

Generating Phantoms
For input into the forward projector we generate two different phantoms or test fields. The first one only has smooth features and is expected to be less sensitive to algorithmic instabilities. The second phantom contains sharp edges and is designed to highlight the limitations of the explicit reconstruction algorithm for discontinuous strain fields.

Phantom 1: smooth
Phantom 1 is constructed from smooth Gaussians which makes it relatively easy to reconstruct. We define a cubic domain r´1, 1s 3 on which the components of f are supported, defined by 3-dimensional Gaussians b α pxq for each of the components f ij according to Table 1, where b α pxq " α expp´50|x´a| 2 q.

Phantom 2: sharp
Phantom 2 is constructed to contain sharp edges, to highlight non-linear strain fields which is not quite compatible with the reconstruction algorithm. As in the smooth case, we define f on the cube r´1, 1s 3 , but set f ij to be the characteristic function of I 1ˆI2ˆI3 , according to Table 2.

Reconstruction Procedure
The recovery of axial components are relatively straight forward as this is just planeby-plane Radon inversion. Hence, we apply a ramp-filter (Ram-Lak) to the data and back-project to achieve the diagonal entries for each rotation axis. Since, we can think of back-projection as the dual operator of ray integration, we reuse the ray tracing code to implement back-projection as the transpose of ray integration.
From (36), we see that the simulated data values J 1 that are collected for each plane need to be differentiated in the p-direction, before any back-projection takes place. To implement this, we perform a regularised derivative, which we carry out 13 in the Fourier domain using a Hamming window regularisation. After performing a one-dimensional Discrete Fourier Transform using the Fast Fourier Transform FFT algorithm, we multiply by´ikwpnq, where wpnq is the Hamming window. For a discrete signal of length N, labelled by n " 0, ...., N´1, the Hamming window wpnq is given by The result is returned to the spatial domain using the inverse FFT algorithm. Following the filter we back-project the differentiated plane by plane data onto the voxel grid and the tangential vector field components (i.e. λ) are calculated using a three dimensional FFT algorithm and application of a ramp-filter in frequency space. Then equations (41) are used to recover the off-diagonal terms in frequency space. The only exception is the voxel py 1 , y 2 , y 3 q " p0, 0, 0q, wheref 12 ,f 13 andf 23 are undefined. Here, the value is set using linear interpolation from nearby voxels. To complete our reconstruction, we employ the three dimensional inverse FFT to recover f ij .

Results and summary
Below we illustrate the results of our implemented reconstructions which clearly show the performance of the algorithm on two different tensor fields. As expected, the smooth (Gaussian) field is reconstructed well. However, when we introduce sharp edges in the components of a field such as a crack, we see that as expected the reconstruction is inaccurate and artefacts are visible. The artefacts increase as more noise is added.    Overall the procedure we describe would provide an efficient method for x-ray diffraction tomography. It handles smooth tensor fields well but for discontinuous strain fields it would be sensible to try different modified ramp filters or explicit regularization methods such as total variation. If would be possible in the case of a tensor field that is the result of an infinitesimal strain to use only two axes of rotation, and recover the displacement field directly. However it might be better in practice to use the general procedure and then verify to what extent the compatibility condition holds on the reconstructed tensor.
We have pointed out that as there are two distinct methods of calculating the diagonal components this provides a consistency condition on the data. As the plane-byplane data is written in terms of scalar, vector and tensor longitudinal ray transforms, and the ranges of these operators can be determined in the plane case as a singular function expansion in a suitable Hilbert space. In fan beam coordinates the singular value decomposition of the ray transform is given by [3]. See also [1] and [2] for a parallel beam formulation.
The Helgason-Ludwig range conditions for the scalar Radon transform in the plane [7,Sec II.4] simply states that the kth moment of the data 8 ş 8 s k Xf pξ, sξ K q ds, when it exists, is a polynomial of degree ď k in ξ. For the LRT of a rank m tensor the condition is simply degree ď k`m. A deeper connection between these conditions and the singular function expansion is given by [6]. Such consistency conditions, characterizing the range of the forward operator are of great assistance in diagnosing errors and unaccounted for physical effects in experimental data.
Another avenue worth considering on the practical side is to develop a reconstruction algorithm involving general (non-orthonormal) axes. In experiments it is often not feasible to rotate the specimen through 90˝and remain in the field of view of the measurement system.
Explicit reconstruction algorithms such as the one we have given are useful practically for data that is complete and uniformly sampled. For partial, sparse or irregularly sampled data representing the forward problem simply as a sparse matrix and solving using iterative algorithms with explicit regularization is generally better, although typically requiring large amounts of memory and parallel processors.