|Home | About | Journals | Submit | Contact Us | Français|
A major source of artifacts in diffusion-weighted imaging (DWI) is subject motion. Slow bulk subject motion causes misalignment of data when more than one average or diffusion gradient direction is acquired. Fast bulk subject motion can cause signal dropout artifacts in diffusion-weighted images and results in erroneous derived maps e.g. fractional anisotropy maps. To address both types of artifacts, a fully automatic method is presented that combines prospective motion correction with a reacquisition scheme. Motion correction is based on the Prospective Acquisition Correction (PACE) method modified to work with diffusion-weighted data (DW-PACE). The images to reacquire are determined automatically during the acquisition from the imaging data i.e. no extra reference scan, navigators or external devices are necessary. The number of reacquired images i.e. the additional scan duration can be adjusted freely. DW-PACE corrects slow bulk motion well and reduces misalignment artifacts like image blurring. Mean absolute residual values for translation and rotation were less than 0.6 mm and 0.5 degrees. Reacquisition of images affected by signal dropout artifacts results in diffusion maps and fiber tracking free of artifacts. The presented method allows the reduction of two types of common motion related artifacts at the cost of slightly increased acquisition time.
Diffusion-weighted magnetic resonance imaging (DWI) has become an established tool in clinical and research settings since a basic diffusion sequence was described by Stejskal and Tanner (1–3). Today, single-shot echo-planar imaging (EPI) is the method of choice for most applications unless very high spatial resolution or very small spatial distortions are required in which case multi-shot methods (4–6) are preferred. Diffusion imaging using EPI is sensitive to several artifacts (7). Some of these artifacts are scanner hardware related like distortions resulting from eddy currents and some artifacts are acquisition related like geometric distortions induced by magnetic susceptibility differences across brain tissue when using an EPI readout. Other artifacts are subject related like physiological motion (cardiac pulsation, respiration) and bulk motion. Left uncorrected, slow bulk motion during the acquisition leads to a mismatch of image data when multiple diffusion gradient directions or averages are acquired. This causes edge artifacts and blurring in the calculated diffusion maps, e.g. trace-weighted images or fractional anisotropy (FA) maps, as well as errors in the tensor estimation in the case of diffusion tensor imaging (DTI). In addition, fast bulk motion during the application of the diffusion gradients can lead to inhomogeneous signal attenuation artifacts (signal dropouts) in the diffusion-weighted images caused by additional phase terms (8). Non-rigid deformations of the brain caused by cerebrospinal fluid (CSF) pulsation lead to non-linear phase errors that also cause inhomogeneous signal attenuation artifacts (9).
Artifacts due to physiological motion are usually minor and can be mitigated by gating (10–12), however, bulk motion artifacts often result in unusable images when left uncorrected. While slow bulk motion can be corrected for retrospectively, real-time motion correction reduces motion-related effects on magnetization history and decreases the variance between volumes compared to retrospective motion correction. If the spatial resolution is anisotropic, as is often the case for clinical protocols, interpolation will also result in sub-optimal image quality. And in case of DTI processing it is necessary to take the apparent changes in diffusion gradient direction into account i.e. the b-matrix has to be corrected before further calculation can be done (13). Otherwise, erroneous tensor and derived values will be measured. Motion can also be corrected online with external motion tracking e.g. using optical methods (14–16). However, these methods require additional external hardware and software as well as time-consuming calibration steps before they can be used.
Images with signal attenuation artifacts cannot easily be corrected retrospectively. To prevent errors in the derived data, these images need to be excluded from processing, leading to a reduction in signal-to-noise ratio (SNR) of the derived maps and a bias in the tensor calculation. To address these problems, reacquisition methods can be used where segments or images that are affected by motion artifacts are reacquired during the same acquisition (6,17–19).
We propose a method that combines prospective motion correction with a reacquisition scheme. Prospective motion correction is used to correct artifacts due to slow bulk motion while reacquisition is used to correct artifacts due to fast bulk motion. Prospective motion correction and reacquisition are based on the image data; therefore they do not require any external hardware or software. The reacquisition scheme is based on the magnitude and phase data of the diffusion-weighted images and is fully automatic. Diffusion maps are calculated on the fly from uncorrected and corrected data.
A twice-refocussed DTI sequence (20) was modified to include prospective motion correction and reacquisition, completely implemented on the scanner giving corrected images and derived images immediately after acquisition. Prospective motion correction was based on the Prospective Acquisition Correction (PACE) algorithm, which uses a least squares cost function for image alignment (21), modified to use separate reference volumes for each b-value (DW-PACE). In short, each volume is compared to the corresponding reference volume with the same b-value and the estimated motion values - based on a rigid body model with six degrees of freedom: three translations and three rotations - are fed back to the scanner in real-time to adjust the acquisition frame of reference before the acquisition of the next volume. To clarify, the PACE algorithm does not try to predict motion but only applies the estimated motion values based on the current head position. The diffusion gradient frame of reference is adjusted at the same time. This ensures that the diffusion gradients are applied in a consistent frame of reference relative to the subject. The remaining motion between the acquired volume and the reference volume is corrected for online during the image reconstruction by resampling. The diffusion gradient vector is adjusted for the estimated motion as well before DTI processing is done. In comparison to the original PACE implementation for blood oxygen level dependent (BOLD) imaging where the reference volume is always the first volume of the acquisition, in our implementation, the reference volume for each b-value is continuously updated by averaging all previously acquired volumes with the same b-value in the current acquisition corrected for residual motion. This increases SNR of the reference volumes and diminishes the effect of diffusion anisotropy, which causes signal intensity changes in each voxel when diffusion gradients in different directions are applied.
Reacquisition was performed at the end of the normal scan time but within the same acquisition. The maximum number of reacquired images can be controlled freely as a percentage of the normal scan time or as number of images to reacquire. To determine which images need to be reacquired due to image artifacts, a score was calculated for each image in real-time during the image reconstruction in a two-step process: First, changes in the magnitude data for each slice from reference volume to currently acquired volume are detected and result in a score in the range of 1.0 to 2.0. Second, the phase data of any image that is considered not affected in the magnitude data is used to calculate a score with range 0.0 to 1.0 based on the slope of the phase data. In both cases, higher values correspond to a higher level of artifacts. The range of scores for magnitude data effects is higher compared to the phase data effects because motion has to be more severe to cause artifacts detectable in the magnitude data. The scores are fed back from the image reconstruction to the gradient control in real-time and the selected percentage of images with the highest scores are reacquired at the end of the normal acquisition. If there are fewer scores greater than zero than the selected number of images to reacquire, the scan will be terminated after all images with score greater then zero were reacquired. Therefore, for subjects without detected motion, acquisition time will not be increased compared to the regular acquisition time. Image reconstruction and calculation of diffusion maps (average non diffusion-weighted images, trace-weighted images, FA, color coded FA) are then performed on the scanner immediately after acquisition. This is done based on the original data as well as on a second set of data where the corrupted images have been replaced by the reacquired ones.
If the number of pixels r(s,v) for slice s and volume v with signal intensities above threshold t(v) of the current image is smaller than f1 times the number of pixels above threshold of the same slice of the first non diffusion-weighted volume r(s,1) it is considered strongly affected by motion and a score S(s,v) in the range of [1.0...2.0] is calculated according to
where f1 was set to 0.7. To account for lower signal intensities due to attenuation caused by the diffusion weighting, the background threshold t(v) for each volume v was scaled according to
where t(1) was set to 100, b(v) is the b-value in s/mm2 for each volume and D was set to 1.0*10-3 mm2/s. The baseline background threshold t(v) has to be adjusted once to the overall signal intensity level of the reconstructed images.
Slices for which the number of pixels above threshold r(s,v) was less than 5% compared to the total number of pixels per slice were ignored for score calculation and their score was set to 0.0. Slices that met this criterion were either top- or bottom-most slices with very little visible anatomy.
If the slice was not considered affected by motion based on the magnitude data, the phase data was examined. The slope of the phase data was determined in the row and column with the maximum number of contiguous pixels above background threshold by performing a linear fit to the unwrapped phase values. The threshold guarantees that the fitted data points are within the brain tissue. A minimum number of 9 contiguous pixels were required for slope calculation. The selection of contiguous pixels ensured that a proper fit could be done. Finally, if the maximum of the absolute values of the phase slopes for row and column p'(s,v) was larger than the threshold for the phase slope f2 = 0.3, the score S(s,v) was computed in the range of [0.0...1.0] as
The threshold f2 was selected to exclude small phase slope values that are caused by eddy currents even without motion.
Two experiments were performed on two volunteers each:
Four scenarios were tested for each experiment: a) no motion & no correction, b) no motion & correction, c) motion & no correction, d) motion & correction, where “correction” implies DW-PACE and reacquisition (Table 1). This resulted in two data sets per scenario with correction: i) original acquired data and calculated diffusion maps; ii) data where corrupted images are replaced with reacquired data and calculated diffusion maps. Diffusion maps were calculated online on the scanner during image reconstruction. Note that motion estimation was always performed even in the “no correction” scenario and estimated motion values were stored in a log file. However, in this case the motion estimates were not fed back to the scanner and not applied during the image reconstruction. This allowed comparison of DW-PACE motion estimates to those estimated by post-processing software. To demonstrate the effects of rotations and translations, both experiments with all four scenarios were performed twice: once with the subjects instructed to rotate their head about the z-axis (paradigm ‘R’) and once with the subjects instructed to shift their head along the z-axis (paradigm ‘T’).
The proposed method was also used to image 23 subjects with recurrent glioblastoma who had experienced prior treatment failure. All patients had histological confirmation of glioblastoma by either biopsy or resection. These subjects were instructed to hold still during the acquisition, as it would be done in a clinical setting.
The MRI scans performed for these experiments were approved by the local Institutional Review Board. All subjects provided informed written consent prior to scanning. All experiments were performed on a 3 T Tim Trio (Siemens, Erlangen, Germany) using a 32-channel head coil. Acquisition parameters were as follows:
Auto-align (23,24) scans were acquired before each scan to ensure the same slice prescription relative to the head position at the start of each scan. Acquisition parameters for the scans on the tumor subjects were as follows: TR 7.9 s, no TR delay, TE 82 ms, FoV 237 mm, 64 slices, 1.85 mm slice thickness, no gap, matrix size 128×128 (image pixel size 1.85 mm × 1.85 mm), bandwidth 1396 Hz/px, 6/8 phase partial Fourier, acceleration factor 2 using GRAPPA, 7 non diffusion-weighted volumes, 42 diffusion gradient directions, b-value 700 s/mm2, 6.2% reacquisitions (maximum of 194 additional images), scan time 6:54 min:s without reacquisition. No triggering was used for any of these acquisitions.
The data from experiments 1 and 2 was processed as follows. Motion correction values as determined with DW-PACE were compared to retrospective motion correction using AFNI (25). AFNI uses gradient descent on a least squares cost function for fast 3D rigid-body image alignment. Mean and standard deviation of the absolute difference between all translation and rotation values as determined by DW-PACE and AFNI were calculated for the cases in which motion correction was disabled i.e. for scenarios (a – no motion & no correction) and (c – motion & no correction). Mean and standard deviation of the absolute residual translation and rotation values as determined by AFNI were calculated for the cases in which motion correction was enabled i.e. for scenarios (b – no motion & correction) and (d – motion & correction). In both cases, smaller values indicate a closer match between the motion calculated by AFNI and the motion calculated by DW-PACE. To assess the performance of DW-PACE, mean squared errors (MSE) between scenario (a – no motion & no correction) and all other scenarios were calculated. This was done on a slice-by-slice basis as well as for the whole volume based on the trace of the diffusion-weighted images. Deterministic fiber tracking and display was performed on all data from experiment 2 using Diffusion Toolkit and TrackVis (26).
The data from the tumor subject was processed as follows. Maximum absolute translations and rotations per TR were calculated from the motion correction values determined by DW-PACE. Maximum absolute residual translations and rotations were determined using AFNI based on the data without reacquisition.
Motion correction with DW-PACE worked reliably in all experiments i.e. very little motion was found when the subject did not move their head and motion correction values estimated by DW-PACE were comparable to those estimated during post-processing using AFNI (Figure 1a, Figure 1b). Table 2 shows the absolute difference between translation and rotation values estimated by DW-PACE and AFNI for scenarios (a – no motion & no correction) and (c – motion & no correction). All mean values are less than 0.9 mm translation and less than 1 degree rotation. In the case where DW-PACE was applied, remaining motion estimated using AFNI was small suggesting that DW-PACE and AFNI co-registration perform similarly well. Table 3 shows the absolute residual translation and rotation values estimated by AFNI for scenarios (b – no motion & correction) and (d – motion & correction). All mean values are less than 0.6 mm translation and 0.5 degree rotation, indicating that DW-PACE motion correction worked well. Differences of the estimated motion values of the two scenarios with motion per experiment and paradigm are because the translational and rotational motion could not be repeated perfectly by the subjects.
Image quality of the calculated diffusion maps was comparable for scenarios (a – no motion & no correction), (b – no motion & correction) and (d – motion & correction) i.e. when either no motion was present or when motion was corrected using DW-PACE (Figure 2). Even without subject motion DW-PACE therefore does not degrade image quality. In experiments with subject motion and without DW-PACE, image quality was degraded, visible as image blurring and edge artifacts in the calculated diffusion maps. In addition, the color-coded FA map shows color changes indicating erroneous tensor values from the misaligned data.
MSE of the trace of the diffusion-weighted data between scenario (a – no motion & no correction) and the other scenarios increased for scenario (c – motion & no correction) compared to scenario (b – no motion & correction) (Table 4, Figure 3). With DW-PACE, MSE decreased for scenario (d – motion & correction) compared to scenario (c – motion & no correction). MSE values include errors from auto-align inaccuracies, subject motion between the auto-align scan and the diffusion scan, subject motion during the diffusion scan, background noise as well as errors caused by DW-PACE. These errors are reflected in the non-zero MSE values for scenario (b – no motion & correction).
Figure 4 shows examples of slices acquired during subject motion as well as the corresponding reacquired images. If severe motion was present during acquisition, signal dropouts in the magnitude data as well as strong phase gradients are visible. If the motion was small, the magnitude data appears unaffected although subtle signal attenuation is often present. However, the phase data still shows a larger phase gradient compared to data without motion, thereby allowing the determination of the severity of motion. Overall, motion detection and scoring worked well to identify images affected by subject motion. Without reacquisition, artifacts are clearly visible in the diffusion maps (Figure 5). After replacing the images with the most severe artifacts - as indicated by higher scores - by the reacquired images, the resulting diffusion maps were comparable to the ones calculated from the scans without motion.
As for experiment 1 above, MSE of the trace of the diffusion-weighted data between scenario (a – no motion & no correction) and the other scenarios increased for scenario (c – motion & no correction) compared to scenario (b – no motion & correction) and decreased for scenario (d – motion & correction) compared to scenario (c – motion & no correction) (Table 4, Figure 3). The magnitude of the MSE values increased for scenario (c – motion & no correction) due to the images affected by motion artifacts. Reacquisition always resulted in smaller MSE values with two exceptions where marginally higher values were found: scenario (b – no motion & correction), paradigm ‘R’ of subject 1; and scenario (c – motion & no correction), paradigm ‘R’ of subject 2. MSE values also include any non-corrected motion artifacts in cases where more images were affected by motion than were reacquired.
Similar effects can be seen in the fiber tracking results (data not shown). If the images with motion artifacts are included in the fiber tracking processing, severe errors in determination of the fiber orientation are visible in the affected slices as well as false propagation of fibers crossing these slices. Once the corrupted images are replaced by reacquired images, the fiber tracking results are comparable to those calculated from experiments without subject motion.
The range of motion calculated from the tumor subject data was small compared to experiments 1 and 2. Maximum absolute translations per TR determined by DW-PACE were 0.84 mm, 0.54 mm, and 0.74 mm for x, y, and z-axis, respectively. Maximum absolute rotations per TR determined by DW- PACE were 0.56°, 1.43°, and 1.79° for x, y, and z-axis, respectively. Maximum absolute residual translations determined by AFNI were 0.46 mm, 0.21 mm, and 0.47 mm for x, y, and z-axis, respectively. Maximum absolute residual rotations determined by AFNI were 0.35°, 0.71°, and 0.76° for x, y, and z-axis, respectively.
Figure 6 (top box) shows an example of motion related signal attenuation in the diffusion-weighted magnitude image from one of the tumor subjects. The phase image shows rapid 2π-transitions in phase. The reacquired magnitude image is free of signal attenuation and the phase image shows nearly flat phase values. The effect of reacquisition is visible in the color-coded FA maps in which the erroneous blue tint disappears. Data from another tumor subject is shown in the middle box. A concentric phase pattern can be seen in the phase image around the brain stem, possibly caused by CSF pulsation. The reacquired phase image appears free of this pattern. Since this effect is very subtle, there is no visible difference between the diffusion-weighted magnitude images and the color-coded FA maps before and after reacquisition. Interestingly, the proposed method also appears to be suited to capture data corrupted by spikes (bottom box). In this tumor subject, spike-related artifacts that were clearly visible in the diffusion-weighted magnitude and phase images affected two slices. These slices were marked for reacquisition and the reacquired data did not show any spike-related artifacts. Subtle differences can be appreciated in the color-coded FA maps.
Motion is a major cause of artifacts in DWI in clinical and research settings. Physiologic motion (respiration, vessel and CSF pulsation) usually causes only minor artifacts in specific locations of the brain and can be mitigated with existing equipment in a straightforward manner by use of gating at the cost of increased acquisition time (10–12). However, there is currently no standard technology available on MRI scanners to correct for slow and fast bulk subject motion causing mis-registration, image blurring, edge artifacts and signal dropout artifacts.
Retrospective correction of slow bulk motion in diffusion data can be done e.g. by using FLIRT (27,28). If motion correction is done offline, care has to be taken to reorient the diffusion gradient direction i.e. b-matrix before processing (13). If the motion correction is done prospectively and the diffusion gradient directions are corrected with the image acquisition as demonstrated here, the errors without reorientation of the diffusion gradients would be comparatively small, limited to the amount of motion from one TR to the next.
In contrast to slow bulk motion, signal dropout artifacts caused by fast bulk motion cannot be corrected retrospectively. To prevent artifacts in the calculated diffusion maps and in fiber tracking, any corrupted images have to be discarded before processing. This leads to loss in SNR and bias in the tensor calculation if the remaining diffusion gradient directions are unevenly distributed. In the case of acquisitions with only one average of six diffusion gradient directions, the data becomes unusable for tensor analysis. To mitigate these effects, specially crafted sets of diffusion gradient directions have been proposed (29). Manual selection of corrupted images is also cumbersome and error prone since artifacts are not always easily visible in magnitude data and phase data is rarely available.
In multi-shot acquisitions, motion artifacts can be corrected using various methods (5,6,30–34). However, these do not apply to single-shot EPI. If a slightly longer acquisition time is acceptable and the percentage of corrupted images is low enough, these corrupted images can be reacquired during or at the end of the acquisition. This requires identification of affected images in real-time as has been shown before for different sequence types (31), single-shot EPI (18), multi-shot methods (17) and readout-segmented EPI with navigators (6). Our proposed reacquisition method is based on magnitude and phase data compared to the echo peak location and magnitude in k-space as described in Shi et al. (18). Shi's method requires acquisition of multiple averages where the first average serves as reference for the corresponding further averages of non diffusion-weighted and diffusion-weighted volumes. In contrast, our method requires only the first non diffusion-weighted volume as reference volume for estimating the corruption score, defining, which images need to be reacquired.
The values of the score calculation parameters were based on a number of retrospectively analyzed diffusion data sets acquired with different protocols (data not shown). The value of the background threshold factor f1 for the magnitude data was selected to capture only the most severe artifacts due to large bulk motion. In this case, the score calculation based on phase data would fail because of rapid phase wrap. The score calculation based on the phase data captures less severe artifacts. However, it was not specifically designed to address subtle signal loss artifacts like those caused by CSF pulsation. In these cases, non-linear phase errors may confound the linear fitting process and a higher order fit may provide more robust estimates. The adjustment of the background threshold from non diffusion-weighted images to diffusion-weighted images using the background scaling factor D may cause images to be labeled as corrupt in subjects with a very large fraction of CSF compared to overall brain volume e.g. in subjects with severe atrophy or hydrocephalus. The value for the phase slope threshold f2 was selected to exclude small phase slope values caused by eddy currents even without motion. The magnitude of the effect varies depending on the diffusion gradient direction and is likely to be different from scanner to scanner depending on the gradient coil design and eddy current compensation being used. More accurate estimates for f2 could be determined with a phantom calibration step. However, depending on the sensitivity to the diffusion gradient direction, this calibration step would have to be repeated for each protocol with a different set of diffusion gradient directions and was deemed impractical. The 5% threshold that excludes slices with very little or no anatomy was also based on the retrospectively analyzed data. Variation of this threshold within a reasonable range will have little effect on the scoring. The three different protocols used in this study demonstrate that the selected score parameters work well for a range of protocols. Overall, the values of the score calculation parameters mainly determine how many images will be labeled as affected by motion and lead to an overall shift in scores while the relative order of scores would be mostly unaffected. A change in the parameters would therefore affect the likelihood of the maximum allowed number of reacquisitions being played out or the reacquisition phase ending early because fewer images were labeled as corrupt.
The proposed DW-PACE method could be improved in several ways. For data sets with many slices and high resolution and therefore long TR, the motion correction update rate can be slow since the update is done only once per TR. Providing feedback up to once per slice once the first two complete volumes are acquired could increase the update rate. Every newly acquired slice would then replace the corresponding slice in the target volume before co-registration and motion correction is applied. The update rate could also be increased by a factor of at least two with the use of simultaneous image refocusing (SIR) (35) or controlled aliasing in parallel imaging results in higher acceleration (CAIPIRINHA) (36), since these methods allow acquisition of multiple slices per diffusion encoding.
Since separate reference volumes are used for each b-value, the first (i.e. reference) volume of each b-value step is not motion corrected. To allow co-registration of volumes acquired with different b-values, a different co-registration algorithm would have to be used that is insensitive to the changes caused by different b-values e.g. a mutual information-based method. Mutual information-based methods might also provide more accurate results when aligning diffusion-weighted data from different diffusion gradient directions. However, the algorithm would have to be fast enough to provide feedback within a short time frame, typically about a few tens of milliseconds. In scans where the data throughput is very high (e.g. many slices with high resolution or a high number of receive coil elements) and image reconstruction is complex, e.g. due to the use of acceleration using GRAPPA or sensitivity encoding (SENSE) (37), the image reconstruction and following co-registration might not be performed fast enough to provide updated values before the next volume. In these cases, motion correction may not be applied to each volume.
Finally, the DW-PACE method is only applicable as long as the SNR of the images is sufficient and the diffusion contrast does not change too much from diffusion gradient direction to diffusion gradient direction. This limits the application of the DW-PACE method to b-values below approximately 2,000 s/mm2. Most clinical application fall in this range, however, Q-space imaging (QSI) (38), Q-ball imaging (39) and diffusion spectrum imaging (DSI) (40,41) are incompatible with this method. For scans with high b-values and therefore little remaining tissue signal i.e. low SNR, motion tracking would have to rely on either separate navigators or external motion tracking devices e.g. by optical means (14–16). In contrast, the proposed method does not require any external hardware or software and acquisition time is only increased by the amount of time needed to reacquire the corrupted images.
The combination of motion correction using DW-PACE with reacquisition worked well for all subjects included in this study. It is expected to be beneficial in patient and other subject populations, as long as motion is slow or limited to short fast motions sufficiently separated in time so that the selected maximum number of reacquired images is enough to replace the corrupted images. Given the limited range of experiments, it should be noted that this is a proof-of-principle study and that further studies are warranted to fully understand the best application of this method. In conclusion, the proposed method allows calculation of artifact free diffusion maps in the presence of slow and fast bulk subject motion at the cost of only slightly increased scan time.
The authors thank Drs. Tracy T. Batchelor and Elizabeth Gerstner for help with patient recruitment. This research was carried out in whole or in part at the Athinoula A. Martinos Center for Biomedical Imaging at the Massachusetts General Hospital, using resources provided by the Center for Functional Neuroimaging Technologies, P41RR14075, a P41 Regional Resource supported by the Biomedical Technology Program of the National Center for Research Resources (NCRR), National Institutes of Health. This work also involved the use of instrumentation supported by the NCRR Shared Instrumentation Grant Program and/or High-End Instrumentation Grant Program; specifically, grant numbers S10RR021110, S10RR023401, S10RR019307, S10RR019254, S10RR023043. This study was supported by the National Cancer Institute, grant number 5R01CA137254-02.