| Home | About | Journals | Submit | Contact Us | Français |
Despite rigid-body realignment to compensate for head motion during an echo-planar imaging (EPI) time-series scan, non-rigid image deformations remain due to changes in the effective shim within the brain as the head moves through the B0 field. The current work presents a combined prospective/retrospective solution to reduce both rigid and non-rigid components of this motion-related image misalignment. Prospective rigid-body correction, where the scan-plane orientation is dynamically updated to track with the subject’s head, is performed using an active marker setup. Retrospective distortion correction is then applied to unwarp the remaining non-rigid image deformations caused by motion-induced field changes. Distortion correction relative to a reference time-frame does not require any additional field mapping scans or models, but rather uses the phase information from the EPI time-series itself. This combined method is applied to compensate EPI scans of volunteers performing in-plane and through-plane head motions, resulting in increased image stability beyond what either prospective or retrospective rigid-body correction alone can achieve. The combined method is also assessed in a BOLD fMRI task, resulting in improved Z-score statistics.
Head motion is a familiar and significant problem for magnetic resonance imaging (MRI) studies such as diffusion, perfusion, and functional MRI (fMRI) that acquire multiple frames of the same brain volume over time. Even a few millimeters of movement is sufficient to cause misalignment of brain volumes over the time-series acquisition (1,2), leading to large, artifactual, signal fluctuations that will obscure the typically small brain activation signal. To maximize the functional sensitivity of these methods, it is therefore necessary to optimize the image stability by realigning the multiple brain volumes before any subsequent statistical analysis.
To achieve this, prospective real-time motion correction is becoming increasingly popular. Several prospective correction strategies have been proposed – including image-based methods (3,4), 3D-navigators (5–7), optical cameras (8–13), and active markers (14–17) – that compensate for motion in the acquisition stage by keeping the scan-plane at a fixed orientation relative to the head. While differing in the proxy used to track the head motion, the focus of the above methods is to correct for bulk motion using a six degrees-of-freedom rigid-body model. Recent works have shown the benefits of prospective methods over more conventional retrospective image realignment techniques, demonstrating improved BOLD fMRI activation (9), diffusion sensitivity (13), and shot-to-shot image stability (17).
While effective in compensating for bulk rigid-body motion, a secondary effect of head motion is non-rigid image distortion caused by changing magnetic field inhomogeneities (18). Local magnetic field inhomogeneities in the brain arise due to internal field inhomogeneities such as around tissue-air or tissue-bone interfaces, as well as imperfections in the external hardware shims. As the head moves and the position and orientation of these local field inhomogeneities change relative to the static B0-field, the effective field experienced within the brain also changes. The single-shot echo-planar imaging (EPI) sequence, preferred for time-series applications, is particularly sensitive to field inhomogeneities due to its low phase-encode bandwidth; as a result, EPI images generally suffer from geometric distortions and intensity nonuniformities along the phase-encode direction (19). Thus, motion-induced field changes will generate time-varying non-rigid geometric distortions in EPI that further reduce image stability, and cannot be corrected for with a rigid-body model.
Several distortion correction methods have been proposed. These methods characterize the field inhomogeneities within the image by means of a field mapping scan, which is then used to correct for EPI distortions in postprocessing. The simplest way to measure field maps is to use the phase differences between images acquired at different echo times (19–21). Alternative reference maps calculated using Fourier methods (22–24), point spread function mapping (25–27), and reversed gradient polarities (28) have also been developed. These field map measurements require additional scan time and/or pulse-sequence modifications, and assume that the subject does not move between the field mapping scan and the EPI time-series. Since distortions vary with head motion, a single field map acquired at the beginning of a scan may not accurately represent the distortions in subsequent EPI time-frames.
In the presence of subject motion, it would be ideal to obtain a field map estimate at each EPI time-frame so that position-specific distortion correction may be performed. A dual echo-time EPI scan has been used to acquire field maps for each time-frame (21), but the added scan time would be impractical for most applications. To avoid this, field maps have been predicted using the rigid-body realignment parameters together with the distortions in the magnitude images themselves (29), as well as from a susceptibility model of the imaged object that requires a priori information (30). A real-time shim update strategy has also been demonstrated (31).
The distortion correction method explored here unwarps all EPI time-frames back to a reference time-frame, as suggested by Jezzard et al. (18). Our motivation for investigating motion-induced distortions comes from ongoing work on prospective motion correction for EPI using active markers (17). Previously, we had observed that despite accurate rigid-body prospective correction during head motions, a residual non-rigid warping of the brain remained along the phase-encode direction. Using only the image phase-information, our preliminary results have shown that it is possible to retrospectively reduce this time-varying warping artifact (17,32).
We propose a combined method that uses prospective motion-correction together with retrospective distortion-correction to obtain images that are compensated for both rigid-body motions and susceptibility-induced non-rigid deformations. First, single-shot EPI time-series data is acquired with slice-by-slice rigid-body prospective correction using active markers (17). Next, non-rigid distortion correction is retrospectively performed using field maps derived from phase-image subtractions. The unwarping algorithm uses the phase information of the actual EPI time-series data, and does not require any additional field mapping scans, susceptibility models, or pulse-sequence modifications. The current work expands on our preliminary distortion correction algorithm (17,32), and uses new EPI time-series data to investigate the relationship between head motion and field changes in the brain. The effectiveness of the combined method is quantitatively evaluated in terms of improvement in image stability, and BOLD fMRI sensitivity in a breath-holding task.
Details of the tracking hardware and scan-plane update scheme (16), as well as their integration into EPI scans (17), have been previously presented and are briefly reviewed here. Motion tracking is performed using three active markers integrated into a headband worn by the subject. Three non-colinear markers are sufficient to fully describe any arbitrary rigid-body head motion. Each marker is a solenoid inductor, tuned and matched to 64.3 MHz, containing a small (3 mm diameter) glass sphere filled with Gd-doped water solution. The markers are attached to a Synergy Multi-Connect (SMC) box (Philips Healthcare, Best, the Netherlands; by IGC Medical Advances), which then connects to the scanner (see ref. (16), Fig. 1).
Real-time, slice-by-slice prospective correction in a single-shot EPI scan is achieved by interleaving a rapid track-and-update module into the imaging sequence before the acquisition of each EPI-slice (see ref. (17), Fig. 1). For every repetition of this module, a short tracking pulse-sequence of orthogonal 1D projection-readouts is used which re-measures the marker positions. During head motion, the rigid-body transformation that realigns the markers to their initial positions (measured at the beginning of the scan) is fed back to adaptively update the scan-plane position and orientation – thus maintaining it fixed relative to the head – before the next EPI-slice is acquired. High accuracy/precision tracking measurements (~ 0.3/0.01 mm (16)) and corrections are accomplished at a temporal resolution (tracking pulse-sequence + calculations + updates ~ 25 ms per module (17)) suitable for real-time application.
A flowchart describing the distortion correction algorithm is shown in Fig. 1. We begin with a prospectively corrected EPI scan of n time-frames, where the nth time-frame may contain motion-induced geometric warping (superscript “w”) relative to the reference (e.g. first) time-frame. The magnitude (Inw) and phase (Φnw) images are reconstructed using the complex data, and phase images unwrapped to remove [−π, π] jump discontinuities (33). Phase-difference maps ΔΦn (radians), which represent the phase accrual between the nth and reference time-frame, are calculated:
Expressing the phase as a frequency-difference map ΔFn (Hz) provides a measure of the relative field inhomogeneity deviation between the nth and reference time-frame due to the different head positions. Since field inhomogeneities cause negligible pixel mis-location in the readout direction for EPI (19), distortion correction then simplifies into a series of 1D pixel-shifts in the phase-encode direction which – due to prospective scan-plane tracking – remains “locked” to the patient orientation. To calculate the necessary pixel-shifts, the local ΔFn is converted into a pixel-shift map ΔWn (pixels) by dividing by the effective spectral bandwidth BWpe (Hz) per pixel Npe in the phase-encode direction.
Unwarping is based on the polynomial-fitting method established by Jezzard et al. (19). The measured field map ΔWn is first fitted with a third order 2D-polynomial Pn. Only pixels in the white matter (WM) and gray matter (GM) of the brain (segmented using Statistical Parametric Mapping v5 (SPM5), Wellcome Trust, London, UK) are used in the least-squares fitting as these pixels represent the tissue of interest, and values where the field is most reliably known. The polynomial surface Pn is masked to yield ΔŴn which defines the area that the distortion transform will be applied. The mask is generated from the GM + WM segmentation, filled, and dilated (disc-element of radius = 2 pixels) to provide distortion correction at the edges of the brain and other pixels where the original data, ΔWn, is not well defined.
The fitted field map ΔŴn is then used to unwarp the nth time-frame Inw back to the reference, yielding Inuw. As both Inw and ΔŴn are estimated in warped space, unwarping resolves into a series of one-dimensional re-gridding operations along the phase-encode direction (y-direction). To calculate each phase-encoding line in Inuw, the signal intensity in Inw at pixel location yw is shifted to yuw =yw + ΔŴn (yw) (19). The Inuw that results is defined on an irregularly sampled grid yuw, which is then sampled onto a regular grid using cubic spline interpolation for sub-pixel shifts. As a final step, image intensity is corrected to first-order using the derivative of the pixel-shift map along the phase-encode direction. For a detailed treatment of the unwarping operation, see references (19,21).
Retrospective distortion correction is performed using two different reference frames. First, Φref = Φ1 (resulting in field map ) is used to unwarp all n time-frames back to the first time-frame. We note, however, that if Φnw is geometrically warped relative to Φ1, a mismatch will be introduced in the pixel-by-pixel subtraction to calculate ΔWn is When is applied in the unwarping operation, a residual error will result in Inuw that is proportional to the degree of geometric distortion in the data. To determine and eliminate this error, we alternately use Φref = 0 (resulting in field map ) to unwarp all n time-frames relative to a common zero-phase reference-frame. No pixel mismatch occurs in this case. The distortion correction algorithm was implemented in Matlab (Mathworks, Natick, MA, USA).
Five EPI studies were performed on four healthy volunteers (one repeat volunteer on separate days). In axial slices, left-right “head-shake” and head-foot “head-nodding” rotations induce in-plane (θip) and through-plane (θtp) motions, respectively, and were designed to investigate the relationship between head orientation and field changes within the brain. For each study the volunteer was trained to reproduce three types of motion throughout the scan: a baseline resting position, a slow but continuous rotation θip ~ ±5°, and θtp ~ ±5°. For each motion, two EPI time-series were acquired – with prospective correction ON and OFF – to evaluate correction quality under similar motion conditions. For scans with correction OFF, all tracking and geometry calculations were performed and logged – but not applied to update the scan-plane – to verify reproducibility of volunteer movement. Volunteers were blinded to whether prospective correction was ON or OFF. All experiments with human subjects were in accordance with local Institutional Review Board (IRB) regulations and informed consent was obtained before each exam.
Experiments were performed on a 1.5T Philips Achieva (Philips Healthcare, Best, The Netherlands). Imaging was performed with a standard quadrature bird-cage coil, and tracking via active-marker headband. Slice-by-slice prospective correction (tracking parameters: TE/TR = 1.8/4.3 ms, α = 4°, resolution = 1 mm) was applied to an axial, single-shot 2D-EPI time-series (imaging parameters: TE/TR = 40/1680 ms, α = 85°, FOV = 192×192 mm, voxels = 3×3 mm, thickness/gap = 5/1 mm, slices = 20, phase-encode direction = right-left, time-frames = 10, BWpe/Npe = 19.5 Hz, shim = first-order (subject-specific), scan time = 38 s). The scanner’s standard post-processing pipeline was used to reconstruct the complex image data. Scan durations increased due to the extensive real-time writing of tracking data to a log-file still required at this developmental stage; with logging disabled, scan times reduce to 23 s. The original scan with prospective correction completely disabled is 18 s.
To compare traditional motion compensation methods with our combined approach, six different processing pipelines were investigated. Three stages of post-processing were performed on scans with prospective correction OFF: N1) no further processing, i.e. completely uncompensated for motion; N2) retrospective rigid-body correction using standard image-based realignment (SPM5); N3) retrospective rigid-body realignment (same as N2)), followed by non-rigid unwarping to the first time-frame using SPM5. Unwarping in SPM5 uses the magnitude images and rigid-body realignment parameters to estimate field maps (29), which are then used to unwarp the time-series. Three post-processing pipelines were also performed on scans acquired with prospective correction ON: Y1) no further processing, i.e. prospective rigid-body correction only; Y2) retrospective non-rigid distortion correction using ; Y3) and using .
To quantify image stability across the time-series, standard deviation maps (σ-maps) were calculated. Binary masks of GM+WM were segmented from the first time-frame of each scan, and applied to the rest of the time-frames so that only voxels containing brain tissue were included. For each voxel, the standard deviation of its intensity across the time-series was then normalized as a percentage of its mean intensity to obtain the σ-map. The average standard deviation ![[sigma with macron]](/corehtml/pmc/pmcents/x03C3x0304.gif)
over all voxels in the σ-map was also calculated. As movement increases the intensity variation of each voxel, the result of perfect motion correction will be reduced values of comparable to resting scans. For visual assessment, difference images were generated by subtracting slices acquired before and after motion, where perfect correction returns a zero difference-image.
An fMRI breath-holding task was performed on three healthy volunteers to further assess the efficacy of our combined method. The breath-hold paradigm was chosen because of its large global BOLD activation, enabling the evaluation of activation statistics in different brain regions. Volunteers performed 10 repetitions of alternating periods of breath holding after inspiration and self-paced breathing. Task blocks were cued by different colored screens in the following order: green screen (14 s) signifying normal breathing; yellow screen (2 s) during which the subjects takes a deep breath in preparation to hold; red screen (16 s) signifying the breath-hold period. After the last red screen, an extra green screen allows the BOLD signal to return to baseline. This task is similar to ref (34). The task was programmed in E-Prime (Psychology Software Tools, PA, USA) and lasted 333 s. As with Experiment I, the subjects repeated slow but continuous “head-nodding” rotations throughout each scan while performing the breath-hold task.
For each subject, the task was repeated 3 times with prospective correction ON and 3 times with prospective correction OFF. In one subject, only 2 scans for both ON and OFF cases were recorded due to hardware malfunction. For each scan, the volunteer’s breathing was recorded using the scanner’s respiratory bellows monitor. The scanning parameters were the same as Experiment I, except slices = 17, time-frames = 111.
The breath-hold task data was put through the standard fMRI-preprocessing stream using SPM. Prospective correction OFF scans are first slice-time corrected, then realigned for retrospective motion correction, normalized into standard space with an isotropic voxel size of 3 mm, and smoothed with a 6 mm FWHM Gaussian kernel. A standard general linear model analysis with a box-car regressor convolved with the canonical hemodynamic response function, as well as the estimated motion parameters as confounds, is used to identify areas of activation caused by the breath hold for each scan. Prospective correction ON scans followed a similar processing stream, except no retrospective motion correction was applied, confounds were generated by the active-marker measurements, and retrospective distortion correction was performed using our combined method.
For each volunteer, a fixed-effects analysis is used to create a single statistical parametric map for each pipeline. Next, a group-level random-effects analysis is used to create the mean Z-score activation map across all subjects for each pipeline. For region-of-interest (ROI) analysis, all ROIs were taken from the Harvard-Oxford Cortical Structural Atlas. For each region, a two-tailed paired t-test, where each voxel pair is a sample point, was used to compare each activation map. Significance for the paired t-test is defined as p < .05/9 Bonferroni corrected, or p < .0056.
Data from Experiment I is summarized in Fig. 2–5a. Figure 2 shows maps at increasing magnitudes of in-plane (Fig. 2a) and through-plane (Fig. 2b) rotation. These maps demonstrate the relative field change, and hence amount of spatial warping, within the brain as a result of motion. For the maximum in-plane rotation shown (|θip| = 5.32°), frequency-shifts did not exceed 5 Hz (or a pixel-shift = 0.75 mm), while for the maximum through-plane rotation (|θtp| = 5.06°) shifts of up to 20 Hz (3 mm) were observed. Figure 2c summarizes the results for all subjects, where each point in the scatter-plot represents a single slice, with the maximum frequency difference across the slice (Δfmax) plotted as a function of rotational magnitude. As rotational magnitude increases, frequency variation across each slice also increases, with the rate of increase significantly greater in the through-plane vs. in-plane direction. This suggests that head-nodding relative to B0 is predominantly responsible for susceptibility gradient variations.
(mean ± SD) between the six pipelines. For in-plane motion (left) and through-plane motion (middle) bar plots, ...Figure 3 shows σ-maps from a single volunteer performing in-plane (row 1) and through-plane (row 2) rotations throughout the time-series. The σ-maps of unprocessed scans (column 2) illustrate the spatial distribution of increased variation caused by the motions. After retrospective realignment (column 3) regions of increased intensity variation remain, particularly at the brain edges and superior slices, and are more pronounced during through-plane motion. Prospective correction alone (column 4) provides significant improvement in image stability. For in-plane motions (column 4, row 1), standard deviations are reduced to spatially-uniform resting levels of ~ 2%, with distortion correction (column 5, row 1) providing no additional benefit. For through-plane motions (column 4, row 2), remaining signal variation in the anterior and posterior regions is still relatively high, but is significantly reduced after distortion correction (column 5, row 2). Percentage improvement in image stability using the combined method vs. prospective correction alone is shown (column 6), with lighter values (relative to the background gray-level) representing pixels where additional distortion correction has reduced σ; regions of improved image stability are visible and consistent with the areas of increased frequency change in Fig. 2a,b, where distortion correction is expected to have the greatest effect.
Difference images from a single time-frame of Fig. 3 are shown in Fig. 4. Difference images were generated by subtracting I1 from images In acquired at a time-frame n near the maximum rotation performed by the volunteer. maps (column 5) at this head position are applied to the corresponding images in column 3, resulting in column 4. For through-plane motions, pronounced non-rigid warping remains in the anterior brain regions after prospective correction (column 3, row 2), emphasized by white/black arrows that denote areas of positive/negative difference, respectively. This pattern indicates an uncorrected non-rigid shear in the right-left phase-encode direction, and is consistent with the regions of increased frequency variation in . The shearing artifact is significantly reduced after distortion correction (column 4, row 2) using .
Image stability from Experiment I, over all subjects and all six pipelines, is summarized in Fig. 5a. For all 10 deliberate motion scans, combined prospective correction + distortion correction with
(Y2) significantly reduced (i.e. improved image stability) vs. traditional image-based retrospective methods (N2, N3) (p < 0.001, two-tailed paired student’s t-test), reducing average over the whole brain by >25% and 36% for in-plane and through-plane scans, respectively; in the frontal pole during through-plane motions where distortions were most severe, ROI analysis showed a reduction of >49%. For all five cases of through-plane motion, Y2 demonstrated superior image stability vs. prospective correction alone (Y1), reducing by ~9% and 26% over the whole brain and in the frontal pole, respectively; for in-plane cases, no significant benefit was found between Y2 vs. Y1. We also note that for all scans Y1 alone significantly reduced vs. either N2 or N3 (p < 0.001) by an average ~28% over the whole brain. No significant decrease in was observed between Y3 vs. Y2 (<5%) as well as N3 vs. N2 (<2%). As expected, N2 and N3 performed poorer in through-plane vs. in-plane cases.
Results from Experiment II are summarized in Fig. 5b. Group-level mean Z-scores are shown for each pipeline. Three ROIs spanning the brain were selected – Frontal Pole (FP), Anterior Cingulate Gyrus (ACG), and Precuneas Cortex (PC), located in the anterior, mid, and posterior regions of the brain, respectively – to compare the group activation-maps. For all three ROIs, Y1 and Y3 both showed statistically significant improvement in Z-scores vs. N2. In the FP and ACG, statistically significant improvement in Z-scores were observed between Y3 vs. Y1, with the FP showing the greatest improvement. In the PC, no statistical significance was found between Y3 vs. Y1. These results are consistent with the image stability analysis in Fig. 5a.
Image stability is critical in EPI-based applications. Motion-induced variation in the time-series data will reduce image stability, and, as a result, reduce sensitivity to signal changes induced by true brain activity. Image stability was improved when prospective rigid-body correction was followed by retrospective non-rigid distortion correction using the phase differences between time-frames. The combined approach was particularly effective in the anterior brain region during through-plane motions, where observed distortions were most severe. These results are consistent with the field changes in (see Fig. 2, ,4),4), and reflect the area of greatest benefit using the combined approach. In this region, improved image stability also translated into the greatest increase in Z-scores in a BOLD fMRI paradigm.
To investigate the pixel-mismatch error in , unwarping was performed to a common zero-phase reference that does not experience this error. Little improvement in image stability was found between these methods (see Fig. 5a, Y3 vs. Y2), suggesting that the error is small. If the phase images Φnw and Φ1 are slowly varying in space, the pixel-mismatch introduced in their subtraction will be small and will closely estimate the actual unwarping map. This was the case for the experimental conditions and ranges of motion tested, as phase-images were smooth and typically δΦnw/δy < 1 rad/pixel, equivalent to a pixel-shift error of <0.2 pixels for every pixel of warping. For greater magnitudes of spatial warping, such as in the presence of larger head rotations, higher image resolution, etc., the zero-phase reference may be more critical in reducing the error. Alternatively, iterative methods to minimize the pixel-mismatch error in may be employed. Note that unwarping using either or only corrects for time-varying field inhomogeneities. To correctly unwarp the static inhomogeneities still present in all time-frames, one more stage of unwarping may be performed on the time-series with a single EPI-volume acquired at a different TE (19).
Other retrospective unwarping schemes have also utilized the EPI phase-information to reduce time-varying distortions (35–38). One set of methods use an “intrinsic” phase-image (36,37) as a reference for unwarping, which requires an additional acquisition at a different TE at the start of the scan. The intrinsic phase can only be measured inside the region of the imaged object at this initial position, relying on extrapolation for voxels outside of the imaged object (and that there are no strong spatial dependencies). This model was shown to be valid for the small range (0.2°) of motion observed (37), but may prove less robust for larger motions relative to the initial position (> 5° in our experiments) where the model would have to rely heavily on the extrapolated phases. An “average” phase-image (38), which does not require an additional TE acquisition, has also been used as a reference for unwarping. Since this reference is generated by averaging all the phase-images in the time-series, large motions would also reduce the constancy of this approach.
A weakness these models (35–37) and ours have in common is that all phase-images in the time-series must be phase-unwrapped to the same DC phase. The need for phase-unwrapping can be reduced (38) if this step is applied after the subtraction step that generates , since Φnw will not normally differ by > 2π from Φ1. These models, and ours, also assume that any phase-contributions to the image, other than B0-field inhomogeneity, are time-invariant. This is shown to be a reasonable assumption in previous studies (35–38). But if, for example, the B1-field phase significantly varies with time and/or subject motion in the coil, then these models will no longer be valid for all time-frames. For our experiments at 1.5T with volume coils, these effects were negligible.
Prospectively corrected data greatly simplifies both the estimate of the unwarping map, and the unwarping operation itself. Since the phase-images have been prospectively realigned, a simple phase-image subtraction provides a reliable estimate of the unwarping map. Data that is not prospectively corrected also faces the challenge that distortions will be in different directions, depending on the current orientation of the head relative to the initially prescribed phase-encode direction. In our approach, since the phase-encode direction prospectively tracks with the head orientation, unwarping simplifies to a series of pixel-shifts in one dimension. The distortion correction algorithm is not limited to an active marker setup, and can effectively augment any prospective correction method given that the complex EPI data is saved.
Grant Sponsor: NIH R21EB006877, NIH/NCRR UL1RR024156, N I H 5R21EB006860, NIH 1R01EB011654
We thank Heiko Schmiedeskamp, Hamed Mojahed, Rafael O’Halloran, Anh Tu Van, and Stefan Skare for helpful discussions.
PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |