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 socalled Cartesian tensor fiber orientation distribution (CTFOD) function, whose maxima correspond exactly to the orientations of the fibers. The fourth order symmetric positive definite tensor that represents the CTFOD 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 stateoftheart 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 Qball 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 qspace. 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 QBIODF 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 CTFOD for Cartesian tensorFOD, whose maxima correspond exactly to the orientations of the underlying fibers. By imposing efficiently the positivity constraint, the CTFOD function appears as an alternative to the FOD described in the SH basis. Furthermore, thanks to its polynomial form, the maxima of the CTFOD function can be easily located. Therefore, in our work we start by estimating the CTFOD function to reconstruct the symmetric high order tensor from the diffusion weightedMRI (DWMRI) 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 CTFOD 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 Zeigenvalues introduced by Qi in 2005 [10]. This maxima localization method, the socalled 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 rank1 tensor has been proposed to the maxima extraction issue in [12], and the lowrank 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 CPdecomposition 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 AdecompSHOT for analytical decomposition of symmetric high order tensor. This AdecompSHOT method was initially proposed by Brachat et al. in 2010 [14] but remains in the theoretical field. However, the AdecompSHOT seems a priori interesting for the fiber orientations issue because, unlike the suboptimal CPdecomposition approximation, the AdecompSHOT is an analytical one and not restricted to subgeneric ranks. Thus, rather than CPdecomposition, AdecompSHOT 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 AdecompSHOT based approach to extract the fiber orientations from diffusion weightedMRI data and we prove through many validations tests the effectiveness of the proposed method. The AdecomSHOT is based on the SHOT; therefore we propose to use the CTFOD for reconstructing the SHOT from the DWMRI signal, since the CTFOD constitutes the stateoftheart. We begin this paper by presenting first the CTFOD algorithm, before explaining the CPdecomposition and providing a detailed version of the AdecompSHOT algorithm; in this section we will also present the results of the intrinsic study done on both of AdecompSHOT algorithm and CPdecomposition. Then we describe our AdecompSHOT 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 CTFOD 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 bvalues, 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 lowrank 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 AdecompSHOT approach we propose to analytically decompose a symmetric tensor of any order and any dimension in a minimal sum of rank1 terms.
2.2.1. Numerical Method: CP LowRank Approximation
The tensor decomposition problem consists in writing a given tensor, in sum of outer product of vectors, that is, rank1 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 rank1 tensors (vectors).
However, determining the rank of high order tensors (order ) is a hard mathematical and NPcomplet problem. Therefore, a lowrank 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 rank1 tensors .
is the subgeneric rank of the tensor , for symmetric tensors . Thus, the CPdecomposition 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 rank2 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 rank1 tensors and the CPdecomposition solutions, and Figure 1(b) gives the rank1 tensors weights found by the CPdecomposition, 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 CPdecomposition 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 CPdecomposition are intrinsically limited. Furthermore, the CPdecomposition has another important limit which is the requirement to predefine the rank of the decomposition.
The use of the CPdecomposition 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 lowrank value (rank ) in order to extract two crossing fiber orientations. The extracted fiber orientations correspond to the rank1 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 CPdecomposition 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: AdecompSHOT. Rather than the CP numerical approach, the AdecompSHOT method gives a minimal decomposition without any apriority. However, to date the AdecompSHOT 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 AdecompSHOT applied to the fiber orientations search in dMRI would be more interesting than the CP lowrank approximation method.
2.2.2. Analytical Method: AdecompSHOT
The AdecompSHOT 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 AdecompSHOT algorithm is able to decompose a symmetric tensor of any order and any dimension, in a minimal sum of rank1 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 rank1 linear forms [14] as follows: with representing scalars weights, rank1 linear forms in with the coefficients , and the minimal rank of , that is, of the tensor.
is a rank1 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 rank4, illustrating how the AdecompSHOT 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 AdecompSHOT algorithm; thus, rank2 fourth order tensors were constructed from a sum of two rank1 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 rank1 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 CPdecomposition, 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 rank1 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 CPdecomposition. 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 rank3 symmetric fourth order tensor constructed from 3 crossing rank1 tensors, according to (7), at angles decreasing from to . Figure 3 shows that the AdecompSHOT method gives once again a correct decomposition with a zero mean error whatever the separation angles between the 3 origin rank1 tensors. Other tests on symmetric fourth order tensors of rank4 and rank5 have been conducted; the results of these experiences show that the AdecompSHOT method still gives the correct decomposition in case of rank4 while for a rank5 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 CTFOD using the AdecompSHOT. 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 rank1 linear forms: with minimal representing the tensor rank; are rank1 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 rank1 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 rank1 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 rank1 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 AdecompSHOT 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 rank1 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 wellknown 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/mm^{2}; 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 CPdecomposition based method and the Zeigenvalues 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 CPDecomposition and ZEigenvalues Based Approaches. In order to compare the angular resolution of our fiber extraction approach to the one of the CPdecomposition and the Zeigenvalues based methods, we have conducted experiments on same noisefree synthetic diffusion dataset. From these data we reconstruct the fourth order CTFODs and then we extract the maxima of the CTFODs using the CPdecomposition, the Zeigenvalues approach, and our proposed AdecompSHOT based approach. As it is required by the CPdecomposition, to decompose the fourth order CTFOD we set the lowrank 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 CPdecomposition our AdecompSHOT 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 CPdecomposition, the Zeigenvectors, and the AdecompSHOT based approach. From Figure 4 we clearly notice that the tensor decomposition methods have a better angular resolution than the Zeigenvalues 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 CPdecomposition and AdecompSHOT based method are equivalent in terms of angular resolution and accuracy. Nevertheless, the CPdecomposition remains limited by the constraint of predefining the rank and by the convergence of the ALS algorithm. The proposed AdecompSHOT 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 AdecompSHOT 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 stateoftheart, we conduct the same experiences for each of the CPdecomposition and the Zeigenvalue 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 AdecompSHOT method to find automatically the number of fibers is not highly sensitive to noise. Thus, our AdecompSHOT 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 Zeigenvalues 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 AdecompSHOT approach, the success rate does not highly differ from the success rate of the CPdecomposition where the number of fibers constitutes an input parameter.
To evaluate the effect of noise on the accuracy of our AdecompSHOT 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 stateoftheart methods, results for each of our AdecompSHOT based approach CP and Zeigenvalues approaches are illustrated on Figures 5 and 6. As in the freenoise 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 freenoise 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 AdecompSHOT and CP methods are equivalent in terms of angular resolution and accuracy if we assume the number of fiber known. Concerning the Zeigenvalues approach, the results on Figure 5 confirm that even in case of a high SNR level the Zeigenvalues 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 CPdecomposition method has a mean of the mean error higher than . These results prove that our proposed AdecompSHOT 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/mm^{2} with 64 acquisitions [21, 22] downloaded from the computerassisted neuroimaging laboratory (LNAO) link: http://www.lnao.fr/spip.php?article112. The fourth order tensor is estimated from the phantom data using the CTFOD 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 AdecompSHOT 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/mm^{2}; the thickness of slice is 2 mm. The CTFODs of fourth order are estimated from these data and decomposed with the AdecompSHOT; Figures 8(a) and 8(b) show the fourth order CTFODs 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 AdecompSHOT based approach to extract the fiber directions from DWMRI data. Till now, the CPdecomposition approach constitutes the stateoftheart of the tensor decomposition applied to the fiber directions search in dMRI. Although the CPdecomposition 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 CPdecomposition 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 CPdecomposition, our proposed approach is able to recover the entire crossing fiber directions whatever its number and that without any assumption. To exploit the AdecompSHOT 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 AdecompSHOT. Indeed, the AdecompSHOT, 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 AdecompSHOT based approach overcomes the limits related to CPdecomposition 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 AdecompSHOT 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 rank4 polynomial of degree 4 in 3 variables is given by the following form: According to (A.1), we generate a order tensor of rank4 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 AdecompSHOT 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 rank1 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 AdecompSHOT 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 rank1. 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