Histogram tomography

,


Introduction
Many tomographic imaging problems fall in to the broad category of rich tomography problems meaning that for each line our data are more than a single scalar.For example in many absorption tomography problems we might record a spectrum, a real valued function on an interval corresponding to the amount of radiation at each frequency or wavelength transmitted along that line.In other examples we might apply a narrow beam of radiation, x-rays, neutrons or electrons, and measure a diffraction pattern, a real valued function of two variables.From this data we expect to be able to recover more than a single scalar image.We might expect a vector quantity such as velocity or magnetic field, a tensor such as strain, or the concentration of more than one chemical.Such techniques have been described as Rich tomography.Alternatively we might expect to recover a single scalar image with fewer projections than would be necessary for traditional scalar measurements.In many cases acquiring a projection from raster scanning a beam is time consuming and reconstruction with fewer projection desirable, in other cases where rotation about an axis is limited it is essential.
In particular numerous problems where the data are spectral allow us to infer the distribution of a scalar variable along each ray from the measurement of the spectrum for that ray, rather than simply its integral.In a practical setting this will often mean a histogram, in which the relative frequency of each value along the line that falls in to each bin of the histogram is recorded.In the limiting case as the size of the bins vanish this is the distribution, explained in the next section.
As a general rule standard tomographic problems involve a transport equation along the ray.In this case the values recorded are cumulative.As a ray travels from one voxel to the next the transport equation transforms the value input to one voxel to that output to the next.In the simplest case of attenuation tomography this results in Beer-Lambert law and the logarithm of the measurement at at the output is a line integral [19].In other non-Abelian tomography problems the transport equation cannot be solved using an exponential and we have a non-linear integral operator along rays.Polarimetric neutron spin tomography provides an example [21].The problems we treat in this paper of a quite a different nature.As an archetypal example consider infrared absorption tomography [3,18] where a certain chemical species absorbs a narrow band of infrared light, but the centre frequency of that band is a known monotonic function of temperature.Assuming the species is uniform distributed over the medium a portion of the spectrum around the shifted bands, when transformed to the temperature variable, directly gives us the distribution of temperatures.We know which temperatures occur along the ray, and how often they occur on the ray, but from one projection we do not know the location where these temperatures occur.Interestingly we do know the maximum and minimum temperatures immediately.
In the next section we review the definition of the distribution of a function with respect to a measure, and the concept of moments of the distribution.This has already been explained very clearly in the case of the Doppler tomography of vector fields by Andersson [4], who used the idea of moments of the distribution.We go on to consider scalar histogram tomography problems and what can be deduced from limited data.We then review concepts from Sharafutdinov [24] of symmetric tensor fields and their ray transforms.We return to the Doppler moment transform and show that from the second moment we can reconstruct a potential field directly.We then review the method of using the neutron transmission spectra near a Bragg edge and show that in principle the histogram of the strain in the ray direction can be recovered from this data.It was mentioned in [17] that as the linear elastic strain is potential the longitudinal ray transform derived from the Bragg edge yields no data on the interior strain.We show that the moment data of the distribution gives rise to non-linear partial differential equations in the strain and give explicitly the second moment in the two-dimensional case.We go on to discuss x-ray strain tomography and the potential for extracting partial histogram data in this case, giving an explicit example of two strain distributions along a line resulting in the same diffraction pattern.
In this paper we have not specified the smoothness assumptions on functions and kept functional analysis to a minimum.This is an effort to make the paper accessible to the diverse applications communities that work with histogram tomography data.The details are easily filled in by anyone familiar with the necessary analysis for the Radon transform [19] or the ray transforms of vector and tensor fields [24].

