Abstract
Diffusion magnetic resonance imaging (dMRI) is the only technique to probe in vivo and noninvasively the fiber structure of human brain white matter. Detecting the crossing of neuronal fibers remains an exciting challenge with an important impact in tractography. In this work, we tackle this challenging problem and propose an original and efficient technique to extract all crossing fibers from diffusion signals. To this end, we start by estimating, from the dMRI signal, the so-called Cartesian tensor fiber orientation distribution (CT-FOD) function, whose maxima correspond exactly to the orientations of the fibers. The fourth order symmetric positive definite tensor that represents the CT-FOD is then analytically decomposed via the application of a new theoretical approach and this decomposition is used to accurately extract all the fibers orientations. Our proposed high order tensor decomposition based approach is minimal and allows recovering the whole crossing fibers without any a priori information on the total number of fibers. Various experiments performed on noisy synthetic data, on phantom diffusion, data and on human brain data validate our approach and clearly demonstrate that it is efficient, robust to noise and performs favorably in terms of angular resolution and accuracy when compared to some classical and state-of-the-art approaches.
1. Introduction
The diffusion magnetic resonance imaging or dMRI [1] is a magnetic resonance imaging (MRI) modality which is particularly suited to study and characterize the white matter neuronal architecture of the brain in vivo and noninvasively. The diffusion tensor model (DTI) introduced by Basser et al. in 1994 [2] was the first technique used for the fiber bundles reconstruction. However, due to the assumption of single fiber bundle per voxel, this technique is ineffective in regions where the fiber bundles intersect. Therefore, the fiber bundles reconstructed by tractography algorithms based on DTI are unreliable. To overcome the limitations of the DTI model, new high angular resolution techniques (HARDI) have been proposed, such as the diffusion spectrum imaging (DSI), the Q-ball imaging (QBI) [3, 4], or the symmetric high order tensors (SHOT) [5]. These techniques allow estimating the diffusion orientation distribution function (ODF) whose maxima are aligned with the orientations of the underlying fibers. The DSI technique provides the real ODF by measuring the diffusion signal on a whole 3D Cartesian grid in the q-space. However, this method is impractical in clinical studies because it requires an important acquisition time due to the huge number of samples and it requires too a very high gradient magnitude. The QBI model approximates the diffusion ODF function directly from the raw HARDI diffusion signal acquired from a spherical sampling of the diffusion space [3, 4]. Although the QBI-ODF contains the angular information by having its maxima aligned on the orientations of the underlying fibers, it has a low angular resolution by failing to reconstruct correctly the fibers crossing with angles less than [3, 6]. A sharper function called fiber ODF or FOD function can be calculated from the ODF by using the spherical deconvolution techniques [7]. The FOD function has also its maxima aligned on the underlying fiber orientations; the FOD allows a gain in the angular resolution up to . Traditionally the ODF and FOD functions are described in spherical harmonics (SH) basis; the angular resolution of these functions directly depends on the order of the SH basis: an acute angular resolution requires a high order of the SH basis. However, in case of higher orders, these ODF and FOD functions are prone to negative lobes due to noise. The high order symmetric tensors or equivalently the homogeneous polynomials have been proposed to model and reconstruct the FOD function [8] called CT-FOD for Cartesian tensor-FOD, whose maxima correspond exactly to the orientations of the underlying fibers. By imposing efficiently the positivity constraint, the CT-FOD function appears as an alternative to the FOD described in the SH basis. Furthermore, thanks to its polynomial form, the maxima of the CT-FOD function can be easily located. Therefore, in our work we start by estimating the CT-FOD function to reconstruct the symmetric high order tensor from the diffusion weighted-MRI (DW-MRI) data. Due to the importance of tractography and its increasing interest in clinical practice, it is important to accurately extract the fiber orientations to perform a reliable and accurate tractography. An efficient and accurate approach to perform this crucial and necessary preprocessing step is to extract the fibers orientations as the CT-FOD maxima. Different methods for extracting maxima of high order tensors exist in the literature; in [9], Bloy and Verma proposed to determine the fiber directions using the concept of Z-eigenvalues introduced by Qi in 2005 [10]. This maxima localization method, the so-called traditional method, suffers from a low angular resolution and does not allow recovering crossing fibers at angles below [11]. The high order tensor decomposition in rank-1 tensor has been proposed to the maxima extraction issue in [12], and the low-rank decomposition approximation method known as CANDCOMP/PARAFAC (CP) [13] was used. It was proved in previous works [11, 12] that the tensor decomposition approaches recover the crossing fibers with a better angular resolution than the traditional methods of maxima localization. However, the CP-decomposition approximation requires predefining the tensor rank, that is, knowing the number of fibers in a voxel, which is impossible apriority. Furthermore, PARAFAC uses the alternating least squares (ALS) algorithm, which is a nonlinear optimization algorithm whose convergence is not guaranteed and depends on the initialization. In this paper, we propose to find the orientations of the fiber bundles from diffusion signals using an analytical decomposition of symmetric high order tensor; for lightness and clarity of the paper we will use the abbreviation Adecomp-SHOT for analytical decomposition of symmetric high order tensor. This Adecomp-SHOT method was initially proposed by Brachat et al. in 2010 [14] but remains in the theoretical field. However, the Adecomp-SHOT seems a priori interesting for the fiber orientations issue because, unlike the suboptimal CP-decomposition approximation, the Adecomp-SHOT is an analytical one and not restricted to subgeneric ranks. Thus, rather than CP-decomposition, Adecomp-SHOT would provide a minimal decomposition; this aspect is particularly interesting in the fiber orientations search since it would give the whole underlying fiber orientations without any apriority. Therefore, in the following we propose an original and efficient Adecomp-SHOT based approach to extract the fiber orientations from diffusion weighted-MRI data and we prove through many validations tests the effectiveness of the proposed method. The Adecom-SHOT is based on the SHOT; therefore we propose to use the CT-FOD for reconstructing the SHOT from the DW-MRI signal, since the CT-FOD constitutes the state-of-the-art. We begin this paper by presenting first the CT-FOD algorithm, before explaining the CP-decomposition and providing a detailed version of the Adecomp-SHOT algorithm; in this section we will also present the results of the intrinsic study done on both of Adecomp-SHOT algorithm and CP-decomposition. Then we describe our Adecomp-SHOT based approach to the fiber orientations search. Finally, we finish by presenting our validations and results on synthetic, phantom, and in vivo human brain data and our conclusions.
2. Materials and Methods
2.1. Symmetric Fourth Order Tensor Coefficients from the Diffusion Data
The diffusion signal corresponding to the acquisition parameters, , , is given by the convolution of the CT-FOD function , modeled by a Cartesian and positive definite symmetric high order tensor of 3 dimensions, with a Watson function [7, 8]: , with being the diffusivity coefficient calculated from a 2nd order tensor, that we estimate from a single fiber response having a high fractional anisotropy (); is the gradient of magnetic fields, stands for the weights b-values, and represents a set of unit vectors sampling the diffusion space. Algorithm 1 describes the estimation of unique coefficients of the symmetric tensor FOD from the diffusion data.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
In the following we are interested to decompose symmetric fourth order tensors of dimension 3 ( and ).
2.2. Symmetric Tensor Decomposition
Symmetric high order tensors appear mostly as multivariate functions (more than two variables), and high order tensors decomposition allows deducing the geometric and invariance properties of a tensor. Therefore, the tensor decomposition raises interest in many practical domains, first in chemometrics [13] and psychometrics fields and then in electrical engineering and electronics [16], in particular for the antenna array processing [17], or else in telecommunication field [18]. Also, tensor decomposition appears very useful in data analysis and in the arithmetic complexity area [19]. Recently the interest in tensor decompositions has expanded to the neurosciences field; we cite among several applications the use of the symmetric tensor decomposition to the problem of extracting the fiber orientations of the white matter in dMRI. However, to date, in dMRI the decomposition problem is still solved with a low-rank approximation method known as CANDCOMP/PARAFAC (CP).
In the following two subsections, we first start by describing the classical CANDECOMP/PARAFAC (CP) method to decompose symmetric high order tensors and then present in detail the Adecomp-SHOT approach we propose to analytically decompose a symmetric tensor of any order and any dimension in a minimal sum of rank-1 terms.
2.2.1. Numerical Method: CP Low-Rank Approximation
The tensor decomposition problem consists in writing a given tensor, in sum of outer product of vectors, that is, rank-1 tensor, and that with a minimal number of terms, the number of terms corresponding to the minimal tensor rank. Considering a symmetric tensor of order and dimension , the minimal decomposition of this tensor should be in the following form: with : the rank of , “”: the outer product, and : the rank-1 tensors (vectors).
However, determining the rank of high order tensors (order ) is a hard mathematical and NP-complet problem. Therefore, a low-rank numerical approximation of the decomposition (4) has been proposed in [13], where the authors approximate the tensor by another tensor whose rank is inferior to the minimal or generic tensor rank; this numerical decomposition method is known as CANDECOMP/PARAFAC (CP): with = 1 and : the weights of the rank-1 tensors .
is the subgeneric rank of the tensor , for symmetric tensors . Thus, the CP-decomposition of the tensor of order and of an unknown minimal rank is done by a nonlinear minimization in of (5) for a given subgeneric rank ; the nonlinear minimization problem is solved using the alternating least squares algorithm (ALS):
Inverse Problem. To evaluate the intrinsic behavior of the CP method, we have simulated an inverse problem by generating fourth order symmetric tensors of rank-2 according to (3). These tensors are constructed from two crossing vectors with variable angles from to and weighted by the same weight . The purpose is the evaluation of both the ability of the method to render the correct solutions and the angular resolution in such ideal case where data perfectly satisfy the decomposition model.
The results of the intrinsic study are illustrated in Figure 1; Figure 1(a) represents the mean error between the simulated vectors or rank-1 tensors and the CP-decomposition solutions, and Figure 1(b) gives the rank-1 tensors weights found by the CP-decomposition, according to the separation angles.
(a)
(b)
From Figure 1(a) we notice that although the tensors to decompose are constructed as to satisfy the decomposition model described in (3), the CP-decomposition begins to give an incorrect decomposition when the solutions are separated with small angles (<); thus, Figure 1(b) shows that only one vector of the two simulated is detected with the weights and for crossing angles inferior to . We conclude that the accuracy and the angular resolution of the CP-decomposition are intrinsically limited. Furthermore, the CP-decomposition has another important limit which is the requirement to predefine the rank of the decomposition.
The use of the CP-decomposition algorithm in dMRI to detect the fiber orientations was proposed in 2011 by Jiao et al. [12]. The authors considered an approximation of the decomposition with a low-rank value (rank ) in order to extract two crossing fiber orientations. The extracted fiber orientations correspond to the rank-1 vectors obtained from the decomposition. Although tensors decomposition is more efficient in terms of angular resolution and accuracy than the traditional maxima localization methods [11, 12], the CP-decomposition could not guarantee the recovering of the whole underlying fiber orientations since the minimal rank of the tensor is not a priori known. Moreover, the convergence of the ALS algorithm is not guaranteed and depends on the initialization.
In 2010, Brachat et al. proposed [14] to solve the symmetric high order tensor decomposition problem analytically; in the remainder of the paper we denote this method: Adecomp-SHOT. Rather than the CP numerical approach, the Adecomp-SHOT method gives a minimal decomposition without any apriority. However, to date the Adecomp-SHOT remains theoretical and not yet expanded to the practical issues or evaluated on physical phenomenons. Due to its ability to render a minimal decomposition, we naturally expect that the Adecomp-SHOT applied to the fiber orientations search in dMRI would be more interesting than the CP low-rank approximation method.
2.2.2. Analytical Method: Adecomp-SHOT
The Adecomp-SHOT method initially proposed by the authors Brachat et al. [14] is a generalization of the Sylvesters theorem [14], initially introduced for binary cases and extended to larger dimensions. Thus, the Adecomp-SHOT algorithm is able to decompose a symmetric tensor of any order and any dimension, in a minimal sum of rank-1 terms.
Consider a symmetric tensor of order and dimension given in the following polynomial form: with a homogenous polynomial of order in variables, and the coefficients of the homogenous polynomial .
An affine decomposition of consists in writing as a sum of powers of rank-1 linear forms [14] as follows: with representing scalars weights, rank-1 linear forms in with the coefficients , and the minimal rank of , that is, of the tensor.
is a rank-1 linear form of dimensions [14]: with for .
An affine decomposition of exists, if and only if an affine decomposition of exists; thus, the decomposition of is equivalent to decomposing ; with being the linear form associated with in the dual space , the coefficients of are then calculated from the coefficients of [14] as follows: is obtained by unhomogenizing according to the variable ; this is done by dividing each monomial of the homogenous polynomial by an appropriate power of . Thus, is a nonhomogenous polynomial of degree .
The necessary and sufficient existence conditions of the decomposition are based on the rank conditions of the Hankel matrix and the commutation properties. The Hankel matrix is a matrix with corresponds to the number of unique coefficients of a -order symmetric tensor of dimensions. The elements of the Hankel matrix are computed from the coefficients of ; we note that the elements corresponding to monomials with total degree higher than the polynomial degree are unknown. We provide an appendix (Appendix section) where we give an example of a step by step decomposition of a 3 dimensional 4th order tensor of rank-4, illustrating how the Adecomp-SHOT algorithm works and particularly we show in Step 2 of the Appendix section how the Hankel matrix is constructed from a homogenous polynomial.
To obtain the points of (8), we calculate from the Hankel matrix the multiplication matrix as described in Step 5 of Algorithm 2, with , for a given rank ; then we resolve the generalized eigenvalue problem and we deduce the weights by simply resolving a linear system. The readers can refer to the example given in the Appendix section for much more details. The critical part of Algorithm 2 is the step of extending the Hankel matrix, Step 6 of Algorithm 2, to verify the stability of the rank in case of higher ranks. Indeed, when the Hankel matrix is not totally defined, the extension requires finding the unknown parameters of the Hankel matrix satisfying the commutation properties of the multiplication matrix with ; this leads to solving a nonlinear equations system. Therefore, for higher ranks the uniqueness of the decomposition is not guaranteed.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Inverse Problem. Once again, we have simulated an inverse problem, this time for the Adecomp-SHOT algorithm; thus, rank-2 fourth order tensors were constructed from a sum of two rank-1 linear forms of order four, as described in (7). We have simulated the same configuration as in Section 2.2.1 by constructing the fourth order tensors from two crossing rank-1 linear forms at varying angles from to and weighted by the same weight . Our results of the intrinsic study of the analytical decomposition method are thus given in Figure 2.
(a)
(b)
Contrarily to the CP-decomposition, the analytical approach clearly renders the correct decomposition whatever the separation angles. As shown in Figure 2(a) the obtained mean error is zero, and the two rank-1 tensors are recovered with the correct weights as represented in Figure 2(b). Furthermore, the analytical decomposition is minimal; that is, the rank of the tensor is automatically found without any assumption such that it is required for the CP-decomposition. To confirm the ability of the method to give always a minimal decomposition regardless of the rank of the tensor, further tests on higher rank tensor have been conducted; Figure 3 shows the results of decomposing a rank-3 symmetric fourth order tensor constructed from 3 crossing rank-1 tensors, according to (7), at angles decreasing from to . Figure 3 shows that the Adecomp-SHOT method gives once again a correct decomposition with a zero mean error whatever the separation angles between the 3 origin rank-1 tensors. Other tests on symmetric fourth order tensors of rank-4 and rank-5 have been conducted; the results of these experiences show that the Adecomp-SHOT method still gives the correct decomposition in case of rank-4 while for a rank-5 fourth order tensors the decomposition is found with insignificant angular error not exceeding . However, by increasing the order of the tensor from 4 to 6 the error drops to zero.
(a)
(b)
2.3. Fiber Directions from Diffusion Data Using the Analytical Decomposition of Fourth Order Tensor
In this section we propose to extract the fibers orientations from the dMRI signal by decomposing analytically the three dimensional fourth order CT-FOD using the Adecomp-SHOT. The coefficients of the fourth order FOD are estimated from the dMRI data as described in Section 2.1. Thus, we are interested to decompose the Cartesian FOD tensor in sum of powers of rank-1 linear forms: with minimal representing the tensor rank; are rank-1 linear forms in 3 variables with the real normalized coefficients , and are the real positive weights.
The normalized coefficients of represent the Cartesian coordinates of the fibers orientations, weighted by the scalar , and represents the number of crossing fibers bundles in each voxel.
However, as described in Section 2.2.2 the affine decomposition of a homogenous polynomial of any order and any dimension is done assuming that all coefficients of the first Cartesian coordinate in the decomposition are nonzero; that is, , for , with the tensor rank. This constraint implies that only fiber orientations whose first coordinate coefficient is nonzero can be detected. To avoid missing fiber orientations, we propose to introduce a coordinate changing in case where maxima are located in undetectable area. Also, not doing this coordinate transformation systematically and imposing a stopping criterion, we propose to make a first exhaustive search on the FOD by discretizing it on unit sphere and then localize roughly its maxima. A given FOD function , with the gradient of the magnetic field, can be discretized on units sphere as follows: with the number of samples; is a rank-1 fourth order tensor constructed in the orientation and represents the value of the FOD in this orientation, and the 3 dimensional unit vectors are obtained by a uniform tessellation of unit sphere. Thus, the are found by solving the linear system , with representing a vector of weights , vectors of 15 unique coefficients of the fourth order tensor FOD and a matrix; contains the 15 unique coefficients of the rank-1 fourth order tensor constructed from unit vectors . Once the FOD is discretized, we check the first coordinate of the FOD values weighted by corresponding to the FOD lobes, and we make a coordinate changing before decomposing the tensor if these coordinates are close to zero. Obviously, after decomposition, the inverse coordinate transformation is required to bring back the resulted rank-1 tensor to the origin coordinate system. In order to preserve the angles between two vectors and their lengths, the transformation matrix should be orthogonal. Therefore, to preserve the fiber orientations we use a rotation transformation; considering the initial coordinates or variables , is obtained from by the following linear transformation: As mentioned, the matrix is an orthogonal matrix or a rotation matrix. This coordinate transformation can be done either on the homogenous polynomial as a change of variable using a rotation transformation or simply by rotating the fourth order tensor coefficients. This coordinate transformation insures recovering the entire crossing fibers whatever its locations in the Cartesian space.
The Adecomp-SHOT as initially described in [14] would provide a minimal decomposition with a minimal analytical rank, without any constraints on the weights or values; in theory this is not a problem, but when we are interested in detecting the fibers orientations, negatives values for the weights or complexes coefficients do not correspond to any physical meaning. To overcome this limit, we propose ignoring the with complexes coefficients found at Step 7 of Algorithm 2 and then solve the linear system (14) with only the linear forms with real coefficients: represents the tensor rank and . This linear system can be written as a matrix equation , with a vector of length containing the fibers weights and is a matrix containing the polynomials coefficients of the rank-1 linear forms of order 4 and is a vector of length 15 containing the coefficients of the fourth order fiber orientation distribution function . To impose the positivity constraint on the weights values, we propose to solve the minimization problem 2.17 using the well-known Lawson and Hansons NNLS algorithm [15]: However, to verify that a set of such that exists we check the norm ; if the residual norm is less than 1, we assert that it is sufficient to consider the resulted decomposition accurately and the weighted by correspond to the maxima of ; else, we propose to relaunch the decomposition by doing a change of coordinates. This trick allows imposing the positivity constraint and increasing the accuracy of the decomposition.
Finally, to take into account the effect of the noise due to the diffusion model, we have introduced a heuristic cleaning that consists in removing all the fiber orientations weighted by and merging fibers separated by angle [20].
3. Results and Discussion
To validate our proposed crossing fibers detection method, we conduct many tests, first on synthetic diffusion dataset simulated with a multitensor model with 60 gradient directions and a -value of 3000 s/mm2; these data represent crossing fibers with variable separating angles from to . On these synthetic diffusion data we have compared our method to other methods in literature such as CP-decomposition based method and the Z-eigenvalues based approach; the results of this comparison are illustrated in Figure 4. Then, the synthetic dataset is corrupted with a rician noise of different signal to noise ratio (SNR = 40, 30, 20, and 10); 100 trials of noise are performed for each SNR level and for each separating angle; the results of the effect of a rician noise are presented in Figures 5 and 6 and Table 1. Finally, our method is tested on phantom and on in vivo human brain diffusion data as illustrated in Figures 7 and 8.
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(a)
(b)
(c)
(d)
3.1. Validations on Synthetic Diffusion Dataset
Comparison with CP-Decomposition and Z-Eigenvalues Based Approaches. In order to compare the angular resolution of our fiber extraction approach to the one of the CP-decomposition and the Z-eigenvalues based methods, we have conducted experiments on same noise-free synthetic diffusion dataset. From these data we reconstruct the fourth order CT-FODs and then we extract the maxima of the CT-FODs using the CP-decomposition, the Z-eigenvalues approach, and our proposed Adecomp-SHOT based approach. As it is required by the CP-decomposition, to decompose the fourth order CT-FOD we set the low-rank of the decomposition approximation to 2; that is, we assume that we know the number of fibers in the voxel. We recall that rather than the CP-decomposition our Adecomp-SHOT based approach does not require predefining the rank of the decomposition, that is, not require knowing a priori the number of fibers in the voxel. Moreover, to insure roughly the convergence of the ALS algorithm used by PARAFAC we initialize it by the eigenvectors of each mode of the fourth order tensor. This initialization is possible only for lower ranks up to 3, because the tensor is of dimension 3 and contains only 3 modes; for higher ranks, the initialization will be random which do not insure the convergence of the algorithm.
Figure 4 illustrates the results obtained by the CP-decomposition, the Z-eigenvectors, and the Adecomp-SHOT based approach. From Figure 4 we clearly notice that the tensor decomposition methods have a better angular resolution than the Z-eigenvalues classical method which is not able to recover crossing fibers orientations at angles below . The results show too that if we assume that the number of fibers is a priori known, then both of CP-decomposition and Adecomp-SHOT based method are equivalent in terms of angular resolution and accuracy. Nevertheless, the CP-decomposition remains limited by the constraint of predefining the rank and by the convergence of the ALS algorithm. The proposed Adecomp-SHOT based approach efficiently solves these problems without loss of angular resolution or accuracy. With an order of tensors not exceeding four and without predefining the number of fibers in a voxel, the fiber orientations are recovered with an angular resolution limit of and a mean error less than up to separation angles of .
The Effect of a Rician Noise. The most important aspect of the Adecomp-SHOT based method is its ability to render the number of fibers automatically without any assumption; therefore, we evaluate in 100 experiences on noisy synthetic data the success rate of our approach in detecting the number of crossing fiber bundles for different crossing angles. The synthetic data are corrupted with a Rician noise with different SNR levels of 40, 30, 20, and 10; the results are summarized in Table 1. Furthermore, in order to compare our results with the state-of-the-art, we conduct the same experiences for each of the CP-decomposition and the Z-eigenvalue approaches. The results in Table 1 show that up to a crossing angle of in case of SNR 30 the success rate of our method is of ; this rate slightly decreases for a crossing angle of where the correct number of fibers is rendered at , and for a crossing angle of the success rate is of . Notice that, even in case of low signal to noise ratio SNR = 20, the correct number of fibers is found with a rate of up to crossing angle of , and for a really low SNR level of 10 the success rate is higher than up to a separation angle of . These results prove that the ability of the Adecomp-SHOT method to find automatically the number of fibers is not highly sensitive to noise. Thus, our Adecomp-SHOT based method is really reliable when we aim to detect the number of underlying fibers, even in case of really low SNR levels. For comparison, results about the number of fibers rendered by the CP and the Z-eigenvalues methods are too represented in Table 1; we can notice from these results that even if the number of fibers is automatically recovered by the Adecomp-SHOT approach, the success rate does not highly differ from the success rate of the CP-decomposition where the number of fibers constitutes an input parameter.
To evaluate the effect of noise on the accuracy of our Adecomp-SHOT based approach we represent on the left column of Figures 5 and 6 the mean and the standard deviation of the mean error between the recovered fiber orientations and the ground truth fiber orientations, corresponding to different SNR levels, and on the right column of Figures 5 and 6 we represent the mean and standard deviation of the fiber weights ; figures on the center column of Figures 5 and 6 represent the mean and the standard deviation of the number of detected fibers. Once again, in order to compare our approach to other classical and state-of-the-art methods, results for each of our Adecomp-SHOT based approach CP and Z-eigenvalues approaches are illustrated on Figures 5 and 6. As in the free-noise case, Figures 4(a) and 4(b), for SNR = 40, Figures 5(a)–5(c), the method recovers effectively the entire simulated fiber orientations with an angular resolution limit equal to and up to a separation angle of the two crossing fibers are detected with a slight increase in the mean error but the mean of the mean error remains less than as it is shown by Figure 5(a). While, for SNR 30, Figures 5(d)–5(f), the angular resolution limit is somewhat reduced and its value is between and with a mean of the mean error not exceeding for a separating angle of as shown in Figure 5(d); we notice that even when the SNR level equals 30, the angular resolution limit is still better than the angular resolution limit of the classical maxima localization methods in the free-noise case. For a low signal to noise ratio of 20, Figures 6(a)–6(c), the method still recovers the fiber orientations with an angular limit between and and the mean error remains inferior to for a crossing angle of . Up to SNR level of 20 the results of each of the Adecomp-SHOT and CP methods are equivalent in terms of angular resolution and accuracy if we assume the number of fiber known. Concerning the Z-eigenvalues approach, the results on Figure 5 confirm that even in case of a high SNR level the Z-eigenvalues method clearly fails to recovers crossing fiber orientations when the separation angle is lower than . Furthermore, additional tests are performed with really low SNR level of 10 to evaluate much more the robustness to noise. The results are represented in Figures 6(d)–6(f) and show that in case of SNR 10 the angular resolution of our method is between and where the correct number of fibers is rendered at more than (Table 1) with a mean of the mean error not exceeding for a separation angle of while for the same separation angle the CP-decomposition method has a mean of the mean error higher than . These results prove that our proposed Adecomp-SHOT based approach is effective and robust to noise.
3.2. Phantom Diffusion Dataset
We conduct other experiments on the FiberCup phantom data acquired at s/mm2 with 64 acquisitions [21, 22] downloaded from the computer-assisted neuroimaging laboratory (LNAO) link: http://www.lnao.fr/spip.php?article112. The fourth order tensor is estimated from the phantom data using the CT-FOD algorithm and represented by spherical function in Figure 7(a). The recovered fiber directions are plotted on Figure 7(b). The results show that our Adecomp-SHOT based approach is able to render the fiber bundle directions.
3.3. In Vivo Human Cerebral Dataset
We conduct further tests, on real dataset obtained from the Stanford University [23] link: http://purl.stanford.edu/yx282xq2090. These data were acquired with 160 gradient directions with a -value of 2000 s/mm2; the thickness of slice is 2 mm. The CT-FODs of fourth order are estimated from these data and decomposed with the Adecomp-SHOT; Figures 8(a) and 8(b) show the fourth order CT-FODs and the fiber directions, respectively. On a coronal slice, Figure 8 shows that the method reliably extracts the Corpus Collum (CC), the Corticospinal Tractu (CST), the Cingulum (CG), and the Superior Longitudinal Fasciculus (SLF). Two regions representing the CC, CST, and SLF intersections are highlighted on Figures 8(c) and 8(d) and show the ability of our approach to extract crossing fibers with a high angular resolution.
4. Conclusion
In this paper, we have proposed an Adecomp-SHOT based approach to extract the fiber directions from DW-MRI data. Till now, the CP-decomposition approach constitutes the state-of-the-art of the tensor decomposition applied to the fiber directions search in dMRI. Although the CP-decomposition method provides a better angular resolution and accuracy than the classical methods of maxima localization, this numerical method is suboptimal and suffers from two important inconveniences: the inability to ensure the convergence of the algorithm and the requirement to predefining the rank of the decomposition, that is, the number of fibers in the voxel. Thus, to overcome the considerable limits related to the CP-decomposition approach in the diffusion MRI, we propose a novel approach based on an analytical decomposition of symmetric high order tensor to extract the fiber directions in dMRI. Unlike CP-decomposition, our proposed approach is able to recover the entire crossing fiber directions whatever its number and that without any assumption. To exploit the Adecomp-SHOT on diffusion dataset and to take into account the ground truth properties of the diffusion, we imposed a real and nonnegative constraint by using the NNLS algorithm and by introducing a change of variables. The change of variables was done by a rotation transformation to preserve the fiber directions. This coordinate transformation permitted to overcom a significant constraint imposed by the original Adecomp-SHOT. Indeed, the Adecomp-SHOT, as initially described, enables to extract directions in the affine space whose first Cartesian coordinates are zero, but by using optimally the transformation coordinates, the entire fibers directions are recovered whatever its positions in the affine space. Different validations tests were conducted on synthetic noisy diffusion data, phantom, and real data. The tests on synthetic dataset have shown three principal advantages. (1) The Adecomp-SHOT based approach overcomes the limits related to CP-decomposition without loss in the angular resolution and angular accuracy. (2) Our approach is efficient, accurate, and robust to noise: for SNR equal to 30 our Adecomp-SHOT based approach has an angular resolution limit < and up to the mean of the mean error does not exceed . (3) The rician noise does not really affect the ability of the method to detect the number of fibers in a voxel; indeed, up to a crossing angle of for a low SNR equal to 20 the correct number of fibers is found at . Finally, the tests conducted on phantom and real data confirm the ability of the method to reliably extract the directions of fiber bundles especially in regions where many fiber bundles intersect.
Appendix
As described by (7) of Section 2.2.2, the affine and minimal decomposition of a rank-4 polynomial of degree 4 in 3 variables is given by the following form: According to (A.1), we generate a order tensor of rank-4 in 3 dimensions from 4 crossing vectors separated with angle of ; these vectors , given below, are uniformly distributed on the unit sphere and weighted with the same weights . Consider The resulted tensor is represented by the polynomial given by In the following, we propose to illustrate how the Adecomp-SHOT algorithm works by decomposing the homogenous polynomial given in (A.3) step by step. will be thus decomposed in a minimal sum of 4th order rank-1 terms as described in (A.1), and at the end we expect to find the vectors given in the equation (A.2) and the weights . An affine decomposition of exists, if and only if an affine decomposition of exists; decomposing is then equivalent to decomposing , which is the unhomogenization of in the dual space; then, the first step of the Adecomp-SHOT algorithm is the computation of the coefficients of from the coefficients of .
Step 1. Consider computation of the coefficients of from the coefficients of .
From the coefficients of we compute the coefficients of as described in (A.5), with the associated linear form of in the dual space .
Consider a general form of a 4th order tensor of dimension 3 given by
From the expression of in (A.3) we deduce the coefficients of and by the relation (A.5) we calculate the coefficients of from the coefficients :
Then, the associated linear form of in the dual space is given by the following expression:
Computing by unhomogenizing according to ,
The nonhomogenous basis of monomial is given by :
And the coefficients of are the same as the coefficients of :
Step 2. Construct the Hankel matrix of size (Table 3) from ; we notice that the elements corresponding to monomials with total degree higher than the polynomial degree 4 are unknown:
with , the number of rows and columns, respectively, of the Hankel matrix , and corresponding to the monomial .
We denote by the number of unique coefficients in the tensor; is a function in the order and dimension and is calculated by
In our case and , then .
Table 2 gives a part of the Hankel matrix constructed from the coefficients of . We give below an example on how to calculate the Hankel matrix elements from the coefficients of the polynomial given by (A.8):corresponding to the monomial ; ; corresponding to the monomial ; ; corresponding to the monomial ; (unknown); corresponding to the monomial ; (unknown).
Step 3. Check if the polynomial is of rank-1. For that, we check if all the minors that we can calculate from are equal to zero. In our case, all the minors of are not zero; then we set .
Step 4. Find the rank of the tensor iteratively.
For , we compute from a square submatrix of dimension corresponding to a monomials basis connected to one of size . Consider
and its extension of dimension corresponding to the monomials basis of size , which is the extension of . Consider
Since and are totally defined, we just have to find the ranks and of each of them and verify if the rank remains stable.
= rank() = 2; = rank() = 3, ; the stability condition of the rank is not satisfied; then, we increment .
Repeat for = rank() = 3; = rank() = 4, ; the stability condition of the rank is not satisfied; then, we increment .
Repeat for rank() = 4; = rank() = 4, ; the stability condition of the rank or equivalently the commutation properties of the multiplication matrix are satisfied.
Step 5. Compute the matrices and corresponding to the monomials basis and , respectively, and the associated multiplication matrix and , for . The bases and correspond to the monomial basis multiplied by and , respectively: Constructing the multiplication matrix,
Step 6. Find the parameters such that and the matrix commute.
The commutation properties lead to 16 nonlinear equations, with only 5 nontrivial equations to solve; a solution of this problem is .
Now, we verify if the rank is effectively the rank of the tensor by verifying the last condition, which is the multiplicity of the eigenvalues of . Thus, if are simple with a random real value of , then is effectively the rank of the tensor.
; the eigenvalues are simple; then the rank of the tensor is .
Step 7. Solve the generalized eigenvalues problem for rank : calculate the eigenvalues of the common eigenvectors of the multiplication matrix , with and :
and are, respectively, the eigenvectors and the eigenvalues of . Consider
and are, respectively, the eigenvectors and the eigenvalues of .
We construct the elements of the vectors from the eigenvalues of the common eigenvectors , with the first element corresponding to the unhomogenization variable equal to 1 as follows.
If , then
For instance, ; then
By normalizing the vectors , we get
The resulted corresponds exactly to the original vectors used to construct the tensor.
Step 8. Then we resolve the linear system in :
The solution of the above system is the following:
Thus, the minimal decomposition of the symmetric fourth order tensor of rank-4 in 3 dimensions associated with the homogeneous polynomial of degree 4 in 3 variables is given as follows:
Abbreviations
| MRI: | Magnetic resonance imaging |
| dMRI: | Diffusion magnetic resonance imaging |
| DTI: | Diffusion tensor imaging |
| HARDI: | High angular resolution imaging |
| DSI: | Diffusion spectrum imaging |
| QBI: | Q-ball imaging |
| SHOT: | Symmetric high order tensor |
| ODF: | Diffusion orientation distribution function |
| FOD: | Fiber orientation distribution function |
| SH: | Spherical harmonics |
| CT-FOD: | Cartesian tensor-FOD |
| CP: | CANDECOMP/PARAFAC |
| ALS: | Alternating least squares |
| Adecomp-SHOT: | Analytical decomposition of a symmetric high order tensor |
| NNLS: | Nonnegative least squares |
| SNR: | Signal to noise ratio |
| CC: | Corpus Callosum |
| CST: | Corticospinal Tract |
| CG: | Cingulum |
| SLF: | Superior Longitudinal Fasciculus. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work is supported by USTHB in Algeria and by INRIA (International Department) in France. The authors thank both institutions for their financial support. They would like to express also their thanks to Bernard Mourrain (Inria, Galaad) and Aurobrata Ghosh (Inria, Athena) for the useful discussions and comments they had during the progress of this work.