recent work 
 shows that diffusion magnetic resonance imaging ( dmri ) 
 can help recovering 
 complex white matter fiber bundles .
however this is 
 still an open problem due 
 to the structural complexity of the fiber bundles , 
 which can have crossing 
 configurations .
diffusion tensor imaging ( dti ) 
  is restricted in these 
 conditions due to the hypothesis that the diffusion within 
 a voxel follows a 
 gaussian distribution , a model that can not 
 model intravoxel crossings .
q - ball 
 imaging ( qbi ) , 
 a recent high - angular resolution diffusion imaging ( hardi ) 
 technique , overcomes this limitation by 
 reconstructing a diffusion orientation 
 distribution function ( odf ) , a spherical 
 function that has its maxima agreeing 
 with the underlying fiber populations . 

the odf reconstruction from qbi is 
 attractive because it is model - free and has 
 been recently shown possible with a 
 regularized and analytical solution 
 , 
 which produces a robust and very fast 
 odf reconstruction .
in fact , the odf estimation is , 
 in practice , as fast as a 
 standard least - square diffusion tensor ( dt ) estimation .
efficient segmentation of fiber tracts in dmri images 
 is an important problem in neuroimaging 
 problem because it has many potential 
 applications .
for example , it could 
 potentially provide important information 
 on diseases that affect fiber tracts . 

alteration of the fiber tracts may 
 provide new biomarkers in white matter 
 pathologies and segmentation of these 
 tracts can also improve our understanding of 
 the functional role these tracts 
 have and the cognitive consequences of their disruption .
the goal of this work is to provide a segmentation 
 method that can separate the main white 
 matter fiber bundles in the brain .
we 
 propose a new method that can segment fiber 
 bundles and deal with fiber crossings 
 while also requiring a minimum number of hypothesis 
 from the data and a small 
 number of algorithmic parameters .
spectral 
 embedding and clustering methods 
 have recently proved to be effective in image segmentation 
 [ 4 , 5 ] . however , 
 classical approaches require restrictive hypotheses 
 that are difficult to meet 
 in real applications .
for instance , n - cuts 
  and 
 laplacian eigenmaps 
 require data within each cluster to be 
 uniformly sampled , which produces 
 artifacts when this hypothesis is not met . 

moreover , classical approaches for 
 image segmentation also assume that the scale 
 within each cluster is the same 
 using a single - scale parameter for the whole dataset . 
 in order to overcome 
 these limitations
this method looses the dependence on the 
 sampling of the elements to 
 cluster .
moreover , we propose to use an 
 adaptive scale - space parameter in order 
 to deal with space - scale differences across 
 different clusters .
finally , we 
 propose two approaches to automatically 
 determine the number of clusters by 
 analyzing the spectra of the image embedding .
another contribution of this paper is to show that the 
 q - ball odf clustering using diffusion maps 
 can reproduce the dt clustering 
 using n - cuts on simple synthetic images 
 without crossings . on more complex data 
 with crossings ,
we show that our method 
 succeeds to separate fiber bundles and 
 crossing regions on synthetic data , where 
 the dt - based methods generate 
 artifacts and exhibit wrong number of clusters . 

finally , we successfully 
 segment the fiber bundles in a real - human 
 brain dataset in different regions 
 with fibers crossing .
the main goal 
 of this work is to produce a segmentation 
 algorithm able to segment white 
 matter fiber bundles from dmri data . 
 in order to represent intravoxel 
 crossings with the odf
, we need at least 15 
 real coefficients when a spherical 
 harmonic basis is used [ 3 , 8 , 9 ] .
this leads to 3d 
 images with a high dimensional 
 element at each voxel .
this high dimensionality 
 makes previous diffusion 
 imaging segmentation approaches based on level 
 set methods such as [ 1012 ] 
 computationally expensive .
moreover , these methods 
 require an initialization 
 step . in order to perform the segmentation 
 in an initialization - free manner and 
 with a lower - dimensionality image
, 
 we use spectral clustering methods 
 [ 4 , 
 5 ] 
 which perform dimensionality reduction 
 before performing the segmentation and 
 do not need initialization . 
 the segmentation
is then performed on the 
 statistics within each cluster and the fiber 
 crossings can be identified . in this section , we present the three main parts of 
 our algorithm .
first , the estimation of the 
 q - ball diffusion odf and its 
 compact representation using spherical harmonics . 

second , the metric used to 
 measure distances between q - ball odfs . 

last , the diffusion aps spectral 
 clustering technique used to segment the 
 odf image into the background and the 
 different fiber bundles .
qbi 
 reconstructs the diffusion odf directly from the 
 hardi measurements on a single 
 sphere by the funk - radon transform ( frt ) . 
 in practice
, the frt value at a given 
 spherical point is the great circle integral of 
 the signal on the sphere 
 defined by the plane through the origin with 
 normal vector .
the odf is 
 intuitive because it has its maximum(a ) 
 aligned with the underlying population 
 of fiber(s ) .
however , computing statistics 
 on a large number of discrete odf 
 values on the sphere is computationally heavy 
 and infeasible to integrate into a 
 segmentation algorithm of the whole brain . 
 a more compact representation of the 
 odf is thus needed . in
[ 3 , 8 , 9 , 13 ] a 
 simple analytic spherical harmonic ( sh ) 
 reconstruction of the odf is proposed . 
 for completeness of the article , we now 
 review and develop the main parts of our 
 regularized analytical odf 
 reconstruction solution .
the idea is to 
 first estimate hardi signal on the 
 sphere with a regularized spherical harmonics 
 approximation and then do a 
 simple linear transformation of the 
 harmonics to obtain the desired regularized 
 odf .