Distribution and moments
The idea of the probability distribution of a random variable is a familiar, but in this paper we will need the concept of the distribution of a function with respect to the standard measure on Euclidean space rather than a probability measure.The (Lebesgue) measure λ on Euclidean space R n is a function that assigns a measure of size to a set [7], for example total area (n = 2), volume (n = 3).For n = 1 the measure of an interval λ[a, b] = b − a is simply the length.
Given a real valued (measurable) function f on a bounded measurable subset Ω ∈ R n we define the push forward measure f * λ on R where for any interval This definition is the same as the probability distribution of a random variable f except that we do not assume a probability measure.We also define the cumulative distribution Φ f (y) = f * λ(−∞, y], the measure of the sub-levelset of f up to y.Even though Φ need not be continuous we have Φ f = φ f in the sense of distributions.See figure (1) for an example of the distribution and cumulative distribution of a function.
In the case where f (x) = c a constant then clearly φ f (y) = δ(y − c)λ(Ω).By contrast for a differentiable function f and (for simplicity dimension n = 1) let y ∈ R be such that f −1 (y) is finite, and f (x) = 0 for all x ∈ f −1 (y) (not a critical value), then As illustrated in Figure 1 isolated critical points of a smooth function result in singularities of the distribution.This means that one can in principle count distinct critical points and know their critical values is the integral of the function.If f were a random variable and the measure normalized to be a probability measure, this would be the expected value.We also define the k-th moment for k ∈ N as While the above is perhaps familiar for non-negative functions [4,Cor. A.1.]proves it for functions that can take negative values.We will need to consider in some cases generalised functions, such as Dirac δ-functions, which are also called distributions.To avoid confusion with the distribution of a function, in this paper we use the term generalized functions.
Let L θ,p = {x : x 1 cos θ + x 2 sin θ = p} = {x : x • Θ = p}, where Θ = (cos θ, sin θ), be a line in the plane.We define the Radon transform on functions on the plane by where Θ ⊥ = (− sin θ, cos θ).We generalise this to the histogram Radon transform the distribution of f restricted to each line.Clearly we can recover the Radon transform by taking the first moment with respect to y Rf (θ, p) = Hf (θ, p, •).
We note also that This shows that for the scalar case each of the moments produces no more information than the Radon transform and for a non-negative function exactly the same data.Indeed a non-negative bounded function is determined completely by its moments [2] so the data are identical.It is interesting to note however that while fitting a function f to its histogram tomography data Hf is a non-linear problem, each of the problems (6) is linear.One of the main questions, for the scalar case in the plane, is to identify interesting subsets of lines for which Hf determines f , without the reduction in stability associated with the limited data problem for the Radon transform.
We will need the Radon plane transform of functions of 3-space later.This is simply the integral over a plane normal to a unit vector Θ ∈ R 3 a distance p from the origin so the notation is the same as for the integral over lines in the plane 3 Recovery of level sets and discrete tomography A key insight that provides the connection between between the histogram tomography transform and the body of classical work on tomography is that the cumulative histogram gives the Radon transform of sub-level sets.To unpack that statement,the sub-level set of f at y is S f (y) = {x : f (x) ≤ y},and is just the total length (measure) of the intersection of L with the set which is also the Radon transform of the characteristic function χ S f (y) .This leads us to the chord length problem, which is exactly to recover a bounded subset S in the plane from this data for some set of lines.For convex sets S Hammer's x-ray problem is to recover S from Rχ S (θ, s) for a specified set of θ and all s.The solution to this problem is unique for four values of θ that are not rationally related [10].
For non-convex sets less is known.But there is another approach.Local tomography methods use limited data for the Radon transform to recover the singular support of a function.In the case of χ S f (y) the singular support is exactly the contour f −1 (y).As an example [9] shows that it is only necessary to know Radon data for all the lines passing through a neighbourhood of the singular support to recover it uniquely.This means for each contour one only needs the subset of the histogram tomography data for lines passing close to that contour y value.
As characteristic functions take only the values 0 or 1, Radon transform inversion for these functions falls into the approach of what is called discrete tomography [15], that is tomography where the range of the unknown function is a finite set.In the {0, 1} case this is also called geometric tomography [11] with the emphasis mainly on convex bodies.As far as numerical algorithms for discrete tomography the Discrete Algebraic Reconstruction Technique (DART) is a popular variation of the standard iterative algorithm used in tomography adapted to functions with a binary range [5].
As mentioned in Sec 2 the distribution of a function along a line has singularities corresponding to critical values, so if there are isolated critical points with distinct critical values these can be seen on the distribution.Suppose now we have a function f on the plane with isolated critical points.Let x c be a critical point with f (x c ) = y c and ∇f (x c ) = 0 then for any line through x c , that is L θ,p such that Θ • x c = p the distribution Hf (θ, p, •) will have a singularity at y c .This could be used to gain information about critical points from a limited data from the histogram Radon transform.
As well as discretization of the range of values of f it is common to discretize the domain in to pixels or voxels.If the pixels are reduced to points and the rays still taken to have no thickness this results in tomography problems on finite subsets of integer lattices, also known as discrete tomography but in this case with more of a number theory emphasis.In a more general setting the function is replaced by a function on the vertices of a graph and a set of paths through edges take the role of rays.When the values of the function can take is a discrete set, in this context refereed to as colours, and the rays are simply subsets of vertices.The number of each colours occurring in each ray is specified we get a discrete tomography graph colouring problem.Mainly the computational complexity of such problems is studied in [6], rather than specific algorithms.It is however notable that in a fully discrete setting this is the histogram tomography problem we describe with labels rather than numerical values.Such problems arise in scheduling of work rosters and have not drawn much attention from those working on spectral tomography methods.

