Open Access
Add to BibSonomyAbstract
Quite recently, a method has been presented to reconstruct X-ray scattering tensors from projections obtained in a grating interferometry setup. The original publications present a rather specialised approach, for instance by suggesting a single SART-based solver. In this work, we propose a novel approach to solving the inverse problem, allowing the use of other algorithms than SART (like conjugate gradient), a faster tensor recovery, and an intuitive visualisation. Furthermore, we introduce constraint enforcement for X-ray tensor tomography (cXTT) and demonstrate that this yields visually smoother results in comparison to the state-of-art approach, similar to regularisation.
© 2015 Optical Society of America
1. Introduction and background
Most applications of X-ray imaging limit themselves to attenuation, providing images with strong contrast between highly absorbing structures such as bones or artificially contrasted arteries, and lowly absorbing background such as tissue. From a physical point of view, however, more information could be acquired, for instance by means of X-ray interferometry. For decades, this technique was limited to laboratory setups, as it requires highly coherent X-ray illumination, and thus a synchrotron or an equivalent source.
Only some few years ago, Pfeiffer et al. [1] have demonstrated a novel setup for X-ray grating interferometry: As shown in Fig. 1, an ‘ordinary’ setup of X-ray tube (T), specimen (S) and detector (D) has been extended by inserting a source grating (G0) after the tube, making the incoming illumination partly coherent, and two more gratings (G1, G2) forming the interferometer itself. Instead of a single image per orientation/angle as would be done in ‘classical’ computed tomography, multiple images are acquired while interferometry grating G1 is shifted sideways, thus producing several images of the same situation with different interference patterns. Based on such images, three signal components can be extracted per perspective: An absorption image as commonly used in X-ray computed tomography, a phase-contrast and a dark-field component [2]. Example images are given in Fig. 2.
Fig. 1 Sketch of an X-ray grating interferometry setup. X-ray tube (T), source grating (G0), shifting interferometer grating (G1), specimen (S), static interferometer grating (G2), and detector (D).
Fig. 2 Three components of a sample projection of a tooth: Absorption, phase contrast and dark field. The contrast of these images has been manually improved for better visibility, and the images have been cropped. The structure in the lower-left quadrant is the sample holder.
Despite certain problems such as phase wrapping, tomographic phase-contrast reconstructions have since become quite popular, particularly due to their improved soft-tissue contrast, and improving them is an active field of research.
A straight-forward three-dimensional tomographic reconstruction of the dark-field signal yields a volumetric representation of the respective X-ray scattering magnitudes. Such information can be of practical value when isotropic scattering is of interest. Previous studies have shown its potential for lung and breast imaging [3–5], micro-bubble contrast agents [6] and material testing [7]. In all these cases, dark-field reconstructions have been shown to provide contrast where the other two signals, phase-contrast and absorption, only yield poor results at best.
Tomographic reconstruction of the anisotropic component of the dark-field signal on the other side has not yet gained much attention. Only recently, Malecki et al. [8] have presented tensor-valued tomographic reconstruction, thus taking directional information into account. In addition to the usual tomographic axis, the sample is also rotated around the other two axes using an Euler cradle. The authors present a first mathematical model and a specialised SART-based reconstruction algorithm, thus obtaining the first scattering tensor reconstructions. Informally speaking, they obtain an ellipsoid at every voxel, and its size and shape hints at the structure of the specimen at this location. In particular, plate-like scattering ellipsoids hint at fibre- or tube-like structures, as sketched in Fig. 3.
Fig. 3 Expected X-ray scattering (blue ellipsoids) at fibre- or tube-like structures (grey). Vice versa, we will interpret the smallest half-axes of reconstructed scattering ellipsoids as fibre directions for visualisation.
The only other work on scattering tensor reconstruction has recently been presented by Bayer et al. [9]. In contrast to the aforementioned approach, this group uses the single tomographic axis of rotation only and relies on extended mathematical modelling. In particular, they apparently do not reconstruct full scattering tensors, but projections of fibre directions onto the tomographic plane, and thus a stack of two-dimensional projections rather than a volume of three-dimensional shapes. Apart from X-rays, another well-known tensor-based imaging modality is diffusion tensor magnetic resonance imaging (DTI) [10] – which only shares visualisation techniques.
In this paper, we use the approach of Malecki et al., and present several algorithmic improvements. Among them are, most importantly, a novel, more generic description for solving X-ray tensor tomography (XTT) problems [11], and two ways to enforce tensor-shapes during reconstruction.
2. Methods
Scattering is not a scalar entity (as X-ray attenuation is), but a tensor-valued one, that is a three-dimensional shape per location. This requires a more complicated mathematical model.
Malecki et al. [8] propose to consider a finite set of K pre-defined, normalised, well-distributed sampling directions ε̂k ∈ ℝ3 rather than a tensor. Here, and in the remainder of the paper, we use the ‘hat’ to denote normalised vectors, that is ‖ε̂k‖2 = 1 ∀k = 1,...,K. For each of these ε̂k and every voxel xi, a corresponding scattering coefficient ζk(xi) ∈ ℝ will be reconstructed, and the tensor itself can later be approximated by voxel-wise fitting of tensor descriptions to these ‘bouquets’ of weighted directions. This entire process is sketched for a single voxel in Fig. 4, a visual account of a larger region is given in Fig. 5. Note that the exact orientation of the sampling directions is rather arbitrary, and the choice does not require any knowledge about the expected reconstruction outcome. The sampling directions are just auxiliary, virtual entities, allowing to avoid an explicit description of the scattering tensor at a given location by considering its projection onto a finite set of well-distributed reference directions instead. They are chosen from a single hemisphere only, as the scattering tensors are supposed to be symmetric.
Fig. 4 Sampling directions, and ellipsoid fitting for a single voxel. Usually, we use K = 13 directions for reconstruction, but limit ourselves to K = 7 in this sketch, for clarity. They consist of the standard base vectors – red, green, blue in (a) for x, y, z, respectively – and diagonals – black. For a given voxel x, reconstruction yields a scattering coefficient ζk(x) for every sampling direction ε̂k. These are indicated by bold black marks in (b), and we mirror them along the negative sampling direction, yielding the small black marks. Finally, an ellipsoid can be fitted to that scaled, mirrored ‘bouquet’ afterwards, see (c).
Fig. 5 Different variants of the same data, a downsampled detail from an actual result. Reconstruction is performed in terms of scaling the sampling directions (a), the scattering tensors are obtained by retrospectively fitting ellipsoids in a voxel-wise fashion (b). The colours indicate the directions of the smallest half-axes, see Fig. 8(a).
Based on this approach, Malecki et al. derive a forward model for X-ray scattering [12–14], resembling the Beer-Lambert law describing X-ray absorption:
Here, dj ∈ ℝ denotes the j’th scalar dark-field measurement, and Lj the corresponding ray with normalised direction Îj ∈ ℝ3. Note that despite acquiring two-dimensional projections, we treat the individual pixels as independent measurements with corresponding rays. Index j = 1,..., J extends over all pixels of all projection images hence. The normalised sensitivity direction t̂j ∈ ℝ3 is orthogonal to the grating lines and parallel to the surfaces of the – mutually parallel – gratings, and thus depends on the perspective during measurement j.Considering Eq. (1) in more detail, the true scattering measurement at every voxel x along the ray is generally modelled as finite sum of K measurements along the pre-defined sampling directions ε̂k. For every such summand, two things need to be taken into account: Scattering must be possible to occur at all, considering the direction Îj of the incoming X-ray with respect to ε̂k. This fact is modelled by the magnitude of the cross-product | · × · |, thus excluding head-on views. Furthermore, if scattering occurs, only components along sensitivity direction t̂j contribute to the measurement as modelled by the scalar product 〈·,· 〉 [13, 14]. Vice versa, in order to have a good sampling of the scattering tensors, more than just the usual tomographic scanning perspectives are taken into consideration while scanning, see Fig. 6.
Fig. 6 Viewing directions and Euler cradle. Only the perspectives marked by red lines are in the ‘normal’ viewing plane as used in standard CT applications. Blue lines mark additional off-plane perspectives. The dashed black line indicates the trajectory of the X-ray source-detector “camera” around the green cube, representing a specimen. The coverage gaps are due to limitations imposed by the Euler cradle.
Returning to Eq. (1), factoring out the (unknown) squared scattering coefficients ηk(x) := ζk(x)2 from the squared scalar product and defining the weight factor
yields a formulation very similar to the well-known Radon transform: Note that the factors defined in Eq. (2) are independent of the unknowns, and can thus be precomputed for more efficient reconstruction.2.1. Reconstruction
Obviously, Eq. (4) is a weighted sum of line integrals. In particular, in a discretised setting, it can be rewritten as scalar product
where ηk ∈ ℝI is the vector of the squared coefficients for the kth sampling direction for all voxels (I is the number of voxels). aj denotes the system matrix row for measurement j, and is thus part of the ‘standard’ matrix as used for an equivalent attenuation reconstruction. This vector contains geometric information about the arrangement of X-ray source, specimen and sensor during measurement j, an information already included abstractly in Eq. (4) as ray Lj. Let A = (aj), A ∈ ℝJ×I, denote the entire system matrix describing all J line integrals. Furthermore, let Dk = diag(vk1, vk2, ...) denote a diagonal scaling matrix containing the weighting factors from Eq. (2) for sampling direction k. Then, using Eq. (5) and defining the measurement vector m = (mj) = (−ln dj), a huge linear system can be derived: This linear system has K times as many unknowns than the corresponding system for computing a ‘traditional’ tomographic attenuation reconstruction: The original system matrix A is of size J × I, and H is of size J × IK. In practical settings, that is when reconstructing a sufficiently sized three-dimensional volume, the system matrix A is already far too large to fit into computer memory. Instead of handling the matrix directly, projector software is typically used that simulates the ray, thus computing the entries of A on the fly [15].When reconstructing scattering tensors, it is desirable to use existing software infrastructure that has been put into place for other tomographic reconstruction scenarios rather than solving Eq. (10) directly. Consequently, in their original work, Malecki et al. suggest a specially crafted variant of SART [16, 17] for reconstructing the unknown squared coefficients ηk.
Using Eq. (5), however, a much more generic approach can be taken, supporting arbitrary iterative solvers such as particularly the better-behaving method of Conjugate Gradients (CG) [18,19]. Let x1 = approximate(A, b, x0) denote an auxiliary function running a single iteration of an iterative linear solver such as CG or SART, thus very approximately ‘solving’ the linear system Ax = b. Then, an iterative algorithm solving the tensor reconstruction problem (5) can be defined: A single iteration q essentially consists of ‘approximately solving’ K modified linear systems
for each sampling direction k using function approximate as previously defined, where right-hand side vector denotes a reduced measurement vector relevant for sampling direction k, based on the estimates of the previous iteration, and an intermediate, temporary vector. The latter is then used to compute the current iterate via relaxation as Altogether, this proposed approach can be thought of as simultaneously solving K linear systems in an interleaved sense. This entire proposed generic tensor reconstruction approach is given in alg. 1. Note that the scattering coefficients need to be extracted after reconstruction of all ηk by component-wise application of the square root: Using the absolute magnitude of ηk is a valid security measure, considering the symmetry of the tensors.
Alg. 1. Generic tomographic X-ray tensor reconstruction. A denotes the system matrix describing the imaging process, Dk a scaling matrix containing weighting factors as defined in Eq. (2). m is the measurement vector, and a vector containing the qth iterate of the voxel-wise squared scattering coefficients corresponding to sampling direction k. approximate is a function running a single iteration of an arbitrary iterative linear solver.
2.2. Ellipsoid Fitting
After reconstruction, that is after executing several iterations as described in section 2.1, voxel-wise scattering coefficients ζk for the sampling vectors ε̂k have been recovered. Malecki et al. propose to fit ellipsoids to these weighted vectors in order to obtain voxel-wise tensors. In particular, they propose to use an iterative ellipsoid fitter [20], apparently intended for rather degenerate cases where ellipsoids need to be matched to just a couple of ill-distributed sample points.
However, as shown in Fig. 4, the sampling locations are well-distributed, and we propose to use principal component analysis [21–23] (PCA) instead of the iterative approach, thus saving considerable computation time.
In more detail, at every voxel xi, a set of 2K direction vectors
can be defined by scaling the normalised sampling directions ε̂k with the positive and negative reconstructed corresponding coefficients ζk(xi) = |ηk(xi)|1/2. Again, this is due to the symmetry of the scattering tensors. Note that the mean of this set trivially equals 0. Let Ci ∈ ℝ3×3 denote the covariance matrix of Si. Then, its eigen-decomposition yields a diagonal matrix Λi ∈ ℝ3×3 containing the three eigenvalues λi,1, λi,2, λi,3 of the set’s covariance matrix Ci, and an orthogonal matrix Vi ∈ O(3) ⊂ ℝ3×3 containing the corresponding mutually orthonormal eigenvectors vi,1, vi,2, vi,3 as column vectors.Let η̄i = ∑k |ηk(xi)|/K denote the average squared scattering magnitude at voxel xi, and λ̄i = (|λi,1| + |λi,2| + |λi,3|)/3 the corresponding average eigenvalue magnitude. Define a size correction factor σi := η̄i/λ̄i for scaling the statistically defined ellipsoid to the point set S. Then, the scattering ellipsoid at voxel xi is defined by half-axis lengths ri,1 = [σiλi,1]1/2, ri,2 = [σiλi,2]1/2, ri,3 = [σiλi,3]1/2 with respect to the orthonormal basis formed by vectors vi,1, vi,2, vi,3.
2.3. Constraint enforcement
As it will turn out in the experiments, unconstrained reconstruction of the scattering coefficients yields useful but noisy results. Instead of just fitting ellipsoids retrospectively, it makes sense to enforce ellipsoidal shapes during reconstruction. Therefore, we propose to post-process the scattering coefficients at the end of every iteration as indicated in the last comment of alg. 1, thus forcing them to evolve in the vicinity of the manifold of ellipsoid-shaped tensors.
2.3.1. Hard ellipsoid constraint
An obvious approach to reach this aim is to fit ellipsoids after every iteration, and to project the reconstructed squared coefficients ηk onto them. We refer to this method as hard ellipsoid constraint as the coefficients will really be forced into ellipsoidal shapes.
In detail, after the computation of all squared coefficients in iteration q, we visit every voxel xi and fit an ellipsoid as described in section 2.2, thus obtaining a coordinate frame Vi = [vi,1, vi,2, vi,3] ∈ ℝ3×3 and corresponding half-axes ri,1, ri,2, ri,3 ∈ ℝ. Then, the reconstructed coefficients are replaced with the projections of the respective normalised sampling directions ε̂k onto this ellipsoid.
In detail, to do this for every sampling direction k and voxel xi, we first rotate the (by design normalised) vector into the ellipsoid’s coordinate frame, thus obtaining a unit vector . By definition, the scaled vector σi,k · [xi,k, yi,k, zi,k]T, σi,k ∈ ℝ, resides on the surface of the ellipsoid if
This Eq. can simply be solved for σi,k2, and considering that squared coefficients are reconstructed, the correctly projected coefficients are for all voxels xi and directions k. See Fig. 7(a) for an illustrative sketch.Fig. 7 Constraint enforcement. In both cases, the reconstructed coefficients are forced to be close to the manifold of valid ellipsoids while iterating. The hard constraint (a) projects back onto ellipsoids directly, where a – in this toy example exaggerated – scattering coefficient (dotted black line) is shortened appropriately (solid black line). The soft constraint smoothes the coefficients with respect to the other sampling directions of the same voxel respecting the angular relation, thus favouring ellipsoids in a relaxed sense. This is done using a Gaussian smoothing kernel based on the scalar product as shown in (b), ranging between 0 (blue) and 1 (yellow), before normalisation.
2.3.2. Soft ellipsoid constraint
Alternatively, we propose to use a softer variant favouring ellipsoids but allowing some more freedom. Essentially, this is done by smoothing the K direction coefficients per voxel, over the ellipsoid. We refer to this approach as soft ellipsoid constraint. Note that we still only consider individual voxels and their K coefficients; we do not take a neighbourhood into account as one would do for smoothness constraints in regularised X-ray attenuation reconstruction.
In detail, we first compute the pair-wise absolute scalar products between all K normalised sampling directions, and obtain a matrix S ∈ [0, 1]K×K, Sk,l = |〈ε̂k, ε̂l〉|. In order to increase the influence of more similar sampling directions, we additionally apply a Gaussian and obtain a matrix G ∈ ℝK×K, Gk,l = exp[−(Sk,l − 1)2/2μ] with some selectable mean μ. The larger μ is chosen, the higher will be the influence of the smoothing.
To process coefficient k, a matching convolution kernel can be defined by normalising the kth row of G to a sum of 1, denoted as , . Collecting all K coefficients into a vector, the regularised kth value for voxel xi computes as
In other words, this approach can be considered as angle-dependent Gaussian smoothing of the squared coefficients ηk(xi) over all directions k for every individual voxel xi. This constraint is illustrated in Fig. 7(b).2.4. Visualisation
A last, vital point for X-ray tensor tomography is a proper visualisation of the results. This is of particular importance as the voxel-wise tensors can not be accurately shown using rendering techniques intended for scalar-valued data. Apart from plotting the ellipsoids themselves, a prominent option is to visualise derived shapes.
Interested in fibre- and tube-like structures, we therefore propose a streamline visualisation based on the classical Runge-Kutta method [24–26] (RK4). This metaphor is well-known in the field of diffusion-tensor magnetic resonance imaging (DTI) [10,27], and commonly referred to as tractography. As sketched in Fig. 3, plate-like ellipsoids hint at fibres in the direction of the smallest half-axes. We employ RK4 to trace streamlines along these smallest half-axes, adapted to account for the symmetric, ‘bidirectional’ nature of ellipsoid axes (compared to a vector field). This is done by inverting directions queried from the raw ellipsoid axes whenever the scalar product with the direction of an incoming fibre is negative, i.e. new directions are chosen to continue the streamline in forward direction. Furthermore, we employ colour coding to give cues about the orientation of these fibres. An example is shown in Fig. 8.
Fig. 8 Colour coding of streamline visualisation. We use colours sampled from a sphere (a) to give visual cues about the orientation of the ‘fibres’. Note that the colour ball is symmetric with respect to the x-y-, x-z-, and y-z-planes. In a sample picture of the carbon knot (b), looking down from a slightly elevated perspective, vertically (green), horizontally (red) and obliquely (orange and yellow) oriented parts can be easily distinguished.
Note that the streamlines do not necessarily correspond to real fibres, and there is no intention for this to be the case. The streamlines are to be interpreted as scattering visualisation, however, using a model which is likely similar to the invisible structure.
3. Experiments
We have implemented the proposed reconstruction algorithm within CampRecon [28], a C++ framework for solving linear inverse problems. The discrete integration along lines Ax and its adjoint AT b are approximated roughly via intersections between ray and pixels as proposed by Siddon [29], using the GPU-based X-ray projector developed by Fehringer et al. [15] in OpenCL. All experiments have been run on a computer with dual Intel Xeon E5-2687W processors and a Nvidia Tesla K20 accelerator.
In general, we have computed reconstructions of three datasets, the carbon ‘knot’, a knotted bunch of carbon fibres embedded in hot glue, the tree ‘branch’, a short piece of raw wood, and a ‘tooth’. Photographies of the samples are given in Fig. 9. X-ray images of all samples have been acquired in our experimental setup, see Fig. 10. We use 732 projections of 321 × 321 pixels for the knot, 551 projections of 301 × 301 pixels for the branch, and 902 projections of 701 × 701 pixels for the tooth, with trajectories resembling the one in Fig. 6. The projections in Fig. 2 have been taken from the tooth dataset.
Fig. 10 Wide-angle view of the actual setup. From left to right: X-ray source (T) with source grating (G0) directly in front of it, at the center the phase grating (G1), then the Euler cradle with sample (S), and adjacent to it the analyser grating (G2, in the cross-shaped holder); behind the latter the X-ray detector (D, hidden).
All samples were measured at an acceleration voltage of 60 kVp. Eight phase steps were recorded per projection with a flat panel detector (Varian) with pixel pitch of 127 μm. The exposure time was 1 s per phase step. A π/2 phase grating for the design energy of 45 kV composed of 8 μm high Ni lines and period of 5 μm was used. The two other gratings were absorption gratings with 170 μm high Au lines and a period of 10 μm. The interferometer was symmetric with both inter-grating distances being 92.7 cm.
Each of the datasets was reconstructed without constraint enforcement, with hard, and with soft ellipsoid enforcement, each time with K = 13 sampling directions. In all cases, we executed 100 iterations of alg. 1 and do not employ any other stopping criterion. We use a single iteration of CG as function approximate. Note that this corresponds to a single Landweber step with respect to the sub-problem. For the soft constraint, we generally use an experimentally established smoothing mean of μ = 0.1, unless stated otherwise. Reconstruction times are about 1.75 hours for the branch (2513 voxels), 2 hours for the knot (2013 voxels), and 15.5 hours for the tooth (301 × 501 × 291 voxels). Relating the number of voxels with the input projections, it is clear that the linear system m = Hs as defined in Eq. (10) is underdetermined in all our experiments. The ratios of measurements to unknowns are 24.3 % for the branch, 70.4 % for the knot, and 77.7 % for the tooth, respectively.
Furthermore, we have obtained a reconstruction of the knot produced with the original SART-variant proposed by Malecki et al., for comparison.
4. Results
4.1. Numerical behaviour
In order to check the behaviour of the three algorithm variants (unconstrained, soft, and hard ellipsoid constraint), we have computed the normalised residual norms
for iterations q ∈ {1,...,100}, and also the normalised mean updates Sample plots of these sequences for the knot are given in Fig. 11. There, the curves for both measures flatten out with increasing number of iterations. The unconstrained variant tends to reach a smaller residual norm, while the update norms appear to oscillate. Vice versa, the constrained variants show larger residual norms, but smooth updates. For other datasets, we obtain similar curves.Fig. 11 Behaviour of the algorithm for the knot. The plots show the normalised residual norm r(q) as defined in Eq. (21) over iterations q (left) and the normalised mean update Δ(q) as defined in Eq. (22) over iterations q (right). The unconstrained version yields smaller residuals, but the updates are noisy.
4.2. Comparison with state of art
Next, in order to compare our method with the original SART-based method of Malecki et al., we compare an unconstrained reconstruction of the knot with a reference result, see Fig. 12. As Malecki’s results can not be considered as ground truth, we restrict ourselves to a visual comparison.
Fig. 12 Comparison of attenuation reconstruction (a), Malecki’s tensor reconstruction (b), and our (unconstrained) tensor reconstruction (c). The first image is included for reference, to show that scattering data is of a considerably different nature than usual attenuation reconstructions. The two scattering reconstructions are largely equivalent, considering the different algorithms with different parameters.
4.3. Knot
Figure 13 shows volume renderings of raw scattering coefficients relating to a single sampling direction, but for the three different constraint enforcement schemes. Most importantly, the unconstrained version shows considerable streak artefacts, and the constrained versions do not.
Fig. 13 Volume renderings of the knot’s scattering coefficients for sampling direction ε̂k = [0, 0, 1]T. The unconstrained reconstruction (a) shows strong streak artefacts, the two constrained versions (b) and (c) are much clearer.
Figure 14 and Media 1 show streamline visualisations of the knot. As can be seen, enforcing the two constraints yield visually smoother fibres that are more densely packed. Again, note that these fibres visualise the scattering ellipsoids, and are not to be considered reconstructions of the raw carbon fibres.
Fig. 14 Streamline visualisation of the knot’s scattering ellipsoids. (Also see Media 1.) The streamlines are supposed to follow the directions of the carbon fibres, see Fig. 9(a), but are not intended to accurately reconstruct individual fibres. The unconstrained version (a) shows noise, the two constrained versions (b) and (c) are visually considerably smoother. The ‘waves’ in the lower right appear to be additional scattering caused by the sample holder.
4.4. Branch
Figure 15 and Media 2 show streamline visualisations of the branch. For this dataset, enforcing ellipsoid constraints has even more prominent effects. Considering the structure of wood, scattering will be caused primarily by the tiny vessels transporting water towards the leaves. That is, streamlines are supposed to run mainly in parallel to the axes of the branches. In the unconstrained case, reconstruction of scattering along the main branch fails almost entirely, but the constrained versions are able to recover more reasonable scattering streamlines there.
Fig. 15 Streamline visualisation of the branch. (Also see Media 2.) Scattering is supposed to be caused by the tiny tubular vessels embedded in wood that transport water towards the leaves. The streamlines are thus supposed to run mainly in parallel to the individual branches. The unconstrained version (a) fails to recover useful scattering along the main branch. The lightly constrained versions (b) and (c) show more reasonable scattering there. The hard constraint overshoots and produces wavy patterns along the main branch, and introduces a ‘combing’ effect for the other branches. The latter can already be seen in (c).
We also modify the smoothing parameter μ for the soft constraint. In particular, a ‘combing’ effect can be observed for stronger constraint enforcement where streamlines are nicely aligned to each other, but ‘combed’ away from their expected location, see the side branches in Fig. 15(d). On the other side, a smaller value of μ will increase the amount of noise but cause the reconstruction to stay closer to the data. Especially the side branches appear to benefit from a smaller parameter, as the ‘combing’ effect is considerably reduced, and the main branch still reconstructs properly.
4.5. Tooth
The reconstruction results of the tooth are given in Fig. 16 and in Media 3. This sample is particularly intriguing: Teeth generally consist of a hard crown of highly mineralised enamel, and a root covered by cementum. Below these layers, the core of the tooth consists of dentine and, embedded within, the pulp chamber containing living tissue, blood vessels and nerves. Dentine is a fibrous material and less mineralised than enamel. Tiny dentinal tubules (about 2 μm in diameter) are passing through it in radial direction, between pulp chamber and the surface (but not through the enamel). These tubules are generally not visible in X-ray attenuation reconstructions of usual resolution, due to their small size.
Fig. 16 Reconstruction results of the tooth. (Also see Media 3.) Volume rendering of X-ray attenuation (a) showing enamel (white) and dentine (gray). Equivalent attenuation volume rendering showing dentine only (b). Scattering streamlines obtained without constraint (c), with soft constraint with μ = 0.08 (d) and μ = 0.10 (e), and with hard constraint (f). Scattering is caused by dentinal tubules, tiny structures in radial direction, as indicated by the streamlines. The pulp chamber was obtained by masking.
Despite this, however, the scattering caused by them can be measured and reconstructed. Note that the results in Fig. 16 have been masked using the reconstruction of the attenuation signal. This is necessary as the contents of the pulp chamber, relicts of the tissue, cause considerable isotropic scattering themselves – clearly visible as distinct black region in the scattering projection in Fig. 2(c).
The influence of constraint enforcement is similar to the one observed for the other samples. However, at the lower parts of the roots, the soft constraint also yields streamlines caused by pulp chamber scattering.
5. Discussion
As shown in the first experiments, the proposed reconstruction algorithm apparently shows reasonable behaviour in our experiments. Plots similar to the ones shown in Fig. 11 can be produced for all experiments that we have conducted so far. Note that these plots only provide hints that the iteration does not become unstable, they neither indicate image quality nor do they prove convergence.
Comparing with the state of art, the unconstrained variant of our algorithm yields a result quite similar to Malecki’s original SART-modification, considering the different nature and parameters of the two methods.
The constrained variants of our method have been shown to produce considerably less artefacts and visually smoother, non-oscillating updates. As most notable example, the main branch of the wood sample inherently requires some kind of constraint enforcement, to push the reconstruction towards a reasonable result. On the other side, exaggerated constraint enforcement may lead to new artefacts, particularly unexpected wavy patterns or ‘combed’ directions, as can also be seen for the wood sample, for instance. It appears that constraint enforcement is necessary, but its influence must be limited. From our point of view, soft constraint enforcement with a conservatively chosen parameter μ is the most promising approach.
In general, a remarkable feature of X-ray tensor tomography is that scattering effects caused by sub-voxel-sized structures can be recovered, thus giving insight into the microstructure. In our case, this holds for both, the vessels in the branch, and the dentinal tubules in the tooth. Of course, important future work will be an exact correlation of X-ray scattering results with suitable high-resolution pictures of the samples, through μCT or microscopy, for instance.
An important consequence of our work is that the assumption of pure ellipsoids appears not to be entirely correct, considering the deterioration with stronger constraint enforcement. Based on sampling directions, our reconstruction is agnostic of that assumption, but future work will need to look into more appropriate scattering descriptions. Furthermore, even after flattening out, the residual norm is still rather large. Usually, this would hint at a mathematical problem. Considering the visual quality of the reconstructions and the curve of the update norm, however, this may also be caused by Eq. (1) only partially describing the dark-field signal. Both (highly related) questions will need to be subject to future investigation from a physics perspective.
Considering the proposed algorithm, many interesting mathematical questions are still open and will need to be tackled in future work. As pointed out by a reviewer, it should be possible to position the proposed algorithm in the context of other, well-established solution algorithms, such as column-based block iterative algorithms [17], for the huge system defined in Eq. (10). Particularly interesting will be an investigation of two points: First, our algorithm computes the weighted forward projections WFPk based on the previous iterates, and thus exclusively uses the previous state to update the solution. This decision has been made from an engineering prespective, and to avoid bias towards the scattering coefficients that are processed earlier. From a mathematical point of view, however, it may be beneficial to only use the previous state where necessary and the already known current state where possible. The second question is the dependence on the exact choice of function approximate. In this context, it will also be useful to attempt to solve the original linear system directly. More advanced mathematical topics include an investigation of the convergence, and the influence of the constraint enforcement schemes onto it. Finally, another important open question is inter-voxel regularisation such as total variation minimisation.
6. Conclusion
We have presented a new processing chain for X-ray tensor tomography. In particular, we have detailed a generic formulation for solving the inverse problem using arbitrary linear solvers, firmly based on the forward model and avoiding the tight coupling to SART, along with constraint enforcement schemes improving the reconstructions considerably. Furthermore, we have presented a way for quick recovery of the ellipsoids, and we have introduced streamline visualisation to XTT.
A lot of work still needs to go into XTT. At this point, the mechanical design of the setup requires small samples, and the radiation exposure is considerable. The latter is due to both the large number of acquisitions (more perspectives, more acquisitions per perspective), and the additional absorption caused by the interferometer grating between sample and detector, thus requiring higher radiation at the specimen in order to obtain good detector readings. Further, more theoretical open questions have already been outlined above. Nevertheless, application in ex-vivo imaging or material testing is well within reach today, and improvements are subject of active research.
Acknowledgments
This work was partially funded by the DFG Cluster of Excellence Munich–Centre for Advanced Photonics (MAP) and carried out with the support of the Karlsruhe Nano Micro Facility (KNMF), a Helmholtz Research Infrastructure at Karlsruhe Institute of Technology (KIT). J. V. thanks Andrea Porsche for skilfully recovering the tooth sample. F. S. is funded by the TUM Graduate School.
We are indebted to the anonymous reviewer who provided an exceptional report including many valuable comments about mathematical and numerical aspects of the proposed reconstruction approach.
References and links
1. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nat. Phys. 2, 258–261 (2006).
2. F. Pfeiffer, M. Bech, O. Bunk, P. Kraft, E. F. Eikenberry, C. Brönnimann, C. Grünzweig, and C. David, “Hard-X-ray dark-field imaging using a grating interferometer,” Nat. Mater. 7, 134–137 (2008). [PubMed]
3. A. Yaroshenko, F. G. Meinel, M. Bech, A. Tapfer, A. Velroyen, S. Schleede, S. Auweter, A. Bohla, A. Ö. Yildirim, K. Nikolaou, F. Bamberg, O. Eickelberg, M. F. Reiser, and F. Pfeiffer, “Pulmonary emphysema diagnosis with a preclinical small-animal X-ray dark-field scatter-contrast scanner,” Radiol. 269, 427–433 (2013).
4. M. Ando, K. Yamasaki, F. Toyofuku, H. Sugiyama, C. Ohbayashi, G. Li, L. Pan, X. Jiang, W. Pattanasiriwisawa, D. Shimao, E. Hashimoto, T. Kimura, M. Tsuneyoshi, E. Ueno, K. Tokumori, A. Maksimenko, Y. Higashida, and M. Hirano, “Attempt at visualizing breast cancer with X-ray dark field imaging,” Jpn. J. Appl. Phys. 44, L528–L531 (2005).
5. S. Sidhu, G. Falzon, S. A. Hart, J. G. Fox, R. A. Lewis, and K. K. W. Siu, “Classification of breast tissue using a laboratory system for small-angle x-ray scattering (SAXS),” Phys. Med. Biol. 56, 6779–6791 (2011). [PubMed]
6. F. Arfelli, L. Rigon, and R. H. Menk, “Microbubbles as X-ray scattering contrast agents using analyzer-based imaging,” Phys. Med. Biol. 55, 1643–1658 (2010). [PubMed]
7. V. Revol, I. Jerjen, C. Kottler, P. Schütz, R. Kaufmann, T. Lüthi, U. Sennhauser, U. Straumann, and C. Urban, “Sub-pixel porosity revealed by X-ray scatter dark field imaging,” J. Appl. Phys. 110, 044912 (2011).
8. A. Malecki, G. Potdevin, T. Biernath, E. Eggl, K. Willer, T. Lasser, J. Maisenbacher, J. Gibmeier, A. Wanner, and F. Pfeiffer, “X-ray tensor tomography,” EPL105 (2014).
9. F. L. Bayer, S. Hu, A. Maier, T. Weber, G. Anton, T. Michel, and C. P. Riess, “Reconstruction of scalar and vectorial components in X-ray dark-field tomography,” Proc. Natl. Acad. Sci. U.S.A. 111, 12699–12704 (2014). [PubMed]
10. A. Filler, “Magnetic resonance neurography and diffusion tensor imaging: Origins, history, and clinical impact of the first 50 000 cases with an assessment of efficacy and utility in a prospective 5000-patient study group,” Neurosurg. 65, A29–A43 (2009).
11. J. Vogel, M. Wieczorek, C. Jud, F. Schaff, F. Pfeiffer, and T. Lasser, “X-ray tensor tomography reconstruction,” Proc. Fully3D (2015).
12. M. Bech, O. Bunk, T. Donath, R. Feidenhans’l, C. David, and F. Pfeiffer, “Quantitative X-ray dark-field computed tomography,” Phys. Med. Biol. 55, 5529–5539 (2010). [CrossRef] [PubMed]
13. V. Revol, C. Kottler, R. Kaufmann, A. Neels, and A. Dommann, “Orientation-selective X-ray dark field imaging of ordered systems,” J. Appl. Phys. 112, 114903 (2012). [CrossRef]
14. A. Malecki, G. Potdevin, T. Biernath, E. Eggl, E. Grande Garcia, T. Baum, P. B. Noël, J. S. Bauer, and F. Pfeiffer, “Coherent superposition in grating-based directional dark-field imaging,” PLoS ONE 8, e61268 (2013). [CrossRef] [PubMed]
15. A. Fehringer, T. Lasser, I. Zanette, P. B. Noël, and F. Pfeiffer, “A versatile tomographic forward- and backprojection approach on multi-GPUs,” Proc. SPIE 9034, 90344F (2014).
16. A. H. Andersen and A. C. Kak, “Simultaneous algebraic reconstruction technique (SART): A superior implementation of the ART algorithm,” Ultrasonic Imaging 6, 81–94 (1984). [CrossRef] [PubMed]
17. H. H. B. Sørensen and P. C. Hansen, “Multicore performance of block algebraic iterative reconstruction methods,” SIAM J. Sci. Comp. 36, C524–C546 (2014). [CrossRef]
18. M. R. Hestenes and E. Stiefel, “Methods of conjugate gradients for solving linear systems,” J. Res. Nat. Bur. Stand. 49, 409–436 (1952). [CrossRef]
19. J. R. Shewchuk, “An introduction to the conjugate gradient method without the agonizing pain,” (1994).
20. Q. Li and J. G. Griffiths, “Least squares ellipsoid specific fitting,” Proc. Geom. Model. Process.335–340 (2004).
21. K. Pearson, “On lines and planes of closes fit to systems of points in space,” Philos. Mag. 2, 559–572 (1901). [CrossRef]
22. H. Hotelling, “Analysis of a complex of statistical variables into principal components,” J. Educ. Psychol. 24, 417–441 (1933). [CrossRef]
23. C. M. Bishop, Pattern Recognition and Machine Learning (Springer, 2006), chap. 12.1.
24. C. Runge, “Über die numerische Auflösung von Differentialgleichungen [About numerically solving differential equations],” Math. Ann. 46, 167–178 (1895). [CrossRef]
25. M. W. Kutta, “Beitrag zur näherungsweisen Integration totaler Differentialgleichungen [Contribution to the approximate integration of total differential equations],” Z. Math. Phys. 46, 435–453 (1901).
26. C. H. Edwards and D. E. Penney, Differential Equations (Pearson, 2007), chap. 4.3.
27. C. R. Tench, P. S. Morgan, M. Wilson, and L. D. Blumhardt, “White matter mapping using diffusion tensor MRI,” Magn. Reson. Med. 47, 967–972 (2002). [CrossRef] [PubMed]
28. M. Wieczorek, J. Vogel, and T. Lasser, “CampRecon — a software framework for linear inverse problems,” Tech. Rep. TUM-I1401, TU München (2014).
29. R. L. Siddon, “Fast calculation of the exact radiological path for a three-dimensional CT array,” Med. Phys. 12, 252–255 (1985). [CrossRef] [PubMed]