the sh , 
 normally indicated by ym (  denotes the order and m the phase 
 factor ) , are a basis for complex functions on the unit sphere .
explicitly , they 
 are given as follows : 
 ( 1)ym(,)=2+14(m)!(+m)!pm(cos)eim , 
 where ( , ) obey physics 
 convention ( [0,],[0,2 ] ) and pm is an 
 associated legendre polynomial . for k=0,2,4,  , and m=k,  ,0,  ,k
, we define the new index j:=j(k , m)=(k2+k+2)/2+m and define our 
 modified basis y with elements yj such that 
 ( 2)yj={2re(ykm)if  km<0,yk0if  m=0,2img(ykm)if  0<mk , 
 where re(ym ) and img(ym ) represent the real and imaginary parts of ym , respectively . 

symmetry is 
 ensured by choosing only even order sh and the ratios in front of each term 
 also ensure that the modified basis is real and orthonormal with respect to the 
 inner product f , g=f*gd , where  denotes 
 integration over the unit sphere and f * is the complex 
 conjugate of f for f and g complex 
 functions on the sphere .
we thus approximate the signal at each of the n gradient 
 directions i as 
 ( 3)s(i,i)=j=1rcjyj(i,i ) , 
 where r=(+1)(+2)/2 is the number 
 of terms in the modified sh basis y of order . letting s be the n  1 vector 
 representing the input signal for every encoding gradient direction , c the r  1 vector of 
 sh coefficients cj , and b is the n  r matrix constructed 
 with the discrete modified sh basis 
 ( 4)b=(y1(1,1)y2(1,1)yr(1,1)y1(n,n)y2(n,n)yr(n,n ) ) . 

we can write 
 the set of equations as an overdetermined linear system s = bc . we want to solve for the sh series coefficients cj , where cj=s(,)yj(,)d. at this point , instead of simply evaluating the 
 integrals directly as done in   or performing a straightforward least - squared 
 minimization as in [ 15 , 16 ] , 
 we add local regularization directly into our 
 fitting procedure .
this is to be able to 
 use a high - order estimation without 
 overmodeling the small perturbations due 
 to noise in the input diffusion mri 
 signal .
we thus define a measure , e , of the deviation from smoothness of a 
 function f defined on the 
 unit sphere as e(f)=(bf)2d , where b is the 
 laplace - beltrami operator . using the orthonormality of the modified sh basis , 
 where we have yi(,)yj(,)d=ij , the above functional e can be 
 rewritten straightforwardly [ 3 , 13 ] as 
 ( 5)e(f)=b(pcpyp)b(qcqyq)d=j=1rcj2(j)2((j)+1)2=ctlc , 
 where l is simply the rr matrix with 
 entries (j)2((j)+1)2 along the 
 diagonal ( (j ) is the order associated with the jth coefficient , 
 that is , for j=1,2,3,4,5,6,7,  (j)=0,2,2,2,2,2,4 ,  ) .
therefore , the 
 quantity we wish to minimize can be expressed in matrix form as 
 ( 6)m(c)=(sbc)t(sbc)+ctlc , 
 where  is the weight 
 on the regularization term .
the coefficient vector minimizing this expression 
 can then be determined just as in the standard least - squares fit ( =0 ) from which we obtain the generalized expression 
 for the desired spherical harmonic series coefficient vector 
 ( 7)c=(btb+l)1bts . 
 from this sh 
 coefficient vector we can recover the signal on the q - ball for any ( , ) as s(,)=j=1rcjyj(, ) .
intuitively , 
 this approach penalizes an approximation function for having   higher - order terms 
 in its modified sh series .
this eliminates most of the higher - order terms due to 
 noise while leaving those that are necessary to describe the underlying 
 function .
however , obtaining this balance depends on choosing a good value for 
 the parameter . we use the l - curve numerical method   and 
 experimental simulations to determine a good smoothing parameter [ 3 , 13 , 18 ] . 
 here
, =0.006 is used as in [ 3 , 13 , 18 ] . the true 
 diffusion orientation distribution function ( odf ) in a unit direction u , (u ) , is given by the radial projection of the probability 
 distribution function ( pdf ) of the diffusing water molecule .
tuch   showed 
 that this diffusion odf could be estimated directly from the raw hardi signal s on a single 
 sphere of q - space by the funk - radon transform ( frt ) ( figure 1 ) . in [ 3 ,
13 ] , we showed how this frt can be evaluated 
 analytically with an elegant corollary to the funk - hecke theorem .
the 
 final odf reconstruction on the sphere then becomes a simple linear 
 transformation of the sh coefficients cj describing the 
 input hardi signal s , 
 ( 8)(,)=j=1r2p(j)(0)cjfjyj(, ) , 
 where fj are the sh 
 coefficients describing the odf  and p(j)(0)=(1)/2(135((j)1)/246(j ) ) because (j ) is always even 
 in our modified sh basis .
we see that the shs are eigenfunctions of the 
 funk - radon transform with eigenvalues depending only on the order  of the sh 
 series .
hence , by using an sh estimation of the hardi signal , 
 we have showed that the qbi can be solved analytically .
this was also showed in 
 [ 8 , 9 ] . an important contribution in favor of our approach
is that this solution 
 can be obtained while imposing a well - defined regularization criterion .
the 
 accuracy of the modified sh series approximation with the laplace - beltrami 
 smoothing was established in   and our regularized odf solution was also 
 shown to have better fiber detection properties and shown to be more robust 
 to noise than similar solutions [ 8 , 9 ] . once the odf 
 are computed , we want to capture similarities and 
 dissimilarities between two 
 odfs , that is , two spherical functions ,s2 that can be 
 represented by real - sh vectors of length r , f={f1,  ,fr } and f={f1,  ,fr}r , as shown in ( 8) in the previous section .
since 
 the odfs come from real physical diffusion 
 measurements they are bounded and 
 form an open subset of the space of real - valued 2 spherical 
 functions with an inner product , defined as 
 ( 9),=(,)(,)d=(i=1rfiyi(,)j=1rfjyj(,))d. 
 again , because 
 of the orthonormality of the spherical harmonic basis , the cross - terms cancel 
 and the expression is simply 
 ( 10),=j=1rfjfj. 
 therefore , the 
 induced 2 norm =, giving us the 
 distance metric between two odfs is 
 ( 11)=j=1r(fjfj)2 .
the euclidean 
 distance was also used successfully for odf segmentation in   and for dti 
 segmentation in   even though more appropriate metrics exist such as the 
 j - divergence [ 11 , 20 ] and riemannian geodesic distances .
for instance , since the odf 
 can be viewed as a probability distribution function ( pdf ) of fiber 
 orientations , one can use the kullback - leibler distance between two pdfs , as 
 done in .
however , in that case the problem quickly blows up computationally 
 because one needs to use all n discrete hardi 
 data on the sphere instead of the r sh coefficients 
 ( rn ) .
we now want to 
 segment white matter fiber bundles in a q - ball image . 

one of the open questions 
 in q - ball image analysis and clustering is 
 that which metric should be used to 
 compare q - ball odfs . here , we describe a 
 clustering algorithm that infers an 
 embedding and a metric to compare odf images . 

we derive an affinity measure 
 incorporating the euclidean distance and the spatial 
 location distance between 
 odfs .
as 
 mentioned in , 
 the euclidean distance within this embedding actually 
 represents an intrinsic metric of the data , 
 which can be used to perform 
 statistics in the embedded space and can thus be used to 
 segment q - ball odf 
 images into white matter fiber bundles . in recent 
 years , spectral manifold learning and clustering techniques [ 4 , 6 , 2123 ] 
 have become one of the most popular modern clustering 
 family of methods .
they 
 are simple to implement , 
 they can be solved efficiently by standard linear 
 algebra software , and they very often outperform 
 traditional manifold learning 
 and clustering algorithms such as the classical principal component 
 analysis 
 ( pca )   and k - means   algorithms . moreover , due to the dimensionality 
 reduction properties , they are especially well suited to work with 
 high - dimensional data .
these techniques have been recently used to cluster 
 various types of images [ 4 , 
 5 ] and white matter fiber tracts . in our case
, 
 we perform the spectral clustering for two different types of elements : the dt 
 and the odf . in the dt case ,
the element is represented by a 6-dimensional 
 vector corresponding to the upper ( or lower ) 
 triangular part of the dt 33 
 symmetric matrix . in the odf case ,
the element is represented by the 15-dimensional 
 vector corresponding to the 4th - order spherical 
 harmonic odf estimation .
each element 
 to be clustered is represented as a 
 node in a graph and the edges joining the 
 vertex are a measure of affinity 
 between the elements .
this affinity measure lies between 0 and 1 , 0 being the less 
 affine case . a spectral decomposition of this 
 graph is taken by calculating the 
 eigenvalue decomposition ( evd ) of the graph 
 laplacian .
then a 
 low - dimensional euclidean manifold embedding is inferred from this decomposition . 

preservation of 
 the distance relationship : after a distance is defined between elements , the 
 learned manifold should preserve the distance relation.uniform 
 sampling of the elements : the density of the extracted elements changes if and 
 only if these elements belong to anatomically different bundles.convexity of 
 the elements : if two elements are in the dataset , almost all of the 
 intermediate tracts obtained by the interpolation that can be inferred from the 
 metric used to build the affinity matrix are in the dataset . 
  preservation of 
 the distance relationship : after a distance is defined between elements , the 
 learned manifold should preserve the distance relation .
uniform 
 sampling of the elements : the density of the extracted elements changes if and 
 only if these elements belong to anatomically different bundles .
convexity of 
 the elements : if two elements are in the dataset , almost all of the 
 intermediate tracts obtained by the interpolation that can be inferred from the 
 metric used to build the affinity matrix are in the dataset .
it is not easy to guarantee that the data to be 
 embedded and clustered will adhere to these hypotheses . 

donoho and grimes , in , 
 analyze when a spectral embedding algorithm is 
 able to recover the true 
 parameterization of a set of images . 

as medical images represent the 
 discretization of a continuous space , 
 hypotheses 1 and 3 are plausible .
however , there is no indication that 
 within a fiber bundle the 
 distribution of the elements ( dt or odf ) are 
 uniformly sampled .
moreover , in it is shown that different sampling 
 frequencies within one cluster leads the n - cuts and 
 laplacian eigenmaps methods 
 to subdivide the cluster in several parts . in order 
 to overcome this limitation 
 and to be resilient to sampling frequency 
 differences within a cluster , we use 
 the diffusion maps   spectral embedding technique .
we now describe the three 
 steps involved in the diffusion maps algorithm in turn .
step 1 ( computing the affinity matrix).letting x represent the 
 set of all odf elements to cluster , the main idea is to look for a 
 representation between the elements of x that is more 
 representative than r ( recall that 
 odfs are r ) and reduces 
 the dimensionality of the problem . with keeping this in mind , 
 a fairly good way of 
 representing any set of elements with an affinity function a : xx>0 , is a weighted graph , g(x , e , w( ) ) , where the 
 weight of the edge between two vertices represents the affinity of the elements 
 connected by this edge .
more formally , for an edge , 
 e=(fi , fj)e , the weight of the edge is w(e)=a(fi , fj ) .
hence , each element of the adjacency matrix of g or conversely 
 the affinity matrix of ( x , a( ) ) is 
 ( 12)aij = a(fi , fj ) . 
 taking this in 
 account
, the weighted graph g(x , e , w( ) ) can be also 
 noted as g(x , a ) .
letting x represent the 
 set of all odf elements to cluster , the main idea is to look for a 
 representation between the elements of x that is more 
 representative than r ( recall that 
 odfs are r ) and reduces 
 the dimensionality of the problem . with keeping this in mind , 
 a fairly good way of 
 representing any set of elements with an affinity function
a : xx>0 , is a weighted graph , g(x , e , w( ) ) , where the 
 weight of the edge between two vertices represents the affinity of the elements 
 connected by this edge .
more formally , for an edge , 
 e=(fi , fj)e , the weight of the edge is w(e)=a(fi , fj ) .
hence , each element of the adjacency matrix of g or conversely 
 the affinity matrix of ( x , a( ) ) is 
 ( 12)aij = a(fi , fj ) . 
 taking this in 
 account
, the weighted graph g(x , e , w( ) ) can be also 
 noted as g(x , a ) .
usually , a distance function d( ) instead of an 
 affinity function is given .
the distances can be easily converted into 
 affinities by applying a kernel to the distance function 
 ( 13)a(fi , fj)=e(d(fi , fj)2/ij2 ) , 
 where  is an adaptive 
 scale - space parameter that may depend on the elements fi and fj . in this work ,
the adaptive scale - space parameter is 
 taken following . 
 a  neighbor - number
 k is given as 
 parameter to the algorithm and then ij2=d(fi , fik)d(fj , fjk ) , where fik is the kth closest 
 neighbor according to the distance function d(, ) of element fi .
this choice of a scaling parameter for each point 
 allows self - tuning of the point - to - point distances according to the local 
 statistics of the neighborhoods surrounding points i and j. as in image segmentation , the spatial position of each 
 element is important , the spatial dependency should be incorporated into the 
 affinity matrix . following [ 5 , 
 31 ]
, we use markovian relaxation to incorporate this 
 information . in order to represent the affinity of all the elements that can be 
 reached within one spatial step , the affinity matrix is modified in the 
 following way : 
 ( 14)a1={aijif  coords(fi)coords(fj)21,0 in  any other case , 
 where coords   ( f ) are spatial 
 coordinates of element f in the image 
 ( 15)p1=1maxld(a1)ll{maxld(a1)lld(a1)iiif  i = j , a1  in  any other case , 
 where d(a1 ) is a diagonal 
 matrix with d(a1)ii=ja1ij , usually called the row - sum or 
 degree matrix of a1 .
then , obtaining the affinities of elements that can be 
 reached within s spatial steps 
 is enough to elevate p1 to the power of s , ps=(p1)s as stated in . 

moreover , s can be chosen 
 to be the smallest positive integer which results in nonzero elements in the 
 whole matrix in order to represent the weakest connected induced graph .
the 
 diagonal adjustment forces the inherent random walk to a uniform steady state , 
 hence every part of the markov field will be explored at the same speed . for 
 the sake of clarity ,
ps will be 
 referred to as affinity matrix a in the rest of the 
 paper .
step 2 ( performing the embedding).the algorithm 
 must embed the elements of x into an n - dimensional 
 euclidean space y(x ) .
this is done by applying eigenvalue decomposition to 
 the laplacian of the affinity matrix .
as in [ 6 , 7 , 27 ] , this is done by performing the 
 spectral decomposition of the graph laplacian of the graph induced by a , 
 ( 16)=d(a)a|x||x| , 
 where |x| is number of 
 elements to be clustered.in order to overcome the necessity of hypothesis 2 , 
 we prenormalize the affinity matrix as done in 
 . 

this is done by normalizing the weight of each edge of the graph , 
 aij , by the probability density of both elements relating 
 through the edge , 
 ( 17)(ap)ij = aijp(i)p(j ) , 
 where p( ) , the probability density function of the elements in x , is not known but can be approximated up to a 
 multiplication factor by(18)p(i)=kaik=kaki.due to the necessity of having a uniform behavior of 
 the clustering algorithm without minding the scale of the affinity measure 
 taken , a doubly stochastic matrix normalization is performed : 
 ( 19)ads = d(ap)1/2apd(ap)1/2 |x||x| . 
 as ads is a double 
 stochastic symmetric matrix , the eigenvalue decomposition 
 of ( 16 ) can 
 be calculated by taking the singular value decomposition ( svd ) 
 ( 20)vsvt = ads|x||x| . 
 finally , the euclidean coordinates yi of an element fix in the n - dimensional 
 embedding manifold are 
 ( 21)y(fi)=yi=1vi0(1vi1,  ,nvin)t ,  fix , 
 where 
 ( 22)v=(v0v|x|1)|x||x| 
 is the 
 eigenvector column matrix and the corresponding eigenvalues are , 1=01|x|10 .
the first eigenvector v0 is not taken into account as a component in the 
 embedding because it is constant and hence meaningless as shown in [ 6 , 7 , 27 ] .
the algorithm 
 must embed the elements of x into an n - dimensional 
 euclidean space y(x ) .
this is done by applying eigenvalue decomposition to 
 the laplacian of the affinity matrix .
as in [ 6 , 7 , 27 ] , this is done by performing the 
 spectral decomposition of the graph laplacian of the graph induced by a , 
 ( 16)=d(a)a|x||x| , 
 where |x| is number of 
 elements to be clustered . in order to overcome the necessity of hypothesis 2 , 
 we prenormalize the affinity matrix as done in 
 . 

this is done by normalizing the weight of each edge of the graph , 
 aij , by the probability density of both elements relating 
 through the edge , 
 ( 17)(ap)ij = aijp(i)p(j ) , 
 where p( ) , the probability density function of the elements in x , is not known but can be approximated up to a 
 multiplication factor by(18)p(i)=kaik=kaki . due to the necessity of having a uniform behavior of 
 the clustering algorithm without minding the scale of the affinity measure 
 taken , a doubly stochastic matrix normalization is performed : 
 ( 19)ads = d(ap)1/2apd(ap)1/2 |x||x| . 
 as ads
is a double 
 stochastic symmetric matrix , the eigenvalue decomposition 
 of ( 16 ) can 
 be calculated by taking the singular value decomposition ( svd ) 
 ( 20)vsvt = ads|x||x| . 

finally , the euclidean coordinates yi of an element fix in the n - dimensional 
 embedding manifold are 
 ( 21)y(fi)=yi=1vi0(1vi1,  ,nvin)t ,  fix , 
 where 
 ( 22)v=(v0v|x|1)|x||x| 
 is the 
 eigenvector column matrix and the corresponding eigenvalues are , 1=01|x|10 .
the first eigenvector v0 is not taken into account as a component in the 
 embedding because it is constant and hence meaningless as shown in [ 6 , 7 , 27 ] .
step 3 ( clustering).once the 
 embedding has been performed , several techniques have been proposed for the 
 clustering step [ 4 , 6 , 32 ] . once the 
 embedding has been performed , several techniques have been proposed for the 
 clustering step [ 4 , 6 , 32 ] .  
the first step in this process is to determine the 
 number of clusters , this can be done in two ways .
the first , as in , is 
 choosing the number of clusters according to the  elbow  .
for instance , if the slope of the eigenvalue plot changes 
 noticeably at eigenvector i , the number of clusters should be i+1 .
the second way is reordering the affinity matrix 
 rows and columns following the second eigenvector as proved in , which 
 shows the block structure of the matrix as squared blocks along the matrix 
 diagonal .
their commended 
 number of dimensions for the embedding is the same as the number of clusters . 

finally , the clustering is performed by running a k - means clustering algorithm 
 on the space spanned by y(x ) .
a formal justification for this approach can be 
 found in [ 6 , 32 ] .
synthetic datawe generate 
 synthetic hardi data using the multitensor model which is simple and leads to 
 an analytical expression of the odf [ 2 , 18 ] . for a given b - factor and 
 noise level
, we generate the diffusion - weighted signal 
 ( 23)s(ui)=k=1n1nexp(b  uitdk()ui)+noise , 
 where ui is the ith gradient direction on the 
 sphere , n is the number 
 of fibers , and 1/n is the volume 
 fraction of each fiber . in practice , we use n=81 from a 3rd - order tessellation of the 
 icosahedron , b=3000  s / mm , and n=1 or 2 .
dk( ) is the 
 diffusion tensor with standard eigenvalues [ 3,3,1.7 ] 10  mm / s oriented in direction  , which agree with reported physiological values . finally , we add complex gaussian noise with standard 
 deviation of 1/35 , producing a signal with signal - to - noise ratio of 35 .
we generate 
 synthetic hardi data using the multitensor model which is simple and leads to 
 an analytical expression of the odf [ 2 , 18 ] . for a given b - factor and 
 noise level
, we generate the diffusion - weighted signal 
 ( 23)s(ui)=k=1n1nexp(b  uitdk()ui)+noise , 
 where ui is the ith gradient direction on the 
 sphere , n is the number 
 of fibers , and 1/n is the volume 
 fraction of each fiber . in practice , we use n=81 from a 3rd - order tessellation of the 
 icosahedron , b=3000  s / mm , and n=1 or 2 .
dk( ) is the 
 diffusion tensor with standard eigenvalues [ 3,3,1.7 ] 10  mm / s oriented in direction  , which agree with reported physiological values .
finally , we add complex gaussian noise with standard 
 deviation of 1/35 , producing a signal with signal - to - noise ratio of 35 .
we generate three synthetic data example , two simple 
 examples : one with a ring of sinusoidal - shaped 
 fibers , one with fibers with 
 different sizes and scales , and the other with complex crossing areas simulating 
 the  u-fibers ( corticocortical fibers ) 
 that can occur in the brain . these 
 synthetic datasets help understand the behavior of the different spectral 
 clustering methods when confronted with simple and complex fiber geometries . human brain datadiffusion - weighted data and high - resolution t1-weighted images 
 were acquired on a 
 whole - body 3 tesla magnetom trio scanner 
 ( siemens , erlangen ) equipped with an 
 8-channel head array coil 
 . 
 the spin - echo echo - planar - imaging sequence , te=100  ms , tr=12  s , 128  128 image matrix , 
 fov = 220  220  mm , 
 consists of 60 diffusion 
 encoding gradients  
 with a b - value of 1000  s / mm . seven images without any 
 diffusion weightings are placed at the beginning 
 of the sequence and after each 
 block of 10 diffusion - weighted images as anatomical 
 reference for offline 
 motion correction .
the measurement of 72 slices 
 with 1.7  mm thickness ( no gap ) 
 covered the whole brain .
random variations 
 in the data were reduced by 
 averaging 3 acquisitions , 
 resulting in an acquisition time of about 45 minutes . 

additionally , fat saturation was employed 
 and we used 6/8 partial fourier imaging , 
 a hanning window filtering , and 
 parallel acquisition ( generalized autocalibrating 
 partially parallel 
 acquisitions , reduction factor = 2 ) in the axial plane .
diffusion - weighted data and high - resolution t1-weighted images 
 were acquired on a 
 whole - body 3 tesla magnetom trio scanner 
 ( siemens , erlangen ) equipped with an 
 8-channel head array coil 
 . 
 the spin - echo echo - planar - imaging sequence , te=100  ms , tr=12  s , 128  128 image matrix , 
 fov = 220  220  mm , 
 consists of 60 diffusion 
 encoding gradients  
 with a b - value of 1000  s / mm .
seven images without any 
 diffusion weightings are placed at the beginning 
 of the sequence and after each 
 block of 10 diffusion - weighted images as anatomical 
 reference for offline 
 motion correction .
the measurement of 72 slices 
 with 1.7  mm thickness ( no gap ) 
 covered the whole brain .
random variations 
 in the data were reduced by 
 averaging 3 acquisitions , 
 resulting in an acquisition time of about 45 minutes . 

additionally , fat saturation was employed 
 and we used 6/8 partial fourier imaging , 
 a hanning window filtering , and 
 parallel acquisition ( generalized autocalibrating 
 partially parallel 
 acquisitions , reduction factor = 2 ) in the axial plane .
the brain is peeled from the t1-anatomy , which was 
 aligned with the talairach stereotactical coordinate 
 system .
the 21 images 
 without diffusion weightings distributed 
 within the whole sequence were used to 
 estimate motion correction parameters 
 using rigid - body transformations 
 , 
 implemented in . 
 the motion correction for the 180 diffusion - weighted 
 images
was combined with a global registration 
 to the t1 anatomy computed with 
 the same method .
the 
 registered images were interpolated to the new 
 reference frame with an isotropic voxel 
 resolution of 1.72  mm and the 3 
 corresponding acquisitions and gradient 
 directions were averaged .
distance functions between elements to clusterin order to 
 implement the diffusion maps 
 spectral clustering method a distance function for 
 each data type is chosen .
this distance 
 function is used to calculate the 
 affinity matrix as expressed by 
 ( 13 ) . 
 in the dt case , following , 
 we use the riemannian tensor distance . in the odf case
we use the distance 
 shown in ( 11 ) . in order to 
 implement the diffusion maps 
 spectral clustering method a distance function for 
 each data type
this distance 
 function is used to calculate the 
 affinity matrix as expressed by 
 ( 13 ) . 
 in the dt case , following , 
 we use the riemannian tensor distance . in the odf case
diffusion maps versus n - cutsthe first 
 experiment shows the difference in performance 
 between the diffusion maps and 
 n - cuts approach . the n - cut algorithm does not 
 perform the sampling - based 
 normalization described by 
 ( 17 ) and is thus 
 sensitive to sampling frequency differences 
 within the clusters . in order to 
 show this sampling hypothesis problem , 
 we used a ring fiber bundle with 
 different sampling frequencies . within the 
 ring ,
the fibers have a sinusoidal 
 shape and the frequency of the modulating sine 
 function is 4 times bigger in 
 the lower half of the ring .
more formally , 
 the fibers follow the angular 
 function o()=+(1/8)sin(),0<2 , where =8 for the upper 
 half of the ring and =32 for the lower 
 half .
the results of 
 both clustering techniques are shown in figure 2 , where the background 
 has been masked out .
figure 2(a ) 
 shows the plot of the first 10 
 eigenvalues for the n - cuts method , shown in 
 figures 2(b ) and 2(c ) , and 
 the slope between the line joining 0 and 1 and the line 
 joining 1 and 2 changes 
 drastically .
figure 2(d ) 
 shows the plot of the first 10 
 eigenvalues for the diffusion maps method whose clustering results are shown in 
 figure 2(e ) .
the n - cuts exhibits frequency - dependent clustering 
 artifacts while the diffusion maps method clearly shows two clusters . in the 
 diffusion maps , the clustering has correctly 
 segmented the background and the 
 ring .
the first 
 experiment shows the difference in performance 
 between the diffusion maps and 
 n - cuts approach . the n - cut algorithm does not 
 perform the sampling - based 
 normalization described by 
 ( 17 ) and is thus 
 sensitive to sampling frequency differences 
 within the clusters . in order to 
 show this sampling hypothesis problem
, 
 we used a ring fiber bundle with 
 different sampling frequencies . within the 
 ring ,
the fibers have a sinusoidal 
 shape and the frequency of the modulating sine 
 function is 4 times bigger in 
 the lower half of the ring . more formally , 
 the fibers follow the angular 
 function o()=+(1/8)sin(),0<2 , where =8 for the upper 
 half of the ring and =32 for the lower 
 half .
the results of 
 both clustering techniques are shown in figure 2 , where the background 
 has been masked out .
figure 2(a ) 
 shows the plot of the first 10 
 eigenvalues for the n - cuts method , shown in 
 figures 2(b ) and 2(c ) , and 
 the slope between the line joining 0 and 1 and the line 
 joining 1 and 2 changes 
 drastically .
figure 2(d ) 
 shows the plot of the first 10 
 eigenvalues for the diffusion maps method whose clustering results are shown in 
 figure 2(e ) .
the n - cuts exhibits frequency - dependent clustering 
 artifacts while the diffusion maps method clearly shows two clusters . in the 
 diffusion maps , the clustering has correctly 
 segmented the background and the 
 ring .
odf versus dt imagesin figure 3 , a single 
 fiber scenario with no fiber crossing is 
 shown .
the dt - based and odf - based image clustering 
 produce the same results . 

hence , odf clustering reproduces dt - based 
 results on a simple fiber population 
 example.finally , figure 4 
			 shows a fiber crossing 
 scenario with two overlapping fiber bundles 
 that have different geometries . 
 segmentation was performed over the dt 
 and the odf image shown in 
 figure 5 . 

moreover , the odf n - cuts 
 segmentation exhibits artifacts not 
 present in the odf diffusion maps segmentation . 

the odf diffusion maps 
 effectively identify the two different 
 fiber bundles as well as the fiber 
 crossing areas . in figure 3 , a single 
 fiber scenario with no fiber crossing is 
 shown .
the dt - based and odf - based image clustering 
 produce the same results . 

hence , odf clustering reproduces dt - based 
 results on a simple fiber population 
 example .
finally , figure 4 
			 shows a fiber crossing 
 scenario with two overlapping fiber bundles 
 that have different geometries . 

segmentation was performed over the dt 
 and the odf image shown in 
 figure 5 . 

moreover , the odf n - cuts 
 segmentation exhibits artifacts not 
 present in the odf diffusion maps segmentation . 

the odf diffusion maps 
 effectively identify the two different 
 fiber bundles as well as the fiber 
 crossing areas .
the real - data 
 experiment presented in this section shows the 
 segmentation and labeling of a 
 cropped axial and coronal slice .
the cropped 
 slices were chosen by an expert in 
 regions of known fiber crossings where the dt 
 model is normally limited .
hence , our segmentation 
 algorithm is confronted with elements that 
 have different orientation and 
 different diffusion characteristics . in order to show that odf data segments the white 
 matter fiber bundles better than the dt data in 
 real cases , we analyze the 
 evolution of the affinity matrix as the scale - space parameter changes in the 
 axial cropped slice shown in 
 figure 6 . 

affinity matrices were computed 
 with varying scale - space parameter between 1/5 , 1/10 , 1/20 , and 1/40 of the quantity 
 of elements ( |x| ) to cluster , respectively . in order to show the block 
 structure of the affinity matrices ,
it can be seen in figure 7 
 that as 
 the scale diminishes , the dt data shows a 
 high correlation between all the 
 elements of the slice .
this makes clustering 
 very difficult because the blocks 
 are small and highly correlated . 
 on the other hand ,
this block structure shows a high 
 correlation of the elements within 
 each block and a low interblock 
 correlation , giving a much better input to 
 the clustering algorithm than the dt 
 data . in figure 6 , the location of the 
 cropped axial slice is shown in the axial slice , figure 6(a ) , and coronal slice , figure 6(b ) .
as it can be seen in the segmented 
 and labeled axial slice , figure 8 , 
 the segmentation also allows to identify and label some of the 
 main white matter structures , corpus callosum ( cc ) , anterior 
 corona radiata ( acr ) , forceps major ( fmajor ) and forceps minor 
 ( fminor ) . in figure 9
, the location of the 
 cropped coronal slice is shown in the axial slice , figure 9(a ) , and coronal slice , figure 9(b ) . as it can be seen in the segmented 
 and labeled coronal slice , figure 9(c )
, the segmentation allows to identify and label main 
 white matter structures : corpus callosum ( cc ) , cingulum ( cg ) , 
 corona radiata ( cr ) , superior longitudinal fasciculus ( slf ) .
note 
 that the segmentation is resilient to crossing areas such as seen 
 at the interface between cr and cc .
we have presented 
 an algorithm to perform q - ball imaging segmentation of 
 white matter fiber 
 bundles .
the proposed method combines state - of - the - art 
 hardi reconstruction and 
 state - of - the - art spectral clustering techniques . 

our algorithm is 
 initialization - free and has only two parameters . 
 a scale - space parameter and 
 the number of regions ( clusters ) are to be found . 
 regarding this number of clusters 
 parameter , we have proposed to estimate it automatically . 

we have introduced a 
 spectral embedding technique that does not 
 require uniform sampling of the 
 elements .
to do so , the affinity measure used 
 incorporates an euclidean distance 
 measure between the spherical harmonic 
 coefficients describing the q - ball odfs 
 and also incorporates the spatial location distance 
 between odfs . the affinity 
 measure and
the metric induced in the embedded 
 space is then used to cluster 
 q - ball odf images into multilabel segmentation 
 representing the fiber bundles . 

spectral embedding has already been applied to 
 dmri ( e.g. , ) . 
 however , to our knowledge , this is the first work 
 using the diffusion maps that avoids the high 
 dependence on element sampling . 

we have illustrated that the odfs are the desirable 
 elements to use for clustering in the white 
 matter because the classical dt 
 model is limited in regions of fiber crossings . 

the odf is even more attractive 
 because of the recent analytical spherical 
 harmonic solution to the odf 
 reconstruction [ 3 , 8 , 9 , 13 ] .
the analytical solution is in 
 fact as fast as a 
 standard dt least - square estimation . in this work
, 
 we believe that we have used 
 the state - of - the - art odf reconstruction method , which is regularized , 
 robust and very simple to implement .
the spectral embedding performed by the diffusion maps 
 technique is at the heart of our segmentation 
 algorithm . whereas other spectral 
 embedding techniques have a tendency to produce 
 artifacts in the presence of 
 different sampling characteristics within a 
 cluster , the technique used in this 
 work greatly reduces this tendency by 
 performing the simple linear algebra 
 calculation shown in ( 17 ) .
spectral embedding techniques produce a representation 
 of the embedded data based on element - to - element affinities . 

this leads to the 
 fundamental issue : how to choose the affinity measure ? 

it is a challenge to 
 find a measure that incorporates similarities 
 between elements as well as the 
 spatial location difference between elements . 
 for similarities between 
 elements , we chose the euclidean distance between 
 spherical harmonic 
 coefficients describing the odfs . 

this approach is simple and very efficient 
 because it allows to process the odfs 
 directly on the sh coefficients .
the 
 euclidean distance has also been used 
 successfully in a level set segmentation 
 framework  
 and it would be interesting to compare our spectral clustering 
 approach against it . for spatial location 
 difference , we chose markovian 
 relaxation in order to be consistent 
 with the graph theoretical representation 
 of the diffusion maps technique . 

although this way of representing the distance 
 involves an artificial elimination of 
 all the nonneighboring relations of the 
 odf elements in the affinity matrix and an 
 adjustment of the diagonal elements , 
 we believe that the resulting affinity 
 relations represent the affinity 
 better .
the affinity of two neighboring 
 elements at the beginning of the 
 markovian relaxation algorithm is 
 represented by a function of the euclidean 
 distance between them .
this affinity 
 can be interpreted as the probability that 
 a random walker has of going from the 
 first element to the second .
the affinity 
 of two elements at the end of the 
 relaxation is the probability of a random 
 walker starting from one element and 
 reaching the second in a certain number of 
 steps .
we believe that 
 there is room for improvement in this last 
 part of the algorithm . in the first 
 place , the k - means algorithm needs an 
 explicit number of clusters to find .
this 
 can be heuristically determined by analyzing 
 the eigenvalue plot or the 
 reordered affinity matrix structure , as shown in 
 this work .
however , an 
 automatic method that could find the 
 number of clusters would considerably 
 improve the algorithm . in the second place , 
 the k - means algorithm and its 
 variants , for instance , k - medians , k - medioids , 
 search for isotropic clusters in 
 the embedding space . 

this could also improve 
 the last clustering phase of our 
 algorithm finally , in order to analyze the importance of the 
 difference between our diffusion maps 
 algorithm and the widely used n - cuts , we 
 used synthetic simulations . 
 in these simulations , we generated a synthetic 
 image with a single cluster within which the 
 sampling of the elements changed . 

we showed that when this sampling changes , the 
 n - cuts algorithm produces 
 artifacts while our diffusion maps method 
 does not . as uniform sampling within 
 a cluster is a difficult property to guarantee in the 
 white matter fiber 
 bundles , our diffusion maps method is 
 better suited for this task .
we thus 
 believe that diffusion maps are the right clustering method to be used on dmri 
 processing problems .
, 
 we have shown that in order to 
 perform spectral clustering on complex 
 dmri with crossing fiber bundles , an 
 hardi technique such as q - ball imaging is better 
 than the classical dti 
 technique .
this is because the odf reconstructed from 
 qbi is able to recover 
 multiple crossing fiber populations . 

second , a diffusion maps - based technique 
 for image segmentation was introduced to 
 reduce artifacts arising from the 
 widely used n - cuts image segmentation . 

we have illustrated the advantages of 
 the odf diffusion maps segmentation algorithm , 
 and showed on a real dataset 
 that our algorithm is able to identify important 
 and complex white matter fiber 
 bundles .
finally , the diffusion maps technique has been shown 
 to be more robust to sampling frequency 
 variations within each object to be 
 segmented . in order to cluster the elements , 

we have used an adaptive 
 scale - space parameter and we have used 
 markovian relaxation in order to 
 incorporate spatial dependencies . 

overall , the approach is theoretically sound 
 with the graph - based representation which 
 lies at the heart of spectral 
 clustering methods
. therefore , we have an algorithm to perform fiber 
 bundle clustering for a single brain . 

it is now important to study the behavior 
 over several subjects in order to assess 
 the reproducibility of the algorithm . 
 in time
, this will enable to perform 
 multisubject statistics within bundles in 
 the embedded space .
this will help characterize 
 the white matter fiber bundles 
 of several subjects and study if the alteration of these segmented tracts can 
 provide new biomarkers for white matter diseases .
 white matter fiber clustering aims to get insight about anatomical structures in order to generate 
 atlases , perform clear visualizations , and compute statistics across subjects , all important and current 
 neuroimaging problems . in this work 
 , we present a diffusion maps clustering method applied to diffusion 
 mri in order to segment complex white matter fiber bundles . 
 it is well known that diffusion 
 tensor imaging ( dti ) is restricted in complex fiber regions with crossings and this is why recent 
 high - angular resolution diffusion imaging ( hardi ) such as q - ball imaging ( qbi ) has been introduced 
 to overcome these limitations . 
 qbi reconstructs the diffusion orientation distribution function ( odf ) , 
 a spherical function that has its maxima agreeing with the underlying fiber populations . in this paper , 
 
 we use a spherical harmonic odf representation as input to the diffusion maps clustering method . 
 we 
 first show the advantage of using diffusion maps clustering over classical methods such as n - cuts and 
 laplacian eigenmaps . in particular , our odf diffusion maps 
 requires a smaller number of hypothesis 
 from the input data , reduces the number of artifacts in the segmentation , and automatically exhibits the 
 number of clusters segmenting the q - ball image by using an adaptive scale - space parameter . 
 we also 
 show that our odf diffusion maps clustering can reproduce published results using the diffusion tensor 
 ( dt ) clustering with n - cuts on simple synthetic images without crossings . on more complex data with 
 crossings , 
 we show that our odf - based method succeeds to separate fiber bundles and crossing regions 
 whereas the dt - based methods generate artifacts and exhibit wrong number of clusters . 
 finally , we show 
 results on a real - brain dataset where we segment well - known fiber bundles . 
1. INTRODUCTION
2. METHODS
3. RESULTS
4. DISCUSSION
5. CONCLUSIONS