Tensor ray transforms
Given a vector field f the longitudinal ray transform (LRT) is defined as while for a rank-2 symmetric tensor field f the LRT is The definition is extended to rank k symmetric tensor fields in the same way with the contraction, with respect to ξ, k times.The LRT has a null space consisting of potential tensor fields.In the case of vector fields this is just the usual definition, f = ∇u for a scalar u.Potential rank-2 tensor fields are those that can be expressed as f = (∇u + ∇u T )/2 for some vector field u.In general, following [24], we define the operator d from rank-k to rank-(k + 1) tensor fields formed by differentiation and symmetrization.
A rank-k symmetric tensor field f is said to be potential if f = du for some rank-(k − 1) field u.The Saint-Venant tensor W (f ) has rank-2k with the property that W f = 0 if, and on a simply connected domain only if, f is potential.For k = 1, (W f ) ij = f i,j − f j,i is the tensor of skew symmetric derivatives of f .Here subscripts after a comma denote differentiation.In dimension 3 the nonzero components can be rearranged as ijk f j,k , which coincides with the components of the curl, where is the permutation symbol and sums are taken over repeated indices.For k = 2, W f is the familiar Saint-Venant compatibility tensor that is well known in elasticity.By a similar rearrangement of components, for the case n = 3, k = 2 the non-zero components of the rank-4 Saint-Venant tensor can be rearranged as a rank-2 Koröner tensor Kf .For the general case (which we should call the Georgievskii-Koröner tensor) see [12].In general a rank-k symmetric tensor can be expressed as the sum f = f sol + f pot of a solenoidal part f sol with W f = W f sol and a potential part f pot = du for some u [24].With imposition of suitable far field or boundary conditions the decomposition becomes uniquely determined.
For n ≥ 2 there is an explicit reconstruction formula for W f (or Kf ) from If of filtered back projection type [24].This determines the solenoidal part of f , so the solution is unique up to the null space of potential fields.
Let P ξ be the projection of a symmetric second rank tensor field on to the plane perpendicular to ξ, then the transverse ray transform (TRT) is defined as We consider the important case of dimension n = 3.For a direction η normal to ξ η so in any plane normal to η this is simply the Radon transform of the component η • f • η.This means there is a simple reconstruction for six suitably chosen [17] directions η.
Both these problems have histogram tomography version for the HLRT the data is the distribution of f (x + tξ) • ξ for vectors and ξ • f (x + tξ) • ξ for rank-2 tensor.By contrast P ξ f (x + tξ), along the ray x + tξ, is a tensor field so the HTRT would consist of the joint distribution of the values along the ray of the components of the tensor normal to the ray direction.
In the case of the HTRT, one special case is that we have the distribution of η • f • η along lines on a plane normal to η.As this is a Radon transform we have reduced to the scalar histogram tomography problem for the normal component on this plane.As the TRT data is consists of three independent components for each ray (equivalent to a symmetric 2 × 2 matrix).The distribution of these values can be considered as a joint distribution.But for the argument just outlined only a marginal distribution, the distribution of the normal component to the plane, is needed.

