diffusion mri has provided a great tool for neuroscientists to understand and analyze in vivo the anatomy of the brain white matter fiber tracts that connect different areas of the cortex .
the diffusion tensor model   has become increasingly popular , and the study of scalar indices derived from it has proved useful in the diagnosis of a wide range of neurological diseases [ 2 , 3 ] . for several specific applications , like fiber tractography , this model is , however , known to be limited , and high angular resolution imaging techniques should be used instead , to reconstruct the model - free ensemble average propagator [ 46 ] or the orientation distribution function ( odf ) [ 710 ] .
the acquisition of high angular resolution diffusion images requires longer time than diffusion tensor imaging .
subjects are likely to move during these acquisitions , and we can identify at least three motivations to develop a proper method for the online detection of motion .
first , images can be registered prior to diffusion model estimation ; however this might increase partial volume effects , because of the relatively low spatial resolution of diffusion - weighted images and of interpolation in the registration procedure .
this also modifies the variance properties of the image , which should be considered carefully for group studies .
when the subject moves during acquisition , a warning could be issued , so as to take a decision accordingly .
depending on the number of images already acquired , the decision could be to restart the scan or acquire a few more diffusion weighted images than originally planned to compensate for the variance increase due to the registration .
second , diffusion acquisitions use a gradient table , which is a set of orientations and b - values and has been designed following an optimal sampling strategy . in q - ball imaging for instance , the set of orientations is designed to sample the sphere in an optimal isotropic fashion [ 13 , 14 ] .
when correcting for motion , the diffusion encoding gradients should be rotated to be consistently defined in a coordinate frame related to the subject [ 1517 ] .
this modification might break the optimal sampling strategy as originally planned and affect the reconstruction of the odf .
finally , in the context of online processing of diffusion images , motion must be detected , so that it can be corrected to continue the incremental reconstruction .  
the authors in   use a camera inside the scanner to detect and evaluate a rigid motion .
their study shows the improvement in odf reconstruction with this prospective approach for motion correction over a classical offline registration .
however , this technique requires additional hardware which is to date not always available on scanners .
other approaches [ 19 , 21 ] are based on the interleaving of echo navigators through the acquisition sequence .
the authors in [ 19 , 21 ] report good results in detecting and correcting motion , but these additional acquisitions affect the overall protocol time .
finally , a recent work   introduces a motion detection and prospective motion correction to account for slow motion artifacts such as image misalignment .
they also reduce fast motion effects such as signal dropout , by identifying the volumes most affected by motion , and schedule reacquisition at the end of the scan .
but the motion detection and correction is performed by comparing a diffusion weighted image to the average of diffusion weighted images with the same b - value .
we believe this is suboptimal , as it does not take into account the direction associated to each diffusion weighted image .
, we propose a diffusion weighted - image - based technique for the online detection of motion in q - ball imaging .
our method does not require new hardware or change in the acquisition protocol and is based on a kalman filter reconstruction of the hardi signal [ 22 , 23 ] .
the first contribution is the adaptation of the kalman filtering framework for online reconstruction of odf in constant solid angle recently introduced in the q - ball imaging community [ 9 , 10 ] .
then , we present a solution to the detection of motion in diffusion images , adapted from the generalized likelihood ratio test ( glrt ) . to overcome certain shortcomings of this method , we introduce star ( statistical analysis of residuals ) , an original method for the detection of motion in diffusion weighted images .
the method is tested under various experimental conditions on semiartificial and on real data and compared to glrt . in the results ,
we report successful detection of small motion ( rotation by angle under 2 ) , even in noisy conditions .
the detection using star outperforms glrt , while star does not need any delay for the detection .
in this section , we review the definition and the expression of the odf calculated in constant solid angle .
it has been shown recently that this mathematically correct definition of odf can be reconstructed in q - ball imaging [ 9 , 10 ] .
we present an incremental reconstruction algorithm for this odf , based on the kalman filter .
we formalize the problem of motion detection only from the observation of the diffusion signal .
we present a brief review of methods for fault detection , in particular glrt , built upon the kalman filter , as first described by . finally , we present star , an original approach based on a statistical modeling of the image .
both algorithms are implemented , and we present at the end of this section the validation methods used to compare them .
the odf is a spherical function , retaining the angular information of the ensemble average propagator , p. when defined as the marginal probability of direction , the odf ,  , is the probability for a water molecule to diffuse along a given direction in a constant solid angle .
, we measure the signal , s(q ) , which is related to the ensemble average propagator p through a fourier transform , under the narrow - pulse condition 
 ( 2)p(r)=q3s(q)s(0)e2iqrdq . under the assumption of a monoexponential decay of the diffusion signal
