Abstract
Generalized diffusion tensor imaging (GDTI) was developed to model complex apparent diffusivity coefficient (ADC) using higher-order tensors (HOTs) and to overcome the inherent single-peak shortcoming of DTI. However, the geometry of a complex ADC profile does not correspond to the underlying structure of fibers. This tissue geometry can be inferred from the shape of the ensemble average propagator (EAP). Though interesting methods for estimating a positive ADC using 4th-order diffusion tensors were developed, GDTI in general was overtaken by other approaches, for example, the orientation distribution function (ODF), since it is considerably difficult to recuperate the EAP from a HOT model of the ADC in GDTI. In this paper, we present a novel closed-form approximation of the EAP using Hermite polynomials from a modified HOT model of the original GDTI-ADC. Since the solution is analytical, it is fast, differentiable, and the approximation converges well to the true EAP. This method also makes the effort of computing a positive ADC worthwhile, since now both the ADC and the EAP can be used and have closed forms. We demonstrate our approach with 4th-order tensors on synthetic data and in vivo human data.
1. Introduction
Generalized diffusion tensor imaging (GDTI) [1–3], was proposed to model the apparent diffusion coefficient (ADC) recovered by diffusion MRI (dMRI) when imaging the diffusion of water molecules in heterogeneous media like the cerebral white matter. Essentially, GDTI uses higher-order Cartesian tensors (HOTs) to model the spherical profile of the ADC. However, although the complex shape of the ADC reflects the complex geometry of the underlying tissue, it is well known that the geometry of the ADC does not correspond to the underlying fiber directions [4]. This can be understood from the -space formalism, where it can be seen that the ADC and the diffusion signal are in the Fourier domain of the diffusion ensemble average propagator (EAP), which describes the probability of the diffusing particles. The geometry of the EAP is a direct indicator of the microstructure of the underlying tissue or fiber bundles.
But GDTI was proposed because it overcomes the limitation of diffusion tensor imaging (DTI) [5], which is inadequate for modelling the signal from regions with multiple fiber configurations. In such regions, the HOT that is used in GDTI can model the signal and the ADC with greater accuracy than the 2nd-order diffusion tensor. Therefore, the GDTI model has been of considerable interest and has seen various developments. In particular, a number of contributions were made to estimate 4th-order HOTs in GDTI under the constraint of a positive diffusion profile since negative diffusion is nonphysical. A ternary quartic parameterization was used in [6–8], while a Riemannian approach was proposed in [9]. Other sophisticated methods were also proposed recently to estimate arbitrary even order HOTs in GDTI with the positivity constraint. In [10] the authors relied on a parameterization based on tensor decomposition into a sum of squares, and in [11], the authors used conic programming approaches to achieve this.
The GDTI model was also used to develop “biomarkers” or scalar indices such as the generalized anisotropy (GA) and the scaled entropy (SE) from HOTs modelling the ADC [3]. Additional scalar measures—in the form of invariants—were also proposed for 4th-order Cartesian tensors in [12, 13]. Overall, the GDTI model in particular and the Cartesian tensors in general have evoked great interest and have been extensively explored in dMRI. Various tensor-based models other than GDTI have been advanced and many methods (including GDTI) using tensors have been proposed to leverage the Cartesian and the algebraic structure of HOTs. A comprehensive review can be found in [14]. These indicate the importance and usefulness of HOTs in dMRI.
However, in spite of the interest in HOTs, to describe complex shaped ADCs, the tissue microstructure can only be inferred from the shape of the EAP. But computing the EAP from the HOT model of the ADC in GDTI is not an easy task. In [2], the authors proposed a numerical fast Fourier transform scheme—to emulate diffusion spectrum imaging (DSI) [15] from GDTI—to estimate the EAP and to recover the underlying fiber directions. However, this method is computationally expensive, and although the numerical Fourier transform can compute the values of the EAP at desired points, it cannot compute a continuous and differentiable function which has great advantages. That is perhaps the reason why the GDTI approach has been overshadowed by other methods that estimate the EAP or its characteristics directly from the signal, such as orientation distribution function (ODF) from Q-ball imaging (QBI), persistent angular structure (PAS-MRI), diffusion orientation transform (DOT), and spherical deconvolution (SD) [16–20]. It is interesting to note that DOT also uses the GDTI model to analytically compute the EAP on fixed shells from the ADC. However, in DOT, the spherical harmonic (SH) basis is used to model the ADC instead of Cartesian tensors.
In this paper, we propose a modification to the original GDTI model under the -space formalism, which allows us to compute a closed-form approximation of the EAP using Hermite polynomials. In this modified model, we still estimate HOTs from the signal, but these HOTs are used to describe the ADC over the entire -space instead of just its spherical profile. We show that the approximated EAP from the modified GDTI model converges well to the true EAP, and that it allows us to recover the fiber directions of the microstructure correctly. Furthermore, in our modified GDTI model, when we use 4th-order tensors, we are still able to apply the constraint of a positive diffusion profile while estimating the HOT from the signal before computing the EAP approximation from the tensor. Finally since the solution is analytical, it is fast and can be implemented efficiently.
We first test this approach on a synthetically generated dataset that simulates crossing fibers. We compare the computation time of this method with a numerical discrete Fourier transform scheme to recover the EAP from the original GDTI model. We show that with our modified approach, we are able to recover the underlying fiber layout correctly and also gain considerably in computation time. This is of great relevance in visualization and in post-processing such as tractography. We also conduct experiments on in vivo human cerebral data to illustrate the applicability of this approach on real data. We find that the peaks of the approximate EAP clearly reveal major fiber bundles and also discern crossings in the white matter.
The rest of the paper is structured as follows. In Section 2, we present the main theory and mathematical formulation of the modified GDTI model and algorithm to estimate the EAP from the GDTI-HOT. We illustrate on 4th-order tensors. In Section 3, we describe the synthetic and in vivo data, describe the experiments, and present the results. In Section 4, we discuss some properties of the approximation and how the approximate EAP behaves with respect to the true EAP. We finally conclude in Section 5.
2. Materials and Methods
2.1. Modified GDTI
We recall that the signal , in GDTI [1], is modelled using a th-order tensor (with being even), which describes the spherical or angular profile of the ADC: where when is the diffusion encoding gradient vector, , and are the components of the unit gradient vector . The second equality is a reinterpretation of the first by a rearrangement of the indices that highlights the polynomial interpretation of tensors [6]. gives the 4th-order diffusion tensor model. We also recall that in the -space formalism, the diffusion signal and the EAP are related by the Fourier transform [21]: The -space formalism entails the “pulsed” gradient condition that , implying that .
For , is the DTI model, whose Fourier transform is well known to also be a Gaussian, which corresponds to the free diffusion propagator. We denote by the EAP computed from the Fourier transform of . However, for general , closed forms for the Cartesian Fourier transform of are hard to compute, since in Cartesian coordinates, is not separable in , the components of . In [2], where a method for recovering the EAP from GDTI is proposed, is computed numerically by evaluating more or less densely in -space and by computing its fast Fourier transform.
In this section, we propose to modify the original GDTI model by making (1) separable in Cartesian coordinates. This is done by realizing that GDTI in fact uses two orders and for the radial and the angular components, respectively, where in GDTI and . In this formulation, is a constant with units that makes the exponent unit free when and . The first equality is written in the components of the unit gradient vector , while the second equality is in the components of the reciprocal space vector . To pass from the first to the second equality, the components of have to be multiplied by , the norm of , raised to the appropriate power, . To leave the equality unchanged, it is therefore necessary to also multiply by , which results in the second equality. This reformulation of GDTI allows to become separable in when .
Although the original formulation of GDTI uses a Cartesian HOT, it was in essence written in spherical coordinates since the HOT was evaluated only along the unit gradient vector . Alternately, the spherical coordinates also become evident from the two separate orders and for the radial and the angular parts. By equating the two orders , our modification converts the signal formulation to Cartesian coordinates and recouples the radial and the angular parts. The HOT is now evaluated over the entire -space. This reformulation allows us to compute an analytical Fourier transform of in Cartesian coordinates.
Interestingly, in spite of this reformulation, the signal in (4) still retains a monoexponential form parameterized by the diffusion HOT , like in the original formulation. In (1), the negative logarithm of the signal is , where . In the modified model, if we denote , then This makes it evident that we can again estimate the HOT from the diffusion signal in such a way that its diffusion profile is positive, that is, . Here, the diffusion profile is no longer a function on the sphere but rather a function on the entire -space. However, from this equation, we also see that when , the methods that were developed to estimate a 4th-order tensor with a positive diffusion profile for the GDTI model in (1), can all be directly applied to the modified model in (4). Therefore, using the modified model in (4), it is possible to estimate a 4th-order HOT from the signal which satisfies a positive diffusion profile, before computing the EAP from this HOT. In this paper, we use the ternary quartic parameterization proposed in [8].
2.2. Fast and Analytical EAP Approximation
Our solution for the EAP from the modified GDTI model pivots around the following property of the Fourier transform: where stands for the Fourier transform. If we employ for , then its Fourier transform is . However, the derivatives of the Gaussian function generate the Hermite polynomials . Therefore, The generalization to 3D is simple since the Gaussian function is separable in the variables.
To take advantage of this property of the Fourier transform, and those of the Gaussian function, for computing a closed-form approximation of the EAP from , that is, its Fourier transform , we propose to expand as a multivariate polynomial multiplied by a 3D Gaussian function: where is a constant with units to render the exponent unit free, and the new coefficients in the polynomial expansion contain the imaging parameter and the coefficients of the HOT . If this expansion was possible, then would become separable in .
Such an expansion can be achieved from a few manipulations and a Taylor expansion: where the summation in the first equality is over such that , as denoted in (4), and . Since is an exponential function , we define as the th-order Taylor expansion of in the variables . Therefore, is a trivariate polynomial of degree plus an error term of degree . Ignoring the error term, has the required form . Therefore, we can define the th-order approximation of the signal: Since is the Taylor’s expansion of an exponential function converges to uniformly over all as is made large. Therefore, converges to uniformly over as is made large.
As is separable in , it is possible to compute a closed form for its Cartesian Fourier transform, which is also separable. Using the property in (7): For large , the approximation converges to the true EAP . In practice, we use .
We thus find a closed-form approximation of the EAP from the modified GDTI model of the ADC using HOTs. The solution is a polynomial multiplied by a Gaussian. Therefore, the polynomial can be interpreted as the correction to the free diffusion Gaussian EAP due to the complex heterogeneous medium.
An alternate interpretation to this method can be found from (10) and (11), which avoids modifying the GDTI model. Equations (10) and (11) resemble closely the formulation of the signal in the expansion of the cumulant generating function, and the approximation of the EAP using the Gram-Charlier series, as proposed in [22]. While in [22], in the cumulant expansion, the signal is expanded in the standard polynomial basis with the cumulants as the coefficients, in (10) the signal is in fact expanded in a subset of the standard polynomial basis. Since the Fourier transform of a monomial multiplied by a Gaussian is a Hermite polynomial multiplied by a Gaussian, in [22] too, the EAP is approximated in the Hermite polynomial basis—again with the cumulants as the coefficients using the Gram-Charlier series. Likewise in (11), the EAP is approximated in the Hermite polynomial basis.
The difference between this method and [22] lies in the fact that while [22] uses the entire polynomial (Hermite polynomial) basis to expand the signal (EAP), this method uses only a subset of these bases. Therefore, the coefficients are no longer the cumulants. Or, in other words, if the entire polynomial basis had been used here, then would have become the cumulants. Also, the coefficients are not estimated directly from the signal, though they can be if (10) were used, but are instead computed from the coefficients of the tensor and the Taylor expansion in (10). Therefore, changing the order of the approximation has an effect on the approximated EAP, since it adds or subtracts terms in (11). However, as shown in the experiments, this does not affect the direction of the peaks of the approximate EAP.
We program an efficient implementation of the proposed method through symbolic computation. Using Maple and assuming (4), we expand into a Taylor series in the variables up to predefined orders . This expansion automatically computes for us the new coefficients from the coefficients of the HOT (10). The EAP approximation , is then generated by again computing the Fourier transform of symbolically. The expansion of the EAP is then converted to C-code using Maple, which is compiled. This routine therefore takes the imaging parameters as input, namely, and the coefficients of that are estimated from the diffusion signal. are taken to be equal to 1.
3. Experiments and Results
Although we developed the theory for arbitrary , for the following experiments we consider , that is, and . This is because, as we have seen, for , we can employ the estimation techniques that guarantee that the 4th-order HOT has a positive diffusion profile, . In all the following the 4th-order HOT is estimated from (3), with using the method described in [8]. The estimation in [8] is described for (1), which depends on the -value, that is, . We adapt this to (3) by replacing the -value by the imaging parameters and the ADC by the . We test the approach first on synthetic data and then on in vivo human cerebral data.
3.1. Synthetic Dataset
To conduct controlled experiments with known ground truths, we use a multitensor approach to generate synthetic DWIs [23]. The EAP corresponding to a single fiber is taken to be an anisotropic free diffusion Gaussian distribution, parameterized by a covariance tensor in its canonical coordinates. is rotated using rotation matrices to orient the fiber in space. We generate the signal DWIs for the single fiber by considering the -space formalism and taking the Fourier transform of the Gaussian EAP, which results in the anisotropic Stejskal-Tanner signal equation. Multiple crossing fibers are simulated by considering an EAP, that is, the weighted sum of free diffusion Gaussians, where each Gaussian represents a fiber oriented in space. The signal DWIs for a multifiber or crossing fiber is derived easily in the same fashion as , such that , with a rotation matrix, represents the DWI along the th gradient direction, and are the number of fibers crossing in the voxel. We use a -value of to generate the signal and corrupt it with a Rician noise with signal to noise ratio (SNR) of 30. The gradient directions are considered isotropically spread out on the sphere along 81 encoding directions. Since the dataset is generated from a fixed -value, we consider the imaging parameter , which allows us to compute .
3.2. In Vivo Human Cerebral Dataset
The in vivo human cerebral dataset, described in [24], was acquired with a whole-body 3T Siemens Trio scanner, with an 8-channel array head coil and maximum gradient strength of 40 mT/m. The DWIs were acquired using spin-echo echo planar imaging (EPI) (time repetition , echo time , image matrix, , 72 slices with 1.7 mm thickness (no gap) covering the whole brain). The diffusion weighting was isotropically distributed along 60 encoding directions, with a -value of . Seven images without any diffusion weightings were placed at the beginning of the sequence and after each block of ten DWIs as anatomical reference for offline motion correction. Random variations in the data were reduced by averaging 3 acquisitions, resulting in an acquisition time of about 45 minutes. The SNR in the white matter of the image was estimated to be approximately 37. The motion correction for the DWIs was combined with a global registration to T1 anatomy images. The gradient direction for each volume was corrected using rotation parameters. The registered images were interpolated to the new reference frame with an isotropic voxel resolution of 1.72 mm.
3.3. Synthetic Dataset Experiment
A first example of the synthetic data and the results of our method are shown in Figure 1. In this proof of concept experiment, we consider two fibers crossing perpendicularly in a voxel with equal weights. On the left, in Figure 1(a), is shown the positive modelled by a 4th-order tensor which was estimated using [8]. In Figure 1(b) are shown the profiles of the analytically approximated EAP () for increasing norms of the vector from 12 m to 20 m. As expected, we observe the desired “sharpening” effect from the different profiles of the EAP as the probability of water molecules diffusing declines sharply along nonfiber directions as the radius of the probability distribution is increased. This provides a strong motivation for estimating positive higher-order diffusion tensors from the (modified) GDTI model, since now from the analytically approximated EAP it is also possible to infer underlying fiber directions.
(a)
(b)
In the main synthetic data experiment, we consider two fiber bundles crossing or overlapping in a way that makes them converge and diverge. This changes their crossing angles gradually in the region where they intersect. The voxels outside the fiber bundles are generated using an isotropic diffusion profile. We set three goals for this experiment. First, we test if our analytically approximate EAP can recover the three types of voxel models from the noisy DWIs, namely, isotropic, single fiber, and crossing fiber voxels. Second, we compare the computation time of our proposed method to the numerical Fourier transform approach. Finally, third, we also conduct tests on the effects of the estimation order on the EAP; these are discussed in Section 4.
The layout of the synthetic dataset fibers and the result of the estimated from 4th-order HOTs and the analytical EAP approximations of order 7, are presented in Figure 2. It shows the two simulated fiber bundles and the three types of voxel models that constitute the synthetic dataset experiment. The effect of the EAP transformation of the are highlighted in the zooms. Although the profiles clearly indicate the regions with complex microstructures, that is, fiber crossings, the geometry of the is not aligned with the fiber bundle directions. The peaks of the EAP on the other hand correctly indicate the underlying fiber directions.
To evaluate the validity of the EAP approximation, we have proposed the angular profiles of the EAP, for fixed that are presented in Figure 3. In the two zooms, we take a closer look at some of the voxels in the crossing regions. In the top zoom, we see that the peaks of the angular profile of the EAP correctly detect the changing angle between the converging or diverging fiber bundles. In the bottom zoom, we see the three different types of voxels recovered by the approximate EAP, namely, the isotropic, the single fiber, and the two fibers crossings. Although the isotropic voxels also have some peaks, the peaks of the EAPs representing crossings are much more sharp, and it is easy to distinguish these two types of voxel geometries. The peaks in the isotropic voxels are caused by the signal noise and augmented order (4th) of the HOT. Nonetheless, as opposed to many other higher-order models such as PASMRI or SD, the voxels not containing fiber bundles are clearly identifiable as isotropic.
Speed is of great relevance in visualization and in processing after local estimation, such as in tractography. The closed-form of makes it computationally efficient, especially since the expression for a fixed can be hard coded and compiled. In the synthetic data experiment, we compare this approach to a numerical Fourier transform of the GDTI model. For visualization and comparison, we consider the whole slice, which is partially seen in Figure 3, with voxels. For the implementation of the numerical Fourier transform, we evaluate the GDTI model (1) on a Cartesian grid. We evaluate the numerically computed EAP on a coarse spherical mesh with 162 vertices. The results are presented in Table 1. The computation time on our computer was 526 s. We then compute , but this time on a finer spherical mesh with 2562 vertices. The computation time on the same computer was 73 sec. Despite the finer mesh, is about seven times faster than the regular discrete Fourier transform. On the coarse mesh with 162 vertices, the computation time for was about 10 s.
3.4. In Vivo Dataset Experiment
For the in vivo human cerebral dataset described above, we make certain assumptions about the imaging parameters since this dataset was acquired using a twice refocused Reese sequence [25], and not a standard pulsed-gradient spin-echo (PGSE) sequence. The gradient durations used in the Reese sequence were , , , and . As suggested in [26], Reese sequence parameters are sometimes adapted to the standard PGSE parameters with , and as the time between the start of and the start of . However, since the application times of were unknown, we assume that , which implies .
We choose a coronal slice from the in vivo dataset where three fiber bundles are known to cross. The resulting analytically approximated EAPs from 4th-order tensors in the coronal slice are shown in Figure 4. In the plane horizontally and diagonally is the corpus callosum (CC), top to bottom is the corticospinal tract (CST), and going through the plane is the superior longitudinal fasciculus (SLF). The 4th-order HOTs were approximated from this dataset. In Figure 4, are shown the estimated order 7 approximations of the EAP. The zooms highlight the crossings between the major fiber bundles. In the main zoom is the region where the three fibers, the CC, the CST, and the SLF, intersect each other. In the upper secondary zoom, the crossing between the CC and the cingulum is highlighted, which occurs due to partial voluming. In the lower secondary zoom, are seen the main voxels with three peaks, which correspond to the crossing between the three fiber families—the CC, the CST, and the SLF.
4. Discussion
From the synthetic dataset and in vivo dataset experiments above, we were able to show that it is possible to recover the microstructure or fiber directions of the underlying tissue from our proposed analytical approximation of the EAP that was computed from 4th-order tensors. However, it is important to realize that the proposed approach only approximates the true EAP up to the truncation order of the Taylor series in (10). In this section, we validate the effects of this truncation or approximation on the fiber directions estimated by the approximate EAP. To do this, we again consider the synthetic dataset and look at a region where all three kinds of voxel models are present. We then compute the analytical EAP for approximation orders of . We compare the angular profiles of these different approximations at a fixed radius to the angular profiles of the fixed approximation at for varying radii.
Figure 5 shows the effect of the Taylor expansion order on the EAP approximation. The six images are zooms into a region where the two fiber bundles converge and cross. In the top row, we present , with evaluated for the probability radius . Increasing adds more terms to the EAP approximation in (11), which adds more corrections to the approximation, making it converge better to the true EAP. As the approximation , is corrected, it shows sharper peaks and narrower crossings for , than , for the same probability radius. However, this also increases the computation time. But the peaks of the lower-order approximations seem to be well aligned with the higher-order approximations. In other words, the peaks maintain their angular alignment, although they lose sharpness, and the EAP loses angular resolution, and narrow crossings become harder to discern. However, the angular resolution can be recovered, and the peaks “sharpened” in the lower-order approximations by increasing the probability radius, which saves computing time. This is shown in the bottom row, where we show for the probability radius varying from – These experiments reveal that the effect of the Taylor expansion order is to underestimate the EAP in the approximation. Therefore, we use the order 7 approximation , as a good trade-off between convergence to the true EAP and computation time.
5. Conclusion
GDTI was developed to model complex ADC profiles which was an inherent shortcoming of DTI. GDTI uses HOTs of order to model a complex ADC geometry. However, the shape of the ADC does not correspond to the underlying fiber directions. The microstructure of the tissue can be inferred from the geometry of the EAP, where in the -space formalism, the EAP and the diffusion signal are related by the Fourier transform. But it is not easy to compute the EAP, , from the HOT model of the signal in GDTI.
We overcome this hurdle by modifying the ADC model of GDTI, which allows us to approximate by a multivariate polynomial approximation and by proposing a novel closed-form approximation of using Hermite polynomials. The solution is a polynomial times a Gaussian; therefore, the polynomial can be interpreted as the correction to the Gaussian EAP due to the inhomogeneous medium. An alternate explanation can be used to explain this method, where the signal is expanded in the polynomial basis, and the EAP is expressed in the Hermite polynomial basis, which establishes the similarity of the proposed method to [22]. Also, since the solution is analytical, it is fast, and the approximation converges well to the true EAP.
In case of an order 4 HOT, this method can be directly adapted to the methods proposed for estimating 4th-order diffusion tensors with positive diffusion profiles. Therefore, it is possible to estimate a 4th-order HOT with a positive diffusion profile using this modified model before approximating the EAP. The experiments show that estimating only the 15 coefficients of a 4th-order HOT are enough to reveal the underlying fiber bundle layout. However, this is dependent on the order of the Taylor expansion used. Although the order of the expansion does not change the angular alignment of the peaks of the approximate EAP, it does affect its angular resolution or its capability of discerning narrow crossings. Increasing of the order increases the corrections to the approximation, which improves this angular resolution. However, it also increases the computation time. The angular resolution can be recovered in lower-order approximations, by increasing the probability radius, which saves computation time. However, this overall effect indicates that the truncation in the Taylor expansion has the effect of underestimating the true EAP in the approximation.
Abbreviations
| CC: | Corpus callosum |
| CST: | Corticospinal tract |
| DOT: | Diffusion orientation transform |
| DSI: | Diffusion spectrum imaging |
| DTI: | Diffusion tensor imaging |
| DWI: | Diffusion-weighted image |
| EAP: | Ensemble average propagator |
| GA: | Generalized anisotropy |
| GDTI: | Generalized DTI |
| HOT: | Higher-(than 2) order (Cartesian) tensor |
| QBI: | Q-ball imaging |
| ODF: | Orientation distribution function |
| PASMRI: | Persistent angular structure MRI |
| PGSE: | Pulsed-gradient spin echo |
| SD: | Spherical deconvolution |
| SE: | Scaled entropy |
| SH: | Spherical harmonic (basis) |
| SLF: | Superior longitudinal fasciculus. |
Acknowledgments
The authors would like to thank Dr. A. Anwander from the Max Planck Institute for Human Cognitive and Brain Sciences, Leipzig, Germany, for providing them with the in vivo human dataset. A. Ghosh was partially supported by the ANR Project NucleiPark and the France-Parkinson Association.