Doppler velocimetry
As Sparr [26] and Schuster [23] explain Doppler velocity tomography data are already understood as the distribution of velocity components along a line, in the direction of a line.The technique is used in laboratory and field studies to determine a velocity field from the frequency shift of a narrow acoustic beam.In our terms, after suitable signal processing, this is the HLRT of the velocity field.It is also called the Doppler spectral transform.The first moment is typically used and of course this gives only the solenoidal part of the velocity leaving the potential part to be determined by other means.Andersson [4] shows that the higher moments determine a system of non-linear partial differential equations for the potential part.In this section we give a simple procedure to recover potential velocity field from the HLRT.
Let f = f sol + f pot be the decomposition of the velocity field in to potential and solenoidal parts.As we can recover the solenoidal part of the vector field from the first moment we can assume that f sol is known, hence we know the distribution φ L,sol of f sol • ξ along a line L in direction ξ.Unfortunately this does not lead to a method to recover the distribution of the potential part as we do not know where along the lines the values arise.In probability theory there is a large literature on what can be determined about the distribution of the sum of two dependent random variables where the marginal distribution is known but not the joint distribution.In general if the marginal distribution is known for a multivariate random variable the joint distribution is defined only up to a Copular (the joint cumulative distribution after a change of variables so marginal distributions are uniform) by Sklar's theorem [25], and the set of possible joint distributions with given marginals is called a Fréchet class.The Fréchet-Hoeffding Theorem gives bounds on the joint distribution [20].While this suggests an inyteresting avenue of future work we will now confine ourselves to the special case of reconstructing a field with zero solenoidal part.
In this case we have f = du for a scalar u.Notice the second moment is nothing but the LRT of du du.From [24] we know we can recover the Koröner tensor [12] K mn = mik nj (u ,i u ,j ) k = mik nj u ,ik u ,j , noting that the third derivative drop out as each have a pair of indices over which a skew symmetric sum is taken.Recall that for a square matrix A, AAdj A = det AI where Adj is the adjugate matrix, the transpose of the matrix of cofactors, and I the identity matrix.For the 3 × 3 case and a simple calculation shows Adj Adj A = (det A)A, and (det A) 2 = det AdjA.So when det A = 0 (8) shows that we can deduce A from Adj A up to sign as . When the determinant is zero in the case rank A = 2 the columns of Adj A span the null space of A. In the case rank A = 1 or 0, Adj A = 0. We see where (d 2 u) ij = u ,ij second derivative matrix.When det K = 0 we can recover d 2 u up to sign, but when the determinant is zero K does not uniquely determine d 2 u.Explicitly we have that is known except for points at which det K = 0. We needs some additional data to determine boundary conditions and the sign of u.Indeed it was pointed out in [26] that for radial rotation invariant vector field f the HLRT of f and −f are identical.Note that trace d 2 u = ∇ 2 u, so we consider the Poisson equation and given boundary data for u (for example outside an object it is zero) we can solve uniquely for u given K, when det K = 0, up to a sign ambiguity.We now have a constructive procedure for recovering u from the second moment of the histogram LRT of the potential vector field du up to sign.Apply the reconstruction formula of [24] to this data to recover K then solve Poisson's equation.
The extension to general velocity fields, with non-zero solenoidal part, is not as simple and we know little more than was already presented by Andersson [4].The Solenoidal part can be assumed known from the first moment data.The Kröner tensor of the second moment is a second order system of partial differential equations for the potential with the solenoidal part as coefficients.We hope to investigate uniqueness of solution in subsequent work.