s , the relation between s(q ) , s(0 ) , and the odf  is given by 
 ( 3)(u)=14+1162frt{b2ln(lnss(0))}(u ) , 
 where frt denotes the funk - radon transform and b the laplace - beltrami operator .
the computation of the odf can be implemented using the modified spherical harmonic basis for real and symmetric functions   to describe the transformed signal y = ln(ln(s / s(0 ) ) ) , as both the funk - radon transform and the laplace - beltrami operations in ( 3 ) have a close - form matrix expression in the spherical harmonic basis
. if c^j are the spherical harmonic coefficients that describe y , then the spherical harmonic coefficients to describe the odf  are 
 ( 4)c^j={12j=118(1)lj/213(lj+1)24(lj2)c^jj>1 , 
 where lj = { 0,2 , 2,2 , 2,2 , 4,4 , 4 ,  } for j = { 1,2 , 3 ,  } is the order associated to the jth spherical harmonic .
the computation of the spherical harmonic coefficients c^ describing y from a series of measurements y[k ] = ln(ln(s[k]/s(0 ) ) ) = ln(ln(s(q u[k])/s(0 ) ) ) , k = 1  n at discrete positions u[k ] on the unit sphere , and a measurement without any diffusion encoding gradient s(0 ) is implemented by minimizing : 
 ( 5)m(c)=(ybc)t1(ybc)+ctlc . 

the second term is a laplace - beltrami regularization constraint on the fitted signal , with the matrix l having diagonal elements lj(lj+1 ) .
the matrix  in the data fitting term of ( 5 ) accounts for the uncertainty in the diffusion - weighted measurements s[k ] as well as for the distortion introduced by the nonlinear transform , which is illustrated in figure 1 . the distortion is higher in high - magnitude and low - magnitude signal modes .
the diagonal elements of  can be approximated through first - order error propagation ; the uncertainty on the transformed signal y is simply 
 ( 6)y=yss=1sln(s / s0)s . 

provided that the error on separate measurements is uncorrelated , the diagonal elements [k ] of  are simply 
 ( 7)2[k]=var(s[k])s[k]2ln2(s[k]/s0 ) , 
 where var(s[k ] ) denotes the variance of the diffusion signal s[k ] and can be estimated once for the whole volume using a method like piesno for instance .
provided that the acquisition sequence is incremental ( in this study we use the incremental point sets as in ) , the energy in ( 5 ) can be minimized incrementally using a kalman filter .
the incremental system adapted to the reconstruction of the odf in constant solid angle is given by 
 ( 8)initialization{c=[c]p~=[(cc)(cc)t]p=(p~1+l)1,update{v[k]=b[k]p[k1]b[k]t+2[k]g[k]=p[k1]b[k]tv[k]1p[k]=(ig[k]b[k])p[k1][k]=y[k]b[k]c[k1]c[k]=c[k1]+g[k][k ] .
the [k ] depend on the data as expressed in ( 7 ) , and the covariance v[k ] of the residual (k ) can no longer be precomputed offline .
the expected covariance of the estimated spherical harmonic coefficients c[k ] is the matrix p[k ] computed by the kalman filter .
then the expected mean squared error on the spherical harmonic coefficients describing the odf is given by 
 ( 9)mse(c[k])=tr(ftltp[k]lf ) , 
 where f is the matrix form of the funk - radon transform and has diagonal elements 2plj(0 ) , where plj(0 ) is the legendre polynomial of degree lj evaluated at 0 , and l is the laplace - beltrami matrix as in ( 5 ) .
the kalman filter was derived with the assumption that the local diffusion propagator does not change in time .
next , we show how we can derive a motion detection algorithm from this kalman filter .
the detection of abrupt changes in dynamical systems has been extensively studied [ 24 , 27 ] ; a very good review of methods and algorithms can be found in .
they propose a classification of change detection problems , together with suggested methods and algorithms to address them . in the previous section ,
we have introduced a kalman filter solution to reconstruct the spherical harmonic coefficients of the q - ball signal .
the state of our system is the vector of spherical harmonic coefficients c[k ] , and a motion of the subject at time  is likely to imply a modification of this state , c[k   ] = c[k <  ] + p. the problem of motion detection reduces to the problem of change detection in this multidimensional system . besides , since both the time  and the magnitude p of the change are unknown a priori , the classification in   suggests to use a generalized likelihood ratio test ( glrt ) for the detection . in the next section
we briefly describe this method and its implementation upon a kalman filter , as originally introduced in .
the kalman filter presented in the first section is built under the hypothesis of no motion .
we can monitor the residuals of this kalman filter for each iteration and test whether the hypothesis is still valid . as it has been shown in ,
the prediction error after a change occurred at time  for subsequent iterations can be decomposed as 
 ( 10)[k]=g(k,)p+1[k ] , 
 where 1 is zero - mean gaussian distributed with covariance v[k ] and g(k ,  ) represents the propagation of a jump at time  , to the prediction error at time k. this can be computed as in : 
 ( 11)g(k,)=b[k](ij=k1g[j]g(j,)),g(k , k)=b[k ] .
the problem of a change detection is to discriminate between two hypotheses : ( 0 ) : no change in the state vector : [j ] = 1[j ] , j = 0  k,(1 ) : at time 0 , the state vector becomes c + p0 .
when p0 and 0 are known , a natural statistic for the detection is the likelihood ratio .
provided that the residuals are linearly related to the change ( 10 ) , the log - likelihood ratio is 
 ( 12)l(k;0,p0)=lnp1([0k])p0([0k ] ) . 

provided that the densities p0 and p1 are gaussian , after simplification this is rewriten as 
 ( 13)l(k;0,p0)=j=0k[j]v1[j]g(j,0)p0 .   no change in the state vector : [j ] = 1[j ] , j = 0  k ,   at time 0 , the state vector becomes c + p0 .   in our case , both p and  are unknown .
the generalisation of the likelihood ratio method suggests to replace 0 and p0 in ( 13 ) by their maximum likelihood estimates : 
 ( 14)^(k)=argmax  l(k;,p^(k;)),p^(k;)=(j=kgt(j,)v[j]1g(j,))1  j=kgt(j,)v[j]1[j ]  ( least  squares  estimate ) . 

finally , the decision is taken by comparing l(k;^,p^ ) to a threshold . this technique works fine , yet suffers from several drawbacks .
first , the calculation of the maximum likelihood estimate of p involves the inversion of a matrix in ( 14 ) which has full rank only if k   > dim(p ) . in other words
, this implies a delay at least equal to the dimension of the problem . as an example ,
when the signal is fitted in the 4th order symmetric spherical harmonic basis , this dimension is 15 .
in addition , the choice of a threshold  was reported to be critical and highly dependent on the delay .
finally , in our situation the state vector represents the diffusion signal locally , and glrt does not say how to combine the statistics of different voxels to calculate a statistic which could be an indicator of motion for the whole volume at once . to address these weaknesses ,
we propose in the next section an original approach without delay , incorporating a statistical model of the image , in order to provide a more suitable detection algorithm .
the reconstructed image is a vector field c(r ) , where c is a vector of spherical harmonic coefficients describing the diffusion signal at voxel position r. we consider the difference p between two such vector fields c1 and c2 , representing the same subject before and after a rigid transform . in what follows ,
we consider p(r ) as a random variable , with unknown covariance matrix c. hence if there is no motion , the residuals for the whole volume will be distributed as (0 , v[k ] ) , where v[k ] is the variance predicted by the kalman filter .
otherwise the overall variance of the residuals will increase as v[k ] + g(k , )cg(k, ) , where g(k ,  ) is the matrix for the propagation of a jump at time  to the prediction error at time k , and the covariance matrix c of p is unknown . based on the previous observations , we design a test for motion detection without delay .
this means that based on measurements up to time k , we are able to decide whether a motion occurred at time k or not
. given a sample of m residuals at time k , at voxel positions x1  xm selected randomly within the brain , the hypotheses that a motion did occur at time  or not are equivalent to  ( 0 ) : [k ] has variance v[k ] , as predicated by the kalman filter ;  ( 1 ) : [k ] has a variance  > v[k ] .
: [k ] has variance v[k ] , as predicated by the kalman filter ; : [k ] has a variance  > v[k ] .
we first calculate the statistic : 
 ( 15)t=j=1m2(xj)[k](1/mj=1m[k])2v[k ] . 
 under the hypothesis ( 0 )
we want to reject the hypothesis with a significance level  : under the hypothesis ( 0 ) , we compute p such that (t > p ) = . the value of p is obtained from the inverse cumulative function of the m1 distribution .  
we implemented the incremental reconstruction using kalman filtering , together with glrt and star for motion detection .
these techniques were tested on real data , and a quantitative analysis of both was performed on semiartificial data , where the motion is simulated by a rigid transform . in this section , we describe how these images were synthesized . the simulation is based on a tensor field reconstructed images of still subject , acquired on a 3 t siemens magnet at the center for magnetic resonance research , university of minnesota medical school , with 200 encoding directions computed following the optimal sampling scheme of , b = 1000  s / mm , isotropic resolution 2.0  2.0  2.0  mm , 25 b = 0 images , 128  128 image matrix , 64 slices , te = 90  ms , and tr = 8500  ms .
we choose a series of diffusion gradient directions { g[k ] , k = 1  n } and a b - value for synthesis .
the rigid motion is specified by an instant  , its rotation component r and its translation vector t. provided that the diffusion encoding gradients should be rotated accordingly [ 1517 ] , the gradient directions used for synthesis are { g , g ,  , g[  1 ] , rg[ ] ,  , rg[n]}. the rigid transform is finally applied to the synthetic diffusion weighted images  ,  , n , after which the images are corrupted by rician noise .
we evaluate the general likelihood ratio , and the residual - based statistics computed for star as a motion detection criterion .
then we compare the sensitivity and specificity of glrt and star , for different values of the experimental parameter . within this section ,
we report the true positive rate ( tpr ) , defined as tpr = # detected positives/#positives , and the false positive rate , define as fpr = # mislabeled negatives/#negatives .  
the kalman filter and star were implemented in python , with the use of the scipy   toolkit , which is an efficient library for scientific computing .
based on this implementation , the reconstruction of the odf field for an image of dimensions 128  128  64 , with 200 diffusion - weighted images , took approximately 29  s on a 4-core intel core computer at 3.20  ghz , with 4.0  gb memory , running linux mint 13 .
this means that each diffusion volume is processed within less than 150  ms , which is short with respect to tr .  
one of the advantages of star outlined in the previous section is that the threshold for the detection can be deduced from the target false positive rate . in practice ,
as m becomes large , we approximate the m1 distribution for the decision test described in section 2.3.2 by a normal distribution : ( t - m+1)/2(m-1)~(0,1 ) , and p is given by the inverse normal cumulative density function . for a false positive rate fewer than 5% ,
. this value is experimentally tested , and the results are presented in the next section .
we report in figure 3 the value of the statistics ( t - m+1)/2(m-1 ) , for a series of 100 experiments without motion and a series of 100 experiments where the volume was rotated after 18 acquisitions by an angle of 2. the threshold was taken as 1.64 , for which the false positive probability is 5% .
the empirical false positive rate we report for these 200 simulations is 4% , while the true positive rate is 90% .
we computed both detection criteria on a series of 100 datasets without motion , and a series of 100 experiments with motion .
we plot in figure 4 the curve tpr versus fpr obtained by choosing different threshold values .
we also evaluate the robustness of glrt and star in various experimental conditions . for a fixed fpr
the experimental conditions include the delay , instant of motion in the acquisition sequence , snr , and motion magnitude .
the results of these simulations are reported in figure 5 . the experimental conditions , unless explicitly modified , were a rotation around the left - right axis by an angle of 3 , snr = 20 , motion instant  = 20 and a delay k   = 3 ( for glrt ) .
the experiment includes 400 negatives ( simulations without motion ) and 100 positives ( simulations with motion ) .
the monitoring of residuals in star and in glrt is limited to 500 voxels randomly selected within the brain to get a computational cost adequate for online treatment .
we also validate our methods on a real dataset , with the same imaging parameters as above . during the acquisition ,
the motion was a posteriori identified as a rotation of 20 about the z - axis ( see figure 6 ) .
the detection algorithms could detect this motion : with a delay of 2 acquisitions for glrt , and with no delay for star .
among the challenges of a motion detection algorithm , we have tested both glrt and star in these conditions : small delay for the detection ( ideally no delay);motion that occurred in the first few iterations of the kalman filter;very noisy conditions ( snr down to 10);small motion.with the help of simulated motion , we were able to present this exhaustive set of tests and subsequently the quantitative results .
in addition to these tests on artificially added motion , we also showed the application of the method in a real setting , with a healthy subject asked to move on purpose during the scan .
small delay for the detection ( ideally no delay ) ; motion that occurred in the first few iterations of the kalman filter ; very noisy conditions ( snr down to 10 ) ; both criteria show good results in detecting motion , even in severe experimental conditions . as expected , star is more robust to noise ( figure 5(a ) ) and performs better in detecting small motion ( figure 5(b ) ) , since it combines natively the residuals from different voxels .
in addition , glrt can not be computed if the number of acquired signals is lower than the dimension of the model , which is 15 in the case of 4 order real , symmetric spherical harmonics .
this impacts the ability to detect motion occurring at the beginning of an acquisition sequence : they are detected by star , while glrt can not be computed ( see figure 5(d ) ) .
in addition , glrt needs a delay greater than 6 to get similar sensitivity as star ( figure 5(c ) ) .
therefore , star is an original method that is best adapted to the problem of detecting motion online from a set of diffusion weighted images .   as we have reported in the implementation section ,
the reconstruction of the odf field using kalman filtering is fast and compatible with online implementation .
in this paper , we have proposed a method for the detection of motion in diffusion mri .
we have developed a kalman filter solution for the estimation of the odf in constant solid angle .
the detection algorithm star is based on the analysis of the residuals of the kalman filter , yet it is general and can be directly applied to any linear diffusion model reconstruction . compared to other techniques for the prospective detection and correction of motion [ 18 , 19 ] , our method does not require any camera or additional device .
once motion is detected by our technique , a decision could be taken by the scanner operator , or the protocol in   could be used for motion correction . to the best of our knowledge ,
our technique is the only image - based approach that clearly takes into account the directional information in diffusion weighted images to detect motion online .
the proposed technique was tested on semi - artificial data as well as in a real data and shows good results for the online detection of motion .
compared to glrt , which is a classical solution for the detection of changes in dynamical systems , star combines the residuals at different voxel positions to compute a statistic , on which the decision is based .
star performs better than glrt in the detection of small motion , motion in noise , or motion occurring early in the acquisition protocol . besides
, star does not need any delay for the detection , which makes it very efficient in practical situations .
 the acquisition of high angular resolution diffusion mri is particularly long and subject motion can become an issue . 
 the orientation distribution function ( odf ) can be reconstructed online incrementally from diffusion - weighted mri with a kalman filtering framework . 
 this online reconstruction provides real - time feedback throughout the acquisition process . in this article , 
 the kalman filter is first adapted to the reconstruction of the odf in constant solid angle . 
 then , a method called star ( statistical analysis of residuals ) is presented and applied to the online detection of motion in high angular resolution diffusion images . 
 compared to existing techniques , this method is image based and is built on top of a kalman filter . 
 therefore , it introduces no additional scan time and does not require additional hardware . 
 the performance of star is tested on simulated and real data and compared to the classical generalized likelihood ratio test . 
 successful detection of small motion is reported ( rotation under 2 ) with no delay and robustness to noise . 
1. Introduction
2. Methods
3. Results
4. Discussion
5. Conclusions