|Home | About | Journals | Submit | Contact Us | Français|
Short-axis PROPELLER-EPI (SAP-EPI) has been proven to be very effective in providing high-resolution diffusion-weighted and diffusion tensor data. The self-navigation capabilities of SAP-EPI allow one to correct for motion, phase errors, and geometric distortion. However, in the presence of patient motion, the change in the effective diffusion-encoding direction (i.e. the b-matrix) between successive PROPELLER ‘blades’ can decrease the accuracy of the estimated diffusion tensors, which might result in erroneous reconstruction of white matter tracts in the brain. In this study, we investigate the effects of alterations in the b-matrix as a result of patient motion on the example of SAP-EPI DTI and eliminate these effects by incorporating our novel single-step non-linear diffusion tensor estimation scheme into the SAP-EPI post-processing procedure. Our simulations and in-vivo studies showed that, in the presence of patient motion, correcting the b-matrix is necessary in order to get more accurate diffusion tensor and white matter pathway reconstructions.
Diffusion-weighted imaging (DWI) is a useful MR-imaging modality that can be used to detect ischemic stroke (1). Moreover, diffusion/perfusion mismatch has been used as a method to identify patients that would benefit from reperfusion therapy (2). An extension to DWI is Diffusion Tensor Imaging (DTI) that describes the anisotropic diffusion behavior of water molecules in biologic tissue, whereby multiple diffusion-weighted images are acquired to derive the full diffusion tensor pixel-by-pixel (3). Thus far, high-resolution DTI has been challenging because of ghosting artifacts associated with multi-shot methods that are warranted to achieve such high spatial resolution. That is, miniscule patient motion (e.g. brain pulsation adjacent to the lateral and fourth ventricles) causes each k-space segment that follows a diffusion-preparation period to have a random and spatially varying phase. These spatially-variant phase-errors cause the ghosting artifacts characteristic for multi-shot DWI. To eliminate these ghosting artifacts, the unwanted non-linear phase should be removed from each shot prior to the reconstruction of a diffusion-weighted image (4). Several navigator-based methods have been proposed to correct for the effects of these phase errors in high-resolution DWI and DTI (5–11).
In addition to the phase errors, another disadvantage of the use of multi-shot imaging for high-resolution DTI is that the overall scan time is longer and thus the chances for patient motion increase. This can become problematic, especially in children, the elderly, and patients with certain medical conditions that prevent them from remaining still during the course of the examination (e.g. stroke or Parkinson’s disease). The main challenge with gross patient motion in DTI is that the following unfavorable effects must be taken into consideration: (i) pixel misregistration between different segments, and (ii) the subject is exposed to an unwanted and slightly different diffusion-encoding for each k-space segment (or blade in the case of PROPELLER). The latter causes the b-matrix to vary over k-space segments. Whilst the correction for these altered b-matrices is straightforward to implement for single-shot data during the tensor processing phase (12,13), matters are more complicated for multi-shot or interleaved acquisitions. Here, reconstruction of the diffusion tensors requires non-linear methods to estimate diffusion tensors directly from the complex k-space data. In this respect, a rigorous description of the latter approach has been shown recently elsewhere (14).
Since PROPELLER variants enjoy increasing popularity – due mostly to their self-navigating property – the purpose of this study was to apply b-matrix correction to high-resolution DTI, in particular to short-axis PROPELLER-EPI (SAP-EPI) DTI (15). Recently, the SAP-EPI approach was suggested as a very effective technique to address issues that otherwise would considerably impair high-resolution DWI and DTI. Namely, compared to single-shot EPI, SAP-EPI generally has low geometric distortion levels and better spatial resolution, and the repetitive sampling of the central k-space region allows one to perform phase and gross patient motion correction without the need for extra navigators (16). Compared to RF-refocused PROPELLER, SAP-EPI has considerably less SAR limitations and is much quicker at the price of slightly increased distortion levels (17). The focus of this study was to augment multi-shot diffusion tensor SAP-EPI by the aforementioned b-matrix correction approach and emphasize its requirement for high-quality DTI and tract tracing.
The post-processing chain for SAP-EPI is given in Figure 1 using R EPI-interleaves per blade to further diminish distortions. First, for each blade the Nyquist ghost parameters (19,20) and GRAPPA (generalized autocalibrating partially parallel acquisitions) weights (21) were estimated from the non-diffusion-weighted images (i.e. b = 0 images). The ghost parameters were estimated by finding the linear phase term between even and odd echoes that minimize the image entropy. Using the image entropy as a cost functions is equivalent to confining the image to as few pixels as possible, which in turn aims to reduce FOV/2 ghosts that arise due to the delays between even and odd echoes. The phase terms that minimize the image entropy were found using the iterative Gauss-Newton technique. For each iteration, the phase map was applied to the image in the x-ky domain, 1D Inverse Fourier Transform was taken in the y direction, and the entropy of the resulting image I was calculated. Following the estimation of ghost and GRAPPA parameters from the b = 0 image, these parameters were applied to each of the R EPI-interleaves of the corresponding blade of the b = 0 and diffusion-weighted images. In particular, R in our study was 3 and equaled the GRAPPA acceleration factor. After GRAPPA reconstruction and POCS (projection onto convex sets) processing, a total of R full blades were generated (one fully sampled blade for each interleave). Each of these R full blades was FFT’ed to image space and then phase-corrected using a triangular window (22). Thereafter, for each of the R blades the resulting coil element images were combined using complex sum-of-squares algorithm (23) yielding R independent complex ‘blade’ images. (Do note that this last step reduces the amount of data to be processed by a factor equal to the number of coils.) Next, the phase and ghost-corrected R blade images were inverse FFT’ed and the resulting raw k-space data were fed into the non-linear conjugate gradient algorithm together with all the data acquired at the other blade angles. This non-linear approach will be explained briefly in the next section.
Rotational motion of patient’s head causes each blade to be encoded with a different b-matrix (Fig. 2). This makes it incorrect to combine different blades to form a single diffusion-weighted image even when conventional motion correction to realign morphologic features, such as suggested in (22), is applied on the data. A possible solution is to estimate the diffusion tensors directly from the complex k-space data (14). In the presence of gross motion, the acquired k-space data is given by the following equation:
Here, k is the k-space location, r is the image space location, d(k) is the acquired data, m(r) is the non-diffusion weighted image, D is the 3 × 3 diffusion tensor, b is the 3 × 3 b-matrix, and R and Δr are the rotation and translation matrices, respectively. The diffusion tensors D and the non-diffusion weighted image m are estimated by minimizing the sum-of-squares difference between the acquired and estimated k-space data:
To solve for D and m, the non-linear conjugate gradient algorithm was used (http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf). Details of this step are described in (14).
Experiments were conducted on healthy volunteers with 1.5 T and 3 T whole-body MRI units (GE Healthcare) and an eight-channel head coil. All experiments were approved by the institutional review board and written informed consent was obtained from each participant prior to the study.
To simulate noticeable patient movement: (i) two DTI acquisitions on a volunteer were carried out; (ii) between these two scans, the volunteer was asked to rotate his head by ~10°; (iii) then, the data from both datasets were mixed so that every second blade came from the second dataset. The following acquisition parameters were used for the SAP-EPI DTI sequence at 1.5T: 9 blades with a blade width of 64, a target in-plane resolution and reconstruction matrix of 288 × 288, a GRAPPA-acceleration factor = number of EPI interleaves = R = 3, partial Fourier encoding with 18 overscans (i.e. partial Fourier factor = 56%), a slice thickness of 5 mm, TR = 3s, TE = 67.8 ms, FOV = 24 cm, 30 slices, total scan time = 24 min. A twice-refocused diffusion preparation was used with a b-value = 1,000 s/mm2. 15 isotropically distributed diffusion-encoding directions and 3 b = 0 images were acquired.
The first out of the two SAP-EPI DTI acquisitions was used as a quality reference for the assessment of the performance of different reconstruction methods. Specifically, the reconstruction of diffusion tensors from the simulated dataset was carried out using 3 methods:
A region-of-interest (ROI) around the splenium of corpus callosum was selected for analysis. The different reconstruction methods were compared based on three metrics: 1) Average FA value inside the ROI; 2) Average mean diffusivity inside the ROI; and 3) Average angular deviation (in degrees) of the primary eigenvectors from the reference orientations (after correction of bulk rotation).
Three DTI scans were carried out on another volunteer. The same SAP-EPI sequence parameters that had been used for simulation studies were used for this experiment at 3T. For the first experiment, the volunteer was asked to perform a single head rotation in the middle of the scan (and remain at this new location), whilst for the second experiment the volunteer performed multiple head rotations throughout the entire scan. An additional third dataset was acquired where the patient was asked to stay stationary. This data set served as our quality reference. Image reconstructions and analysis of reconstructed data were carried out using the three metrics described for the simulation study.
Two separate scans were carried out on a third volunteer. For the first scan, the volunteer was asked to stay still, whilst for the second scan the volunteer was asked to rotate the head and remain at this location once every 3 min. The first dataset was reconstructed with no correction and served as the reference data for further analysis. The second dataset was reconstructed using the three reconstruction methods explained above. Specifically, the following parameters were used for the SAP-EPI sequence: five blades with a blade width of 64, a target in-plane resolution and reconstruction matrix of 192 × 192, a GRAPPA-acceleration factor R = 3, NEX = 3, partial Fourier encoding with 18 overscans (i.e. partial Fourier factor = 56%), a slice thickness of 3 mm, 38 slices, TR = 4 s, TE = 75.9 ms, FOV = 24 cm. Twenty-four isotropically distributed diffusion-encoding directions with b = 1,000 s/mm2 and 3 b = 0 images were acquired. The total scan time was 27 min.
Fiber tractography was performed using the visualization and tracking software called SmartTrack that was built in-house (24). Euler’s method was used for fiber tracking (step size = 0.6 mm, curvature threshold = 40°, FA threshold = 0.15). A seed region consisting of ~100 seeds was chosen in the splenium of corpus callosum. Fiber tracking (using always the same seed region) was done on the motion-corrupted data set using the three reconstructions (i.e. Methods A, B, and C) and on the reference data set (without motion). Do note for the reference data set Method A has been used for all analysis. No significant differences between FA and eigenvector values have been found between Method A and C in the absence of motion. For each seed, fibers obtained from the three different reconstructions of the motion-corrupted dataset were pairwise compared to corresponding reconstructions using the reference data set without motion by employing two metrics:
Results of simulations are shown in Figure 3 and Table 1. The FA maps reconstructed without correction (Method A) suffered from detrimental motion artifacts (Fig. 3b). These artifacts were mostly removed after motion-correction (Method B) or after our new method (Method C) had been applied (Fig. 3c). Likewise, the FA maps derived from reconstructions performed with motion correction only (Method B) were qualitatively similar to the ones reconstructed without motion and b-matrix correction (Method C) (Fig. 3c, 3d). However, it was also observed that the mean FA value inside the ROI was slightly increased and closer to the reference value when b-matrix correction was applied (Table 1). Furthermore, in terms of eigenvector orientations, it was observed that the orientations of the eigenvectors reconstructed by the application of b-matrix correction (Method C) were closer to the true orientations compared to the case where no b-matrix correction (Methods A and B) was applied (Fig. 3f, 3g). Table 1 also shows that the mean deviation of the eigenvectors from their true orientations was reduced.
Similar results were observed for the in-vivo studies. Without any correction (Method A), significant motion artifacts were observed on FA maps, which were removed after motion-correction had been applied (Methods B and C) (Figs. 4b, 4c, 4h, 4i). However, for both discrete and multiple head rotations, the application of motion and b-matrix-correction (Method C) reduced the deviation of the primary eigenvector orientation between the reference and reconstructed major eigenvectors (Figs. 4f, 4g, 4l, 4m), whilst using motion-correction only (Method B) left a considerable orientational mismatch. Again, for the case where the subject performed a single motion, after the application of b-matrix correction (Method C) an increase in the average FA value was observed (Table 2) that was also closer to the reference data. In the case of multiple head rotations, no change in average FA was observed (Table 2).
The results of fiber tracking experiments are shown in Figure 5. The top row shows the visualization of the fiber tracts. The tracts obtained from the non-motion corrupted dataset (i.e. the reference data) are shown in gold. The tracts obtained from the datasets reconstructed either with no correction (method A), motion correction (method B), or motion and b-matrix correction (method C) are shown in blue. It can be seen that the fibers reconstructed with motion and b-matrix correction (Fig. 5c) agree best with the fibers from the reference data. The second row shows the deviation of the fibers reconstructed with Method A, B, and C from the fibers obtained from the reference dataset. The fibers obtained from the dataset with b-matrix correction (Method C) show less deviation compared to other methods (Fig. 5f). The mean Hausdorff distances over all fibers between the fibers reconstructed from the reference data and the motion corrupted dataset were found to be 11.2 ± 4.5, 15.0 ± 4.5 and 5.3 ± 1.4 mm for methods A, B, and C, respectively.
Due to the prolonged acquisition time and the increased chance of involuntary patient motion, correction of motion artifacts in high-resolution DTI is paramount for good image quality. So far, efforts to correct for motion in DTI have mainly been limited to compensate for random phase terms due to microscopic motion, as the random phase errors had the most dramatic effects on image quality for multi-shot imaging. The effect of macroscopic subject motion – especially its effect on the effective diffusion-encoding direction (i.e. the b-matrix) for a particular segment – has usually been neglected. This is critical, as it has been shown that a drifting motion of the head is often observed for long DTI scans, even with cooperative patients. Even for small drifts, not correcting for the b-matrix can cause erroneous estimation of fiber tract pathways (13). For high-resolution DTI, the correction of the b-matrix due to motion is challenging due to the fact that every segment (i.e. blade) in k-space is encoded with a different b-matrix. This violates the encoding-consistency requirement that fundamentally underlies the MR acquisition process and makes it impossible to estimate a diffusion-weighted image that corresponds to a single diffusion-encoding direction. The net image is more a mixture of different diffusion-encoding directions.
The effects of performing b-matrix correction on low-resolution DTI scans using variable-density spiral readouts were demonstrated previously (14). That particular study demonstrated the application of non-linear estimation of diffusion tensors on spiral scans in order to correct for the rotation of b-matrix.
PROPELLER-type acquisitions are useful in high-resolution DTI due to their self- navigating capabilities and fast acquisition speeds. SAP-EPI is a variant of PROPELLER and it has been introduced to perform high-resolution DTI/DWI with little geometric distortion and excellent self-navigation capabilities, the latter due to the repetitive sampling of the central portion of the k-space with each blade. SAP-EPI’s low geometric distortions are due to the considerably higher readout bandwidth (compared to single-shot EPI) in the phase-encode direction together with the application of parallel imaging. An important advantage of SAP-EPI is that geometric distortions and final image resolution are decoupled from each other. For SAP-EPI, the blade width determines the distortion level and the blade length the final resolution. For ss-EPI, which is the most commonly used pulse sequence for DTI, higher resolution also results in higher geometric distortions.
As with most high-resolution techniques with reduced distortions, there will be a penalty in scan time. With SAP-EPI, the scan time will be increased by a factor that is approximately equal to the number of blades (in our case it is 5–9 blades). The SNR is also dependent on other factors such as blade width and the number of slices that can fit in a given TR. However, due to the pseudo-EPI nature of the trajectory, the scan time increase for SAP-EPI is smaller than for other high-resolution FSE-based techniques such as Radial FSE or PROPELLER FSE. Moreover, distortion correction techniques can be used for SAP-EPI, which do not require extra scan time. This allows the blade width to be increased, with slightly increased distortions in favor of a higher SNR efficiency.
As shown in this paper, the use of SAP-EPI with b-matrix correction represents the next logical step to accomplish very high-resolution, nearly distortion-free, motion and b-matrix corrected DTI. The most striking results can be seen from our simulation studies that showed that the application of b-matrix correction improved the accuracy of the orientation of the major eigenvectors (Method C vs Methods A and B). We also observed a slight increase in FA values when Method C was used. This finding is consistent with previous studies (26) and is evident when one considers that neglecting the b-matrix will average the tensor sensitivity and thus lead to a more isotropic tensor estimate. In the simulation studies, the fully reconstructed blades acquired with two distinct head orientations were intermixed after performing ghost-correction and GRAPPA reconstruction. It is important to notice that since there was no intra-acquisition motion, the effects of motion between the calibration scan (i.e. the first T2w image) and the DWI scans were eliminated. This was an advantage for the simulation studies since the effects of b-matrix-correction were completely isolated from any errors which may arise from motion occurring between the calibration scan and diffusion weighted scans. For all studies it was also assumed that there was microscopic motion but no macroscopic motion during the diffusion-encoding period. Lastly, our approach was restricted to in-plane motion, although the b-matrix correction can be expanded to all three spatial dimensions. Thus, the next logical step should be the expansion of the gross motion correction from in-plane correction to full 3D correction.
Results similar to those from the simulation study were also observed from the in-vivo experiment, but resembled a more realistic scenario in the sense that the subject performed either single or multiple head rotations throughout the scan. Again, a more accurate major tract orientation was obtained with b-matrix correction (Method C vs Methods A and B). One issue with the in-vivo experiments was that there was motion between the b = 0 scan – from which the FOV/2 ghost and GRAPPA calibration information was obtained – and the DWI readouts to which the calibration was applied. Do note that the geometric mismatch would be a serious challenge for SENSE-type reconstructions (27). However, for SAP-EPI this was not a huge problem. First, the FOV/2 ghost parameters do not to change with motion (as they are determined by the gradient set and the orientation of the blade). Second, GRAPPA is robust to motion – at least for the size of coil elements we used in this study (27). One way to fully correct parallel imaging calibration parameters would be to apply SENSE reconstruction instead of GRAPPA in a second pass following correction of coil sensitivities for motion (28).
Despite the fact that motion between the b = 0 acquisition and the diffusion weighted acquisitions is not a limitation, motion within the b = 0 image would result in erroneous estimation of GRAPPA and ghost parameters. A solution to this problem is to find the ‘best’ b = 0 image among all those acquired in the full DTI dataset, and calculate the GRAPPA and ghost parameters from that image (29).
In our SAP-EPI post-processing, the ghost correction corrects for hardware delays between odd and even echoes, but does not correct for any delay in the gradient amplifier that may arise from eddy current effects in the DW images. In this work, we used a twice-refocused diffusion preparation in order to eliminate the eddy currents – thus we did not observe image degradation. However, if one was to do single-refocusing (i.e. not correct for eddy currents in the sequence), the ghost correction technique we use here would still be applicable, and one could use a post-processing technique based on (30), to correct for the affine distortions.
For both simulation and in-vivo experiments, some residual error remained between the actual and reconstructed fiber orientations. This residual error can be attributed mainly to noise and through-plane motion between blades as well as a slight error in the estimation of the motion parameters.
Our fiber tracking experiments (Fig. 5) showed that the correction of the b-matrix is essential to get more accurate white matter pathways. We observed that the fibers obtained from the dataset reconstructed from the dataset with b-matrix correction (Method C) agreed more closely with the fibers obtained from the non-motion corrupted dataset compared to the case whereby no b-matrix correction was applied (Methods A and B). It should be noted that fiber tracking is very sensitive to very small changes in the DTI data (31,32), so some residual deviation remained even after b-matrix correction. This is most likely due to through plane motion and the mismatch between the seed locations.
Previous results on spiral trajectories demonstrated that the reconstruction time (per slice) for spiral DTI data acquired with an 8-channel head array with parameters 128 × 128 resolution, 8 interleaves and 50 iterations was around 4 hours (288 s per iteration) (14). In this case, going up to a resolution of 192 × 192 would result in (192/128)4 times increase in reconstruction time (the number of interleaves should also be increased to match the off-resonance distortion), which would result in a reconstruction time of 20 h (1458 s per iteration). In our case with SAP-EPI, the processing time (on a server with 2 × Intel Xeon 3.2GHz (dual core), 1MB cache, 8GB RAM running linux) to estimate the tensor information was currently on the order of 7 h for a single slice with resolution 192 × 192 and 200 iterations (126 s/iteration). One advantage of using SAP-EPI post-processing compared to the previous study on spiral trajectories (14) was that it was possible to process and subsequently combine the coil data prior to the non-linear estimation of the diffusion tensor elements, which reduced the processing time by a factor equal to the number of coils. It is also possible to reduce the processing time even further by leveraging on the Cartesian properties of SAP-EPI, which will reduce the amount of time spent on performing the gridding and the FFT. This is also an advantage of SAP-EPI over the spiral trajectory. Specifically, for the 288 × 288 resolution dataset given in the text, the use of spiral trajectory would require an FFT size of 288 × 288, which takes 25 ms on our server. Using the SAP-EPI readout with a blade width of 64 requires an FFT size of 288 × 64, which is approximately 5 times faster. The processing time can be reduced even further by using multiple CPUs, computer clusters (as currently used by our vendor for reconstructing clinical images), or graphical processing units (GPUs) (33). While the current reconstruction times of our Matlab-based reconstruction would certainly be prohibitively long for clinical applications, the aforementioned speed-ups should be quite helpful. Moreover, the majority of DTI studies are currently preclinical, thus fast processing can be traded for better data quality by performing reconstructions as overnight batch-jobs.
Short-Axis propeller-EPI was combined with b-matrix correction for diffusion tensor imaging. SAP-EPI post-processing helped us to simplify parts of the reconstruction chain and reduced the amount of data to be processed, which would otherwise become rather time consuming. In simulations and in-vivo studies it was shown that the change in the direction of effective diffusion-encoding direction caused by patient head motion results in the erroneous estimation of major fiber tract orientations and fiber pathways in DTI. This was circumvented by using a novel non-linear tensor estimation scheme. The combination of SAP-EPI readouts in concert with non-linear tensor estimation can provide very high-resolution DTI with minimal geometric distortions and greater accuracy in eigenvector estimation and white matter pathways which otherwise cannot be derived from regularly processed DWIs.
This work was supported in part by the NIH (2R01EB002711, 1R01EB008706, 1R21EB006860), the Center of Advanced MR Technology at Stanford (P41RR09784), Lucas Foundation and Oak Foundation.
We would like to thank Heiko Schmiedeskamp for helpful discussions.
This work was presented in part at the 17th Scientific Meeting of the ISMRM, Honolulu, HI, 2009.