Neutron Bragg edge strain tomography
Consider a sample of solid material that consists of identical crystals with a uniform random distribution of directions.Each crystal has crystallographic planes denoted by Miller indices, which are triples of integers.A narrow beam of neutrons can be produced with a known distribution of kinetic energies.From a neutron scattering point of view the neutrons behave as waves with a wave length inversely proportional to the energy, and this wavelength can be made to be close to the separation of the crystallographic planes.An incident wave is scattered by parallel crystalographic planes with separation d if the wavelength is a multiple of 2d sin θ.If one measures the neutrons that are transmitted, the scattered ones are lost.Looking at a transmission spectrum, such as Figure 2, one sees a dramatic jump in the count of transmitted neutrons at a wavelength corresponding to θ = π/2 for some crystallographic plane; this is called the Bragg edge.When the material is subjected to a linear elastic strain ε the separation between planes with normal vector η undergoes a relative change proportional to the strain η • εη component.See [22] for details and further references.If the strain was uniform along the path of the neutrons this would simply move the Bragg edge.
Under ideal conditions a Bragg edge can be modelled by linear function of wavelength multiplied by a Heaviside step function.Now consider a distribution of strains along the neutron path with a given histogram, the superposition of these Bragg edges is simply the sum of Heaviside functions, and hence the cumulative histogram, again multiplied by a linear function.On the enlarged plot of the Bragg edge corresponding to the Miller index 110 in Figure 2 we see that this has the appearance of a sigmoid curve as expected of a cumulative distribution.To be precise we have to remove the linear component, which can be done by differentiation and then adjusting an added constant.See figure 3 for a cartoon.We conclude that from this part of the spectrum we recover the HTLR of the strain tensor.
Previously [1] investigated the possibility of reconstructing strain from an average of the shift of a Bragg edge for each ray.L. and Withers pointed out [17] that this amounts to the LRT of tensor field ε = du that is potential, except for discontinuities at the entry and exit of the ray from the material.In this case u is the displacement vector field.In practice this means that the overall change in dimensions of the material along the ray can be recovered from this data and not strain components in the interior of the material.This does at least provide boundary data to solve the Lamé equations of linear elasticity, a second order elliptic system for u, by finite element or difference methods.Provided the Lamé coefficients are known.This approach was Figure 2: Neutron spectra from [22] showing the 110 edge on a magnified scale.The wavelength is given in Angstrom = 10 −10 m.Note that the 211 edge is closer to the ideal behaviour of a step in an otherwise linear trend Figure 3: Shifted Bragg edges produce a saw tooth effect but the derivative is a histogram up to an added constant taken for axisymmetric objects by Gregg et al [13], in the general case for synthetic data by Wensrich et al [27], and experimentally for a general case by Henriks et al [14].
The strain is assumed to be potential and we have the distribution we known its the k-th moment and this is the LRT of the symmetric powers du • • • du= du k .These are in general not potential so the moments will be non-zero.So we can recover the Saint-Venant tensor of these powers.As outlined in the next section this results in a nonlinear system of partial differential equations for each k.Unfortunately we know of no trick like that employed for vector fields that reduces these to a linear system.Essentially we are in the same situation described by [4] for the vector case before this paper.While solving a nonlinear system is undesirable, and we expect more computationally demanding than solving the linear Lamé system, it does offer a new possibility at reasonable cost.The results of [14] could be verified by computing these Saint-Venant tensors for several values of k and comparing with the symmetric powers of the reconstructed strains.As there are few options for measuring the strain independently this may offer a valuable independent check on the results and hence the validity of the assumptions.The details of the calculations are in the next section.

Second moment in the rank 2 case
In general the Georgievskii-Koröner tensor for a rank m tensor f in dimension n has rank m(n − 2).In the case m = 4, n = 3 it is given by where we have followed Georgievskii [12] in labeling indices with a pair of subscripts.Now consider the case f = du du where u is a vector field with round brackets indicating symmetrization over the indices.In this case K is a fully symmetric rank 4 tensor.In general the terms result from differentiating a product of derivatives of the u i four times.
From the product rule we would expect terms including fifth derivatives of the components of u but the skew symmetrizing over pairs of indices that appear as derivatives removes these terms.In practice the components of K consist of sums of pairs of third partial derivatives, and products of 4th and 2nd partial derivatives of u.For dimension 2, or restriction to a plane, K is a scalar consisting of a fourth order non-linear PDE for u.
Unless we have additional information about u this is insufficient to determine u uniquely.However in the three dimensional case data for three or more families of planes normal to three independent vectors has the possibility of yielding a unique solution.Unlike the case of scalar u we do not know a neat transformation that eliminates the non-linearity.

