Abstract:
A computer-implemented method to correct motion and interpolation effects for functional magnetic resonance imaging (fMRI) analysis is provided. The method estimates the motion on every voxel of the data and removes those effects to leave a residual signal that can be analyzed to high accuracy. The estimation of the motion includes solving a normal matrix equation based on the local translational motion of each voxel of the head, and a regularization parameter that depends on the local spatial structure of the head. Methods to regularize a matrix from the normal equation using the regularization parameter are also provided. A rolling filter implementation for real-time processing and motion correction is provided.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
     This application claims priority from U.S. Provisional Patent Application 61/000,140 filed Oct. 23, 2007, which is incorporated herein by reference. 
    
    
     STATEMENT OF GOVERNMENT SPONSORED SUPPORT 
     This invention was made with Government support under contract K25 MH077309 awarded by National Institute of Mental Health. The US Government has certain rights in the invention. 
    
    
     FIELD OF THE INVENTION 
     The invention relates generally to functional magnetic resonance imaging (fMRI). More particularly, the present invention relates to motion correction for fMRI. 
     BACKGROUND 
     Clinical subjects often have large and rapid movements during a functional magnetic resonance imaging (fMRI) scan session, and both of these types of errors need to be avoided or corrected for a successful analysis. Many current clinical data sets cannot be analyzed because of subject motion during the scan session. Head motion causes failures in fMRI analysis and a great deal of effort has been expended to train subjects to be motionless, or constrain their heads in the scanner via packed pillows, bite bars, or custom head cast. These physical constraints are unpleasant at best, and not possible to use on many young, old, or disturbed subjects of interest. 
     Methods have been developed to attempt to correct for subject motion without solely relying on physical constraints. The first step in fMRI processing is to realign the data so all the images are registered to each other. Even after realignment and interpolation, there is unavoidable residual time series errors caused by the subject&#39;s movements. These errors occur regardless of the interpolation algorithm that is used because the spatial high frequency structure of the brain cannot be adequately captured with any interpolation algorithm using data samples that are spatially far apart. In other words, the errors are inherent due to limited spatial resolution of the data samples compared to brain structures of interest. Techniques using motion regressors in the design matrix during fMRI analysis have been developed to correct for motion correction. However, these existing techniques have not been effective against typical subject motions in high motion clinical data sets. In fact, those techniques are completely unusable for many high motion data sets. Technical papers relating to these techniques include Friston et al., “Movement related effects in fMRI time-series,” Magnetic Resonance in Medicine, Vol. 35 (1996), pp. 346-355 and Grootoonk et al., “Characterization and Correction of Interpolation Effects in the Realignment of fMRI Time Series”, NeuroImage 11, 49-57 (2000). Grootoonk et al. utilizes a principle component analysis (PCA) that approximates the variations in intensity due to subject motion by using two functions: sine and cosine. However, the method of Grootoonk et al. fails to transition smoothly in the limit of low subject motion. This failure is unavoidable due to the functions that naturally arise as a consequence of the PCA in the algorithm of Grootoonk et al. Furthermore, the method of Grootoonk et al. is based on faulty assumptions by addressing the issue of poor interpolation algorithms and fails to appreciate the issue that actual fine-grain brain structure makes even perfect interpolation algorithms fail. 
     Existing techniques for motion correction in fMRI also suffer from difficulties with regularization. In particular, many existing techniques are over-regularized, causing errors for the corrected motion. Current fMRI motion correction algorithms also rely on poorly aligned images, thereby introducing unreliable results. In addition, many existing algorithms require many reads of the images to be corrected, which is computationally and time inefficient. 
     SUMMARY OF THE INVENTION 
     The present invention is directed to a method of correcting motion for functional magnetic imaging (fMRI). The method includes having a plurality of aligned images, where the voxels of each image corresponds to the voxels of another image. The alignment provides one or more translational distances between the voxels in each image with the corresponding voxels in another image. Alternatively or additionally, the method can include calculating the translational distances for each voxel. A time-series vector is defined for each voxel to represent the value of the corresponding vectors over a number N of images. A motion correction estimator is calculated for each voxel based on the time-series vector and a plurality of correction parameters, such as nonlinear parameters. The correction parameters are based on the translational distances, are approximately periodic in the size of the voxels, and approach zero when the translational distances approach zero. Each voxel is corrected based on the motion correction estimator. 
     In an embodiment of the present invention, one or more of the aligned images are removed and excluded in the calculation of the motion correction estimator. Images can be removed based on a shift of one or more voxels from one image to another. In particular, an image can be removed when the shift is greater than a threshold shift. In an embodiment, the threshold shift is iteratively changed based on the number or percentage of images to be removed. 
     In an embodiment of the present invention, the translation distances used in calculating the motion correction estimator approximates the translational and rotational motion of the whole head and are in one more translational directions. In a preferred embodiment, the calculation of the motion correction estimator for one of the voxels includes solving a normal equation having one or more matrices. Optionally, at least one of the matrices of the normal equation is regularized by setting a minimum value for one or more of the diagonal elements of the matrix. In a preferred embodiment, the regularization is based on one or more spatial gradients or curvatures of the images. In particular, regions of an image having higher spatial gradient have a lower minimum regularization value than regions of an image having lower spatial gradient. 
     A rolling filter implementation is also provided, wherein each of the images are read at most twice. In the first reading, the motion correction estimator is calculated for all of the voxels and the corrected images are written during the second reading. The rolling filter implementation facilitates motion correction for partial data sets. 
    
    
     
       BRIEF DESCRIPTION OF THE FIGURES 
       The present invention together with its objectives and advantages will be understood by reading the following description in conjunction with the drawings, in which: 
         FIG. 1  shows a flow chart of an example method for correcting motion for fMRI images according to the present invention. 
         FIG. 2  shows an example set of fMRI images with subject motion according to the present invention. 
         FIG. 3  shows matrices and equations used in an embodiment of the present invention. 
         FIG. 4  shows a flow chart of a rolling filter embodiment according to the present invention. 
     
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     Large and rapid movements by a subject during a functional magnetic resonance imaging (fMRI) scan session cause errors that need to be corrected for a successful analysis. The first step in fMRI processing is to align the data so all the images are registered to each other, However, even after alignment and interpolation, there are unavoidable residual time series caused by subject motion. It is important to note that these errors occur regardless of the interpolation algorithm that is used because the spatial high frequency structure of the brain cannot be adequately captured with any interpolation algorithm. The present invention is directed to a robust method for correcting residual effects from large motions. 
     An fMRI analysis utilizes a set of successive images, with all images having a large number of voxels. A typical fMRI image has approximately 80,000 voxels. The fMRI images are generally aligned after they are captured. Alignment includes corresponding the voxels of one image with another image to derive one or more translational distances and/or rotations between corresponding voxels of multiple images. Typically, the translations and rotations are provided with respect to a reference image. 
       FIG. 1  shows a flow chart of a motion correction method of the present invention. In the embodiment shown in  FIG. 1 , the method requires alignment of the raw images that are captured. In a preferred embodiment, the motion correction method of the present invention uses an algorithm that takes as input the alignment parameters and the aligned, sliced images. In other words, the method can be applied to images that have already been realigned. In either case, the algorithm estimates and removes residual motion effects for each voxel, and writes a new set of images where the motion effects have been suppressed 
       FIG. 1  shows an optional step of removing one or more images for which the alignment parameters are likely to be incorrect. This step is important for clinical subjects because the large rapid motions of clinical subjects can cause errors in the alignment calculation and the errors are likely to propagate to the estimate of motion correction. In an embodiment, images are removed based on the shift or root mean square (RMS) deviation of one or more voxels from one image to another image. In particular, the determining characteristic is whether voxels near the outside of the brain shifts by more than a threshold distance from one image to the next or another image. In a preferred embodiment, the RMS of the motions of six voxels at 65 mm distance from the origin is computed and compared with a threshold shift of 0.5 mm, wherein if the RMS is greater than the threshold shift, the image is removed. It is noted that any number and selection of voxels, and any threshold shift can be used to determine removal of images. In an embodiment, the shift threshold is iteratively changed based on a number or percentage of images to be removed. For example, if more than 25% of the images are recommended for exclusion, the specified distance threshold is increased until fewer than 25% of the images will be excluded. 
     After removal of the poorly aligned images, N images remain for use in the motion correction algorithm. For each voxel, a time series vector u is defined, where the elements of u represent the value of the corresponding voxels over the N remaining images. In an embodiment the time series vector u is defined with respect to a baseline value u o . A motion correction estimator g is then calculated for each voxel, where the calculation is based on the time series vector u of the same voxel and one or more correction parameters. 
     It is important to note that the correction parameters used to calculate g are approximately periodic in voxel size and approach zero in the zero motion limit. Using correction parameters with these properties allows for accurate motion correction that extrapolates correctly in the limit of small motion. In a preferred embodiment, the values of the voxel in every image, including the removed images, are then corrected based on the motion correction estimator g. New images are then written with the corrected voxels. 
       FIG. 2  shows an example of multiple aligned fMRI images  210 - 230 . The alignment corresponds the voxels of one image to another image. For example, corresponding voxel  200  is highlighted in images  210 - 230 .  FIG. 2  also shows example plots displaying residual motion in the images  210 - 230 . In an embodiment of the present invention, the translations and rotations of the whole head are locally approximated as translational distances on the voxel. For example, x 2 ′ is the approximate translational distance (in mm) in the x direction for the second image and y 2 ′ is the approximate translational distance in the y direction. Though  FIG. 2  shows only two-dimensional plots, embodiments of the present invention are preferably directed to, but not limited to, three-dimensional motion. 
     As mentioned above, a time series vector u is defined for each voxel, where u is an N by 1 vector as shown in  FIG. 3 . In a preferred embodiment, the motion correction estimator g is a 7 by 1 vector and is the solution to the equation: u=Mg, where M is an N by 7 matrix with i th  row: M i =[sin(x i ) 1−cos(x i ) sin(y i ) 1−cos(y i ) sin(z i ) 1−cos(z i ) 1], and as shown in  FIG. 3 . The elements of matrix M include the correction parameters. In the embodiment shown in  FIG. 3 , the elements of matrix M include correction parameters that are functions of the translational distances. For example, x N  in elements  310  are functions of x N ′ by the following transformation: x N =2πx N ′/v x , where v x  is the voxel size in the x direction. In other words, x N  is the distance in the x direction for the N th  image relative to the voxel size. Translations in the y and z directions are expressed similarly. It is noted that x i , y i , and z i  can be different for every image and different for every voxel. 
     A “normal” equation 320, g=(M T M) −1 M T u, used to calculate the motion correction estimator g in an embodiment of the present invention. In a preferred embodiment, the matrix M T M is regularized before inversion. In one embodiment, regularizing includes setting a minimum value for one or more of the diagonal elements of matrix M T M The minimum value must be small so that important gray matter voxels near the surface of the brain can be adequately corrected even for large motions, yet it cannot be too small, or the solution will become unstable. For example, regularization values within 0.001 and 0.01 give acceptable results with 64-bit arithmetic. In an embodiment, all of the diagonal elements of matrix M T M are set to have a minimum value of about 0.005. Alternate embodiments could use other minimum values or alternate regularizing schemes, such as adding a small value to the diagonal elements of M T M, can be employed. 
     In an embodiment of the present invention, one or more matrices used in calculating the motion correction estimator is regularized based on spatial gradients of the images, spatial curvatures of the images, one or more alignment parameters, or any combination thereof. For example, an image has multiple regions, one of which has a higher spatial gradient or curvature than another region. The minimum value used in regularizing is set to be lower for the high gradient region than for the low gradient region. In an embodiment, the minimum value ranges from 0.005 to 0.05 for a high gradient region and is greater than 1 for a low gradient region. 
     In an embodiment, the motion correction estimator g can be stored as seven three-dimensional images for each voxel. The stored images are related to the seven elements of the vector g and can be stored for user review. 
     The value of each voxel is corrected based on the motion correction estimator g associated with the same voxel. In a preferred embodiment, the motion corrected time series v of a voxel is an N by 1 vector and is calculated from the equation: v=u−Mg+u o , where u, M, g, and u o  are as described in the preceding paragraphs. It is noted that in a preferred embodiment, the motion correction is applied to all images, including those excluded from use during the calculation of the motion correction estimator g. The corrected images can be written for use by succeeding stages of fMRI processing, such as artifact repair, normalization, smoothing, and estimation. 
     Motion correction generally requires substantial amount of time and computational power since the calculations described above are performed for every voxel and there are typically about 80,000 voxels per image. To alleviate this large time-consumption, an embodiment of the present invention is directed to a rolling filter implementation where the images are only read at most twice, as shown in  FIG. 4 . In the embodiment of  FIG. 4 , all of the voxels in an image are processed in parallel or at approximately the same time. The first pass  410  through the images accumulates all of the terms needed to solve for g. Using this approach, it is only necessary to save the 49 images of M T M and the 7 images of M T u to solve for g at every voxel simultaneously. Optionally, the matrices M T M can be regularized. The second pass  420  through the images applies the motion correction estimator g for every voxel and writes the motion corrected images. In this embodiment, motion correction takes less time than a typical realignment step in fMRI processing. In addition, the rolling filter implementation can be applied to partial data sets, hence, the approach could be used progressively from image to image as a basis to provide real-time image correction for fMRI. 
     As one of ordinary skill in the art will appreciate, various changes, substitutions, and alterations could be made or otherwise implemented without departing from the principles of the present invention, e.g. the present invention can be applied to any set of images of any dimension. Accordingly, the scope of the invention should be determined by the following claims and their legal equivalents.