X-ray diffraction strain tomography
In the neutron Bragg edge case we ignored the scattered radiation, but it is possible instead to measure scattered particles/waves.In the Bragg edge case it was possible to use a parallel beam of neutrons, as only the transmitted radiation was recorded.But to capture a diffraction pattern for the material along a ray the simplest approach is to use a single narrow beam, or possibly a small number of beams far enough apart that the diffraction patterns have large regions in which they do not overlap.In [17] we considered the case of a monochromatic x-ray beam and a polycrystalline material.We showed that by taking certain radial moments of the diffraction pattern we could recover the TRT of the strain.In this section we consider what histogram data is available form the diffraction pattern in this case.Intuitively the diffraction pattern is a function of two variables and yet the transverse strain has three independent components so we would not expect to recover the full HTRT of the strain.Consider one specific ray in the direction ξ and the projection of the strain in the plane normal to ξ a 2×2 matrix valued function A(s), as in [17].The diffraction pattern for the voxel at s is an ellipse determined by the matrix A(s), and the total diffraction pattern the sum of these where y is the coordinate in the plane of the screen where the diffracton pattern is formsed and φ is the density function, which is a generalised function supported on the curve A(s).Consider a point y = r(cos θ, sin θ) on the screen.The range of integration in ( 12) is then the points (a 11 , a Note however that the family of planes is given by the two parameters r, θ where as the full set of planes in three dimensional space is three dimensional, so we do not expect to be able to recover a general φ from this.Of course the distribution φ is not completely general but is a generalized function supported only on a curve, so one might think the situation is somewhat better.Using the Fourier slice theorem for the ray transform we see that the data g determines φ on a cone.Let α = (α 11 , √ 2α 12 , α 22 ) be the Fourier transform variables then the components determined by g(y) are Rφ(α/|α|, |α|) = φ(α) with α on the cone α 11 α 22 = 4α 2  12 .We see that if the difference between two distributions has a Fourier transform vanishing on this cone then they will produce the same diffraction pattern g.
To make an explicit example we consider biaxial strain where the principle axes of the strain rotates along the ray though half a turn.The diffraction pattern g will be supported and positive on an annulus with radii corresponding to the principle strains, and rotationally symmetric.Now consider a uniaxial strain in the transverse plane, in this case the bulk strain ranges between the principle strains of the previous case.This will also be positive and supported on the same annulus.To make the diffraction patterns the same we have to calculate the value of g in the first case as a function of r and arrange the bulk strain in the second example to produce the same answer.
Our case is Changing coordinates to u = a 11 + a 22 , v = a 11 − a 22 , w = 2a 12 we see that A 1 , in these coordinates, is u = 3, v = cos 2s, w = sin 2s, that is a circle in the u = 3 plane for s ranging from 0 to π. Calculation of g is now the integral over the plane n θ •(a 11 , a 12 /2, a 22 ) = r −2 of the generalized function δ(u−3)δ( √ v 2 + w 2 −1).In the new coordinates U = (u, v, w) the plane under consideration in ( 14) is where the vector on the left has unit length.This intersects the u = 3 plane in a family of lines a distance 2 √ 2r −2 − 3 from the u.The Radon transform of the unit density on the unit circle in (v, w) space (see Appendix A) To evaluate g we need to take account of the Jacobian for the change of variables to (u, v, w) and the direction cosine of 1/ √ 2 of the plane to the u = 3 plane giving where r 0 = 2/ 2 + 3 √ 2, r 1 = 2/ 4 + 3 √ 2. We simply have to construct an A 2 (s) with a density supported on the u axis that produces the same g.Each plane in (14) intersects the u axis at an angle π/4, at the point U = (2 √ 2r −2 , 0, 0) so we just assign the density φ 2 = g 1 (u 1/2 2 −3/4 )δ(v)δ(w).To find an example of an A 2 (s) with this density we simply form the cumulative density in the u direction and use this as the bulk strain in the transverse direction, with a rescaling so that s takes values in the same interval as for A 1 .
One can repeat this construction for any pairs of curves with the property that any plane with normal in the cone passes through one point on each curve an by choice of the density function make the diffraction pattern the same.
In general the support of the joint distribution of a smooth transverse strain along a line will be an immersed curve in three dimensional space.When the planes defined above intersect this curve non-transversely a singularity will result in the marginal distribution that is given by any of the plane transforms in the data.In the case detailed above a singularity occurs when the plane is tangent to the circle for the A 1 case and at the ends of the range on the A 2 case.
Clearly we can recover the first moment and hence the TRT data as detailed in [17].We can also recover the distribution of the normal component to a plane from the diffraction data for that plane.To see this suppose a 11 is the normal component.The directions (1, 0, 0) satisfies the cone condition in the frequency domain, and so the Fourier transform along this direction determines the marginal distribution for a 11 .Thus all results for the histogram Radon transform can be used plane by plane with this data.Similarly we know the marginal distribution of a 22 .It is not yet clear how to use the additional data in an efficient way.In particular we would prefer a reconstruction method with fewer rotation axes than the three given in [8] for a general strain, or two for a potential (elastic) strain.We have exhibited a wide range of tomographic problems in which a distribution can be measured for each ray.While in general fitting a distribution (or in the discrete case a histogram) of the data is generally a non-linear problems we saw that in the case of a scalar unknown both the moment data and the bins in the cumulative histogram reduce our problem known linear tomographic problems.Further work is needed on practical useful algorithms for scalar histogram tomography.For near infrared spectral tomography of chemical species, attention needs to be directed to the recovery of multiple parameters including the density of the chemical species and the pressure as well as temperature.
For a vector field, we were able to give for the first time an explicit reconstruction for the potential case, using the second moment.The general case is an obvious candidate for further study, and especially the problem of limited data.
In the case of Bragg edge strain tomography we showed how the histogram could be extracted from spectral data with a high resolution in energy near an edge.We were not so lucky with the theory as in the vector potential case, as we were only able to find a non-linear system of PDEs for the unknown displacement field derived from moment data.It still remains to be seen in practice what distributional data can be recovered from neutron spectra.
In x-ray diffraction tomography of strain we were able to give an example of strain distributions resulting in the same diffraction pattern for one ray, so that in this case we do not recover the histogram transverse ray transform.We have yet to see experimental data demonstrating what can be recovered from x-ray diffraction data.
We did not treat in this paper the case of strain tomography from diffraction data from a single crystal.In contrast to the polycrystalline case the unstrained diffraction pattern for a pencil beam is a two dimensional lattice of small diffraction spots, and the lattice is distorted by a strain in the transverse direction.One might expect under ideal conditions and small strain to see each small spot in the lattice spread to a density function when the strain varies along the beam.The mean of these distributions gives the average displacement of the spot in reciprocal space.From this one would expect to be able to derive the TRT data of the strain.As the strain is potential, rotations about two axes are sufficient for reconstruction [8].A practical case where this problem arises is electron strain measurement in silicon for electronic devices.The measurement technique is already established [16], however it is unlikely that measurements can be taken for a complete rotation in a sample of practical significance.This leads again to the consideration of reconstruction fro ma limited range of projections.Like the polycrystalline case, this technique would not yield the full HTRT.Each spot gives a marginal distribution of a vector, but it is not clear what information about the joint distribution of the transverse strain tensor can be recovered.
In many practical cases it is likely that fitting algorithms to the non-linear problem using all available data, rather than simple explicit reconstruction formulae.Nevertheless theoretical methods at least give a useful insight in to the sufficiency or insufficiency of specific sets of measurements.
In practical problems there is always the response of the instrumentation and blurring effects from other physical phenomena to consider.One can generally expect the distribution that is measured to be convolved with a point spread function and some calibration and deconvolution is required to give a good estimate of a histogram.As histogram tomography is investigated experimentally this is one of the challenges that must me overcome.In many cases the first step will be to validate and calibrate the measured histogram from a sample with a known distribution.
The more complex cases have forward problems involving restriction, distributions, linear projections and integrals.There are likely also to problems where a weighted integral or histogram is taken, for example when scattered radiation is measured but it is also attenuated.Further work could look at a general theoretical framework that includes many practical cases.This might be expected to draw further on work from probability theory on recovering a joint distribution from marginal and conditional data.
We look forward to practical realisations of these techniques as well as further theoretical developments and hope this paper will inspire such work.

Figure 1 :
Figure 1: Illustration of distribution and cumulative distribution of a smooth function f .The distributions are shown with the frequency on the horizontal axes for easy comparison with the graph of f .Note that critical points of f along the line result in singularities of the distribution.In general the distribution gives how often each value in the range occurs, but not where they occur.