PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptNIH Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Magn Reson Med. Author manuscript; available in PMC Nov 1, 2009.
Published in final edited form as:
PMCID: PMC2642633
NIHMSID: NIHMS65179
Sliding Window SENSE Calibration for Reducing Noise in fMRI
Christine S. Law,1,2 Chunlei Liu,1 and Gary H. Glover
1
1Stanford University School of Medicine, Department of Radiology, Center for Advanced MR Technology at Stanford, Stanford, California 94305-5488
2Department of Electrical Engineering, Stanford University, Stanford, California, U.S.A.
Address correspondence to: Christine S. Law, Stanford University School of Medicine, Department of Diagnostic Radiology, Lucas MR Center, Stanford, CA 94305-5488., (650) 723-7577 FAX: (650) 723-5795, cslaw/at/stanford.edu
Functional magnetic resonance imaging (fMRI) at high magnetic field with parallel imaging (PI) has become increasingly popular for high resolution imaging. We present a method of self-calibrated PI-fMRI in which sensitivity profiles are calculated using a sliding window of fully-sampled multi-shot imaging data. We show that by updating these sensitivity profiles in a sliding fashion, thermal noise is reduced in the reconstructed image time series. This is accomplished by retaining thermal noise in the sensitivity profiles; no spatial smoothing is performed. These noisy profiles actually provide a closer match to those required for thermal noise-free reconstruction than conventional sensitivity map generation. Our proposed technique is especially applicable in acquiring high spatial resolution images, where thermal noise exceeds physiological noise. With conventional sensitivity calculation, PI-fMRI sensitivity is preserved only when using a voxel size large enough such that physiological noise predominates. With small voxel size, our technique reveals activation from visual stimulation where conventional sensitivity calculation techniques falter. Our technique enhances fMRI detection especially when higher spatial resolution is desired.
Keywords: Functional magnetic resonance imaging, CG-SENSE, noise, parallel imaging
As the magnetic field is increased, in principle, blood-oxygen-level-dependence (BOLD) imaging can have improved spatial specificity in functional magnetic resonance imaging (fMRI) because more signal originates from tissue and small capillaries rather than from large veins. (1) The fundamental tradeoffs, between spatial and temporal resolution, present many challenges for high resolution and high field fMRI. Higher spatial resolution requires longer readout duration and postponed echo time (TE) that respectively result in increased sensitivity to off-resonance and non-optimum BOLD contrast. Readout duration and TE can be shortened by multi-shot acquisition, but temporal resolution will consequently decrease. Despite reduction of signal-to-noise ratio (SNR), parallel imaging (PI) methods in fMRI are known to be advantageous for reducing signal loss in susceptibility-induced signal dropout regions of the brain. (2) Parallel imaging allows increase in spatial and temporal resolution with shorter TE and is particularly advantageous at higher-field.
Weiger et al. first demonstrated use of a spiral sequence together with parallel imaging. (3) They applied parallel imaging to reduce readout duration in the spiral trajectory, which lessens signal dropout in frontal regions. Consequently, they successfully detected activation in a tasting task stimulating the orbito-frontal cortex. To calibrate coil sensitivity, required by the sensitivity encoding (SENSE) reconstruction algorithm, they gathered an extra fully-sampled image prior to their fMRI experiment. We refer to their technique for calculating sensitivity maps as the reference frame method.
An alternative method for deriving coil sensitivity, from a time-series without a separate image acquisition, was presented by Kellman et al. (4) There, the UNFOLD method, “Unaliasing by Fourier-Encoding the Overlaps Using the Temporal Dimension,” together with SENSE reconstruction was employed to remove spatial and temporal aliasing. This Adaptive Sensitivity Encoding Incorporating Temporal Filtering method (TSENSE) was demonstrated for real-time nonbreath-held cardiac imaging with an acceleration factor R = 2. (Acceleration is a measure of undersampling in k-space.) They acquire odd and even k-space lines in sequential images. Sensitivity profiles are calculated by temporally filtering the aliased images of the same undersampled data. Kellman shows that, with TSENSE, sensitivity profiles can be calculated adaptively to track signal changes due to breathing or motion. In order to obtain accurate sensitivity profiles, images used to generate them must be alias-free. The amount of residual artifact after temporal filtering is related to a subject’s dynamic nature, number of temporal frames used for filtering, and filter bandwidth. Depending upon the application and subjects (e.g., task design in fMRI and children vs. adults), filtering parameters may require monitoring and adjustment.
Preibisch et al (5) used an EPI-SENSE sequence to investigate SNR vs. fMRI sensitivity at two spatial resolutions. At low spatial resolution, parallel imaging was used to reduce readout duration; surprisingly, SNR and fMRI sensitivity were hardly reduced. At high spatial resolution, parallel imaging was used to increase spatial resolution. Physiological noise was reduced and fMRI sensitivity was similar to that of the regular low-resolution single-shot acquisition. They concluded that parallel imaging has negligible effect on fMRI sensitivity despite decrease in SNR.
The lack of reduction in fMRI sensitivity with decreasing SNR in parallel imaging has been observed and explained by several authors. (5,6) fMRI sensitivity is inversely related to temporal signal instability originating from physiological and thermal noise. Parallel imaging increases thermal noise but not physiological noise. Thus, when physiological noise dominates, fMRI sensitivity is hardly affected. Therefore, as long as physiological noise is much greater than thermal noise, PI fMRI sensitivity is comparable to non-PI techniques. Nevertheless, a major incentive to use parallel imaging at high field is high spatial resolution; i.e., small voxel size. As voxel size decreases, thermal noise dominates. At present, the limiting factor in achieving high spatial resolution with parallel imaging (at high field and with sufficient fMRI sensitivity) is the fact that thermal noise leads to temporal signal instability, thereby degrading fMRI sensitivity.
Herein we present a PI technique for improving fMRI sensitivity at high spatial resolution when thermal noise predominates. We propose measurement of coil sensitivity by combining fully sampled data from a number of adjacent frames in an interleaved acquisition. In a two-shot interleaved acquisition, for example, k-space is fully sampled at every two repetition times (TR). Sensitivity profiles are then calculated from these full-resolution alias-free images and updated synchronously with image acquisition. This technique provides the convenience of self-calibrating sensitivity maps. In addition, we demonstrate that the technique provides inherent noise filtering and significantly enhances fMRI activation measurement. Improvement is especially significant at higher resolution and thinner slices where thermal noise exceeds physiological noise. For example, the technique detects activation at high spatial-resolutions (e.g. slice thickness of 2mm) better than previous methods. The technique results in higher temporal resolution and slightly more or comparable activated volumes compared to a standard two-shot reconstruction in fMRI experiments.
The technique may be incorporated into any interleaved acquisition. We demonstrate an implementation with two-shot interleaved spiral-in/out trajectory using an acceleration factor R = 2. The spiral-in/out trajectory is known to be effective in reducing signal drop-out from susceptibility-induced field gradients. (7) We compare the proposed sliding window PI technique with conventional (non-accelerated) two-shot reconstruction, as well as with the conventional PI reference frame method and an all-frame method. In the all-frame method, data acquired over the whole time-series is averaged together to calculate coil sensitivity. In all three comparisons, the proposed technique demonstrates enhanced BOLD activation. The improvement is most significant for high resolution, thin slice, and low SNR cases.
Sliding Window method
In the sliding window method, several adjacent two-shot images are selected for calculation of sensitivity profiles. These profiles are then applied exclusively to reconstruct data under that window. The minimum sliding window width that will provide fully-sampled data, with which to calculate sensitivity maps, is the acceleration factor R ; which is taken to be 2 here for all demonstrations and further discussion. The technique eliminates the temporal filtering step required by TSENSE. Since the sensitivity profiles are calculated from two-shot images, they are fully sampled, unaliased, and free of other post-processing artifacts. Moreover, these two-shot images are reconstructed from functional data itself; so no extra scans are needed. We shall show that our proposed coil sensitivity calculations reduce image noise when compared with all-frame or reference frame methods both in theory and in vivo, especially when thermal noise equals or exceeds physiological noise.
The original TSENSE diagram from Kellman (4) and a diagram of our proposed sliding window technique are shown in Figure 1 for comparison. The main difference between our technique and TSENSE is that TSENSE calculates its sensitivity maps by temporally lowpass filtering time frames of undersampled k-space data. Those time frames have alternating k-space undersampling patterns so that aliasing is modulated at the Nyquist frequency and then removed via lowpass filtering. The minimum number of time frames required for lowpass filtering is the number of alternating k-space sampling patterns (same as reduction factor R). In practice, the number of time frames used for lowpass filtering is generally greater than R for better separation of desired and aliasing frequency components. In our sliding window technique, sensitivity maps are calculated from multi-shot reconstructed images. The number of time frames used in our multi-shot reconstruction is equal to the number of alternating k-space sampling patterns (exactly R), so that the calibration maps are fully sampled, thus distortion free. The impact of averaging more time frames will be discussed later where we will find that reduction in thermal noise is best when R frames are used.
Figure 1
Figure 1
(a) TSENSE diagram adapted from Kellman (4). (b) Proposed sliding window technique (R=2) as comparison.
Effect of noise
We express the acquired k-space data as
equation M1
[1]
where m is the object image vector, p is a physiological noise vector, ε is a vector of thermal noise, S is a matrix of coil sensitivities, and E is the Fourier kernel matrix. A reconstructed image with estimated sensitivity maps Ŝ is expressed
equation M2
[2]
We show, in Appendix A, how an image reconstructed under the sliding window method can be expressed
equation M3
[3]
where [epsilon with overline] is a vector of thermal noise averaged over a number of consecutive temporal frames, and Q−1 is a normalization matrix to the estimate of coil sensitivity.
When ε=[epsilon with overline], the image estimate m is free of thermal noise. But that would occur only when R=1, which is not of interest since no readout acceleration is provided. When sliding window width is minimal (R=2), then [epsilon with overline] is the closest approximation to thermal noise ε achievable and so m will contain the least amount of noise. As the sliding window width becomes wider, [epsilon with overline] resembles ε to a lesser degree which results in m having more thermal noise. When the all-frame method is used, on the other hand, [epsilon with overline] is averaged to zero under a Gaussian noise assumption; thermal noise sampled in k-space is therefore propagated to m. Smoothing sensitivity maps would also result in thermal noise propagation to m.
The difference of the noise terms in Eq.[3] leads to a reduction of noise because thermal noise in ε is correlated with that in [epsilon with overline]; that is, the case of a sliding window method with a narrow window width. Otherwise, the subtraction would not necessarily lead to a noise reduction.
Data Acquisition
We demonstrate the technique with a two-shot spiral-in/out trajectory providing an acceleration factor of R=2. For all even time frames, the trajectory is rotated by 180 degrees so that k-space is fully sampled every two TRs. Experiments were performed with four slice thicknesses (2mm, 3mm, 4mm, and 5mm) in order to alter the relative contributions of thermal and physiological noise, and the number of slices gathered was 9, 7, 7, and 5 respectively to cover approximately the same volume of brain tissue. Slices were separated with a gap of 1mm. All experiments were performed on a GE 3T whole-body scanner (GE Signa, GE Healthcare, Waukesha, WI) equipped with a maximum gradient amplitude of 40 mT/m and a slew rate of 150 mT/m/ms. An 8-channel head coil (MRI Devices Corporation, Pewaukee, WI) was used for all image acquisitions. The TE was set to be minimal allowed by the spiral-in trajectory (35.4 ms), TR/α/matrix size/FOV = 2s/70°/128×128/20cm. The scan time was 248 seconds (data from the first 8 seconds discarded). T2-weighted fast spin-echo (FSE) scans were obtained for anatomic reference (TR/TE/ETL=4000ms/68ms/12). Six healthy subjects participated after giving informed consent in accordance with a protocol approved by the Stanford Institutional Review Board.
The fMRI task consisted of 6 cycles of an on/off block design having a period of 40s. During the on-block, subjects saw a checker board flashing at 8Hz. Subjects were told to stare at a fixation-cross to reduce eye movement during the off-block. To keep subjects’ attention throughout scans, a red cross appeared at random times; twice for 0.5s during each on-block. Subjects were instructed to press a button as soon as they saw the cross. This experiment was performed once for each slice thickness; the order of the slice thickness scans was randomized among subjects.
Subjects 1–3 also participated in noise measurement scans in separate sessions. Noise measurement scans are similar to functional scans above except that the flip angle is set to 0 so that only thermal noise is gathered. Only 3 slices and 64 time frames were collected. Two scans were performed: flip-angle = 70° and 0°. A fast T1 mapping scan (8) was utilized to obtain images for later gray matter segmentation; one for each slice thickness. Subjects were instructed to relax while no stimulation or task was given.
Data Analysis
Sensitivity profiles are calculated using the functional data itself. No extra calibration scans are needed; this reduces inconsistency between a reference frame and the functional data to be reconstructed. Standard spiral trajectory gridded-reconstruction is used to make two-shot images from each coil's data. (9) For every two TRs, in other words, two fully-sampled two-shot images (one spiral-in and one spiral-out) are reconstructed for each coil. The sensitivity profile of each coil is the ratio of a fully sampled two-shot image from that coil to the square root of pixel-wise sum of squares of all coil images. Separate sensitivity profiles are generated for spiral-in and spiral-out data. For the all-frame method, all fully sampled two-shot images of each coil over the entire time series are first averaged for sensitivity profile calculation. The resulting sensitivity profiles are used for conjugate gradient (CG-)SENSE reconstruction (10,11) of the entire time series. The number of iterations in CG-SENSE reconstruction is 7 regardless of sensitivity map calculation method. For the sliding window method (window width = 2 TRs), a new fully sampled two-shot image from each coil is formed after every TR. The sensitivity profile for each coil is then updated at every TR and used for CG-SENSE image reconstruction at that TR. The window width (number of TRs included in the average for SENSE map calculation) is assumed to be 2 TRs for the sliding window method unless otherwise noted. Concomitant field effects correction and navigator correction are first performed on the raw data. (12) The same CG-SENSE program (written in MATLAB, Mathworks R13, Natick, MA) is used for reconstruction regardless of sensitivity profile calculation methods; it generates one spiral-in image and one spiral-out image for every TR.
Data gathered for noise measurement (flip angle 70° and 0° with no task performed) is reconstructed with all-frame and sliding window methods described above. A gray matter mask is segmented manually from a T1 map for each slice and each slice thickness. This mask is used to extract gray matter pixel time-series from data at flip-angle=70°. The temporal standard deviation of each gray matter pixel is recorded and then averaged over three slices; this gives the total noise σ. The thermal noise σ0 is calculated by the same procedure using images with flip-angle = 0°, and the resulting standard deviation is multiplied by 1.527 to correct for the Rayleigh distribution of Gaussian noise in magnitude images. (13,14) Physiological noise σp is then calculated from (15)
equation M4
[4]
The method for calculating t-score is adapted from Lee et al. (16) (also see 1719). The steps involved are summarized here: each fMRI time series is first fitted to a second order polynomial for detrending. It is then correlated with sine and cosine functions at the fundamental task frequency (higher harmonics are negligible because the task blocks are short). The temporal phase of those voxels, with correlation coefficients exceeding a threshold of 0.2, is computed from the arctangent of sine and cosine correlations and averaged across each slice. The correlation is then projected onto this phase, per voxel, to provide a maximally in-phase correlation map. This model-free (Fourier) method eliminates bias from inaccuracy in assumed hemodynamic response function. Finally, a Fisher transform converts the correlation coefficient to t-score according to each voxel’s degrees of freedom. A sigma filter (12,20) is applied to these maps to cluster pixels in a 3×3 region, thereby reducing single-voxel false positives.
Figure 2 shows gray matter signal magnitude, physiological noise, and thermal noise averaged over three subjects using the all-frame and sliding window methods of sensitivity map calculation at each slice thickness. As expected, signal magnitude increases with slice thickness. Physiological noise, being proportional to signal magnitude, also increases with slice thickness while thermal noise is fairly constant. It is evident that the sliding window method effectively reduces thermal noise while signal magnitude and physiological noise are independent of sensitivity map calculation methods.
Figure 2
Figure 2
(top) Mean and standard deviation of signal magnitude, (middle) thermal noise, and (bottom) physiological noise in gray matter from three subjects. Solid lines represent results using the all-frame method. Dash lines are results using sliding window.
Figure 3 shows typical images obtained using all-frame and sliding window methods of sensitivity map calculation at 2mm and 5mm thickness. Images using the sliding window method are less noisy than those using all-frame. Standard deviation (SD) images, signal-to-fluctuation noise (SFNR) images, and activation maps all benefit when sliding window is used. This is more evident in thin slices (2mm) than in thick slices (5mm). These results are consistent for both spiral-in and spiral-out trajectories. The sliding window width default is 2 TRs. However, noise measurement data from subjects 1–3 were further reconstructed using the sliding window method with window widths of 4, 8, and 16 TRs. Figure 4 shows the standard deviation for the thermal noise σ0 vs. sliding window width. The noise reduction is only apparent for window width of 2 TRs. Beyond a window width of 2 TRs, sliding window reconstruction becomes the all-frame method in the limit. For 2mm to 5mm spiral-in and spiral-out images, Figure 5 shows ratios of activation volumes using the sliding window method to those using the all-frame method. The fat bars are averaged ratios of all six subjects. Standard deviations among subjects are indicated by the skinny error bars. For comparison, ratios of activation volumes using 2-shot reconstruction to those using all-frame are also shown. The activation threshold is set to p<0.005 in all cases.
Figure 3
Figure 3
Typical Spiral-in and Spiral-out images of 2mm and 5mm thickness. Four types of images are shown in each case: (a) single time frame images, (b) standard deviation (SD) images, (c) signal-to-fluctuation (SFNR) images, and (d) activation maps (p<0.005). (more ...)
Figure 4
Figure 4
Thermal noise using various widths of sliding window from Subjects 1–3 over 3 slices each with 64 frames.
Figure 5
Figure 5
Ratios of numbers of activated pixels from using sliding window to those from using all-frame. Also shown are the ratios using 2-shot reconstruction to all-frame (p<0.005). Advantage of using sliding window over all-frame is more apparent in thinner (more ...)
Denoising with Sliding Window
We hypothesize that the sliding window reconstruction removes image noise because the noise information is retained in the corresponding sensitivity maps. To demonstrate that the denoising effect (shown in Figure 3Figure 5) derives from incorporation of temporally localized fluctuations in sensitivity maps, 2mm slice thickness data from subject 2 was again reconstructed with sliding window methods (width = 2 TRs) but with a window positioned 2 TRs in advance of its optimal position (concurrent with the time frame being reconstructed). When using sensitivity maps from temporally advanced data, reconstructed images are visibly noisier; as illustrated in Figure 6.
Figure 6
Figure 6
(top) Spiral-in and (bottom) Spiral-out images of 2mm thickness from volunteer 2. Images on the left use correct (current) sliding window sensitivity maps. Those on the right use sensitivity maps generated from images 2TRs ahead. The increase in noise (more ...)
As further demonstration of the hypothesis that sliding window sensitivity maps reduce image noise after CG-SENSE reconstruction, the sensitivity maps used in sliding window reconstruction were spatially smoothed by spline interpolation for Subject 4. Figure 7 demonstrates unsmoothed and smoothed sensitivity maps from Coil 1 as well as images reconstructed using both types of sensitivity maps. We expect that smoothing sensitivity maps would result in noisier images than without smoothing. It is apparent that the observed increase in image noise, particularly in Figure 7(c–d), is due to sensitivity map smoothing. Hence, smoothing the sensitivity maps is not desirable for the sliding window method. But that is not generally true for sensitivity maps corrupted by noise less correlated with images under reconstruction, as occurs with the conventional all-frame or reference frame methods.
Figure 7
Figure 7
(top) Spiral-in and (bottom) Spiral-out images with 2mm thickness from subject 4. (a) Images reconstructed with sliding window sensitivity maps of which those corresponding to coil 1 are shown in (e). (f) Sliding window sensitivity maps are spatially (more ...)
Comparison with Reference Frame Method
As mentioned earlier, a widely used technique for gathering sensitivity profiles is to obtain a full resolution scan prior to engagement of the subsampling methods. We define our reference frame method as averaging four image time frames from the beginning of a time series and then smoothing with spline interpolation. Figure 8 shows ratios of numbers of activated pixels using the reference frame method to those using the all-frame technique for 2mm and 5mm spiral-in and spiral-out images. The bars are the averaged ratios of all six subjects. Standard deviations among them are indicated with lines drawn above each bar. For comparison, ratios of numbers of activated pixels using sliding window to those using all-frame are also shown. Activation threshold is set to p<0.005 in all cases. As shown in Figure 8, the sliding window technique performs significantly better than both reference frame and all-frame techniques, while the latter two perform comparably. The improvement was especially more significant at the 2mm thickness.
Figure 8
Figure 8
Ratios of numbers of activated pixels from using sliding window to those from using all-frame. Also shown are the ratios from using reference frame to all-frame (p<0.005). Reference frame performance is very similar to that of the all-frame technique (more ...)
The SNR for each subject’s reconstrcutions using sliding window, all-frame, and reference frame methods was calculated. Signal measurements were made by averaging magnitude intensities within a rectangular ROI (region of interest) inside the brain at each time frame. For a given slice, the same rectangular ROI, which included gray matter, white matter, and CSF were used in sliding window, all-frame, and reference frame reconstructed data. A similar rectangular region drawn outside the brain was used to estimate the noise from the standard deviation among those pixels at each time frame. Time series SNR was calculated by dividing the signal mean by noise standard deviation. Averaged mean and standard deviation of this time series SNR over all slices are listed in Table 1. Sliding window gives 15% to 70% increase in SNR over all-frame or reference frame methods. SNR is generally higher at thinner slice thickness for which thermal noise predominates over physiological noise.
Table 1
Table 1
SNR measurements of each subject using sliding window (width = 2 TRs), all-frame, and reference scan methods. Numbers inside parentheses indicate standard deviations over all slices and time frames.
Sensitivity Maps from Off Blocks and On-Off Blocks
To assess impact of signal intensity change (due to the task paradigm) on sensitivity maps and consequent image quality, two additional methods of generating sensitivity maps were studied:
In the first method, sensitivity maps were generated by averaging full resolution images only during the off blocks (no task) of task paradigm. These sensitivity maps are used for CG-SENSE reconstruction of the entire time series data. This is referred to as the off blocks method.
The second method is referred to the on-off blocks method; it makes two distinct sets of sensitivity maps: one for on blocks (task) and another for off blocks. The on blocks sensitivity maps are made from full resolution images gathered during task, and they are subsequently applied to CG-SENSE reconstruction of on blocks images only. The off blocks sensitivity maps calculated are the same as those generated from the off blocks method but are applied only to reconstruct off blocks images. (One image immediately adjacent to block boundaries is excluded in the sensitivity map calculation to avoid reaction delays and to ensure proper comparison.)
Figure 9Figure 11 show SD images, SFNR images, and correlation coefficient maps using off blocks and on-off blocks methods as well as those from all-frame, reference frame, and sliding window methods. From these figures, it is clear that neither the off blocks or the on-off blocks method provide lower image noise or higher SFNR values than the sliding window method.
Figure 9
Figure 9
Standard deviation images from subject 3 reconstructed using sensitivity maps from (a) off-blocks, (b) on-off blocks, (c) all-frame, (d) reference frame, and (e) sliding window. All images are scaled to the same level for display. Sliding window sensitivity (more ...)
Figure 11
Figure 11
Activated pixels (p<0.005) overlaid on anatomic images corresponding to those in Figure 9. Number of activated pixels are indicated next to each image. (a) Off-blocks, (b) on-off blocks, (c) all-frame, (d) reference frame, and (e) sliding window (more ...)
Comparison with Theory
Equation 3 implies that reduction of thermal noise in the reconstructed time frame occurs when thermal noise ε in a given time frame is close to that in the sensitivity map, [epsilon with overline] (Appendix A). Intuitively, peak correlation between thermal noise in sensitivity map and thermal noise of the raw data is inversely proportional to the width of sliding window, i.e., narrower window width corresponds to higher peak correlation. It is shown in Appendix B that the correlation coefficient of raw data thermal noise and sensitivity map thermal noise is approximately equation M5, where W is width of the sliding window (in number of TRs).
To quantitatively verify our theory, we gathered data containing only thermal noise in a separate experiment with echo planar imaging and gradient recalled echo (GRE EPI) (matrix size = 64 ×64, FOV=50cm, 8-channel head coil, no ramp sampling). The reason for using EPI instead of spiral trajectories is: with EPI we were able to acquire the fully sampled noise at each time frame and conveniently generate undersampled data by throwing away odd or even k-space lines. Although we were also able to collect a fully sampled k-space with a single-shot spiral, we were not able to generate undersampled data easily on a two-interleaf spiral without using data interpolation (which might further complicate the noise statistics).
In this GRE EPI experiment, a total of 60 time frames were acquired with the RF transmitter turned off. The first four time frames were discarded. Data from each coil at each time frame was reconstructed by Fourier transform. For each time frame, one thermal noise image was calculated by combining each coil image using square root of sum of squares. These images contain raw data noise ε. Sensitivity maps with window width of 2 TRs were formed using data from odd k-space lines of time frame 1 and data from even k-space lines of time frame 2. Sensitivity maps with window width of 4 TRs were formed using averaged data from odd k-space lines of time frame 1 and 3 and averaged data from even k-space line of time frame 2 and 4. Twenty eight sets of sensitivity maps (widths = 2, 4, 6, …,56 TRs) were formed. Temporal correlation coefficients were calculated between each set of sensitivity map and thermal noise images, as illustrated in Figure 12a. Correlation coefficients between the sensitivity map with a width of 2 (calculated using time frames 1 and 2) and the noise images at time frames 1 and 2 are approximately 0.5, whereas correlation coefficients between this sensitivity map and later time frames are close to 0. Correlation coefficients between the sensitivity map with a width of 4 (using time frames 1–4) and the noise images at time frames 1–4 are about 0.35. Figure 12b plots the peak correlation coefficients as a function of window width. The peak correlation coefficients for each window width are calculated as the mean coefficient over the entire plateau region of each corresponding curve shown in Figure 12a. The number of points used for calculating the peak coefficient is the same as sliding-window width. The theoretically predicted peak correlation coefficients are also plotted in Figure 12b, which agrees will with the experimentally measured values.
Figure 12
Figure 12
(a) Correlation coefficient, over time, between thermal noise data and sensitivity maps of different sliding window widths. (b) Averages of the plateau of each curve in (a) are plotted as circles. These results agree well with theory (square). First time (more ...)
Figure 12b also plots the noise standard deviation of CG-SENSE reconstructed images as a function of window width. A total of 28 images are reconstructed using the same data set (odd k-space lines of the first time frame) but with sensitivity maps of various sliding-window widths. Since k-space data is acquired on a Cartesian grid, reconstruction is carried out using simple least-squares matrix inversion. Standard deviation of each image is calculated and normalized by that using the widest window width. As expected, standard deviation decreases with decreasing window width. This behavior agrees with observations from the fMRI study shown in Figure 4.
We have shown, both experimentally and theoretically, that thermal noise is significantly reduced by using sliding window sensitivity maps for CG-SENSE reconstruction in fMRI. This is especially important when imaging at high resolution where thermal noise predominates. With small voxel size, the technique reveals activation from visual stimulation where other conventional sensitivity calculation techniques falter. With a slice thickness of 2mm, for example, the proposed technique results in improvement by a factor of 4 when compared to conventional PI techniques (Figure 5 and Figure 8). Further, activation volume is statistically equivalent to conventional 2-shot reconstruction with half the temporal resolution; which provides an incentive to perform fMRI with parallel imaging (CG-SENSE). Use of parallel imaging for fMRI has been debated in the literature with regard to loss of SNR and activation detectability. (3,2123) Utilizing the sliding window CG-SENSE approach, temporal resolution is enhanced and acquisition readout time is reduced (less off-resonance effects) with little loss in activation detectability.
Admitting noise into sensitivity maps appears to contradict conventional wisdom that sensitivity maps should be smooth since receiver coil sensitivity is slowing varying in space. It is true that coil sensitivity is inherently smooth and noiseless. That coil sensitivity which is measured by the sliding window technique not only provides an estimation of intrinsic coil sensitivity, but it also provides an estimate of thermal noise. Once this thermal noise estimate is obtained, the noise in the final image can be reduced by using these noisy sensitivity profiles in the image reconstruction process. By reducing thermal noise while maintaining physiological fluctuation induced by functional activation, we are able to improve the fMRI detection. In other words, receiving data using multiple coils simultaneously provides a means to separate physiological noise and thermal noise. The ability to separate these two noise types should increase with increasing number of coils. Further, by improving temporal resolution with parallel imaging, our technique also improves activation detection when compared to standard multi-shot technique.
Since the sliding window method uses the current and adjacent time frames for calculating sensitivity maps, which are then applied to CG-SENSE reconstruction, the temporal features of the time series can be slightly smeared from the increased autocorrelation. In particular, degrees of freedom among time frames are reduced. We calculated degrees of freedom for each pixel using a first order autoregressive model, proposed by Worsley et al. (24) and Kruggel et al. (17). Degrees of freedom are compared in Figure 13 using data (subject 1, 2mm spiral-out, slices 1–9) reconstructed via sliding window, all-frame, and reference frame methods. Histograms show distribution of degrees of freedom using these reconstruction methods. Increase in correlation and reduction in degrees of freedom using the sliding window method is small, and cause little adverse effects in activation statistics. That is because the t-distribution approaches a normal distribution, when degrees of freedom are greater than 30, for which degrees of freedom are irrelevant to the critical value for a p-statistic. (25,26) For example, for probability of error p<0.005, degrees of freedom 115 corresponds to critical t-score of 2.619, while degrees of freedom 90 corresponds to critical t-score of 2.632. Applying those t-scores to volunteer 3, for example (2mm spiral-out sliding window), the number of activated pixels is 136 and 133, respectively. For reference frame method, the number of activated pixels is 31 and 30 respectively. For all-frame method, the number of activated pixels is 31 in both cases.
Figure 13
Figure 13
Histograms showing degrees of freedom in a pixel by pixel basis using same data reconstructed by all-frame and sliding window (width = 2 TRs) methods. Sliding window method slightly reduces the degrees of freedom in time series than the all-frame method. (more ...)
Though demonstrated with fMRI, we expect this proposed sliding window technique to directly benefit other dynamic SENSE applications such as cardiac and flow imaging. It is also likely to reduce motion artifacts when scanning young children if motion is slow relative to sliding window widths, since sensitivity maps are updated over time and follow the motion.
Figure 10
Figure 10
SFNR images corresponding to those in Figure 9. (a) Off-blocks, (b) on-off blocks, (c) all-frame, (d) reference frame, and (e) sliding window sensitivity maps are used. All images are scaled to the same level for display. Sliding window sensitivity maps (more ...)
Acknowledgement
This work is supported by NIH RR09784, NIH-1K99NS057943-01, Lucas Foundation, and GEMS. The authors would like to thank Drs. Roland Bammer, Peter Kellman, Gunnar Krueger, and Tie-Qiang Li for helpful discussions.
Appendix A
Data received by the ith coil at a k-space location is
equation M6
[A1]
where r represents image index, and
d i(k):Sampled k-space data of ith coil
M(r):Original image
p(r):Physiologic noise
εi(k):Thermal noise
si(r):Sensitivity of ith coil
With Κ k-space samples and N object samples, define the following vectors and matrices with corresponding sizes:
equation M7
[A2]
equation M8
[A3]
equation M9
[A4]
equation M10
[A5]
equation M11
[A6]
Matrix E0 is a subsampled (fat) 2D DFT matrix constituted by entries ejK1r1 = ej(kx1x1+ky1y1); it is some phase of a complete DFT matrix. In the case of reduction factor R, the modulo frame index is defined l[set membership]{1,2…R}. For reconstruction of any particular frame, one of R phases of the DFT matrix is used depending on l. Matrix E0 may therefore be considered constant in all the following equations.
equation M12
[A7]
For a number c of coils (c ≥ 1), we have
equation M13
[A8]
equation M14
[A9]
equation M15
[A10]
equation M16
[A11]
The acquired data can be written as
equation M17
[A12]
Estimated coil sensitivity ŝi(r) is the ratio of the 2-shot reconstructed image from the ith coil to the square root of sum of squares of 2-shot images from all coils:
equation M18
[A13]
equation M19
[A14]
equation M20
[A15]
Here, [p with overline](r) and ei(r) are the noise averaged over a number of M consecutive temporal frames. Without frame averaging the corresponding variables are denoted
equation M21
[A15a]
When the estimated sensitivity maps are used in the image reconstruction, the sensitivity matrix becomes,
equation M22
[A16]
where X = e(M + P)−1 is assumed relatively small because e is small compared to M + P.
equation M23
[A17]
equation M24
[A18]
equation M25
[A19]
equation M26
[A20]
equation M27
[A21]
Using the estimated coil sensitivity, image is reconstructed as
equation M28
[A22]
Under the assumptions that X is small and using first-order Taylor series approximation for an inverse matrix, for
equation M29
[A23]
equation M30
[A24]
Using first-order Taylor series approximation for matrix inversion and assuming X is small,
equation M31
[A25]
equation M32
[A26]
equation M33
[A27]
Since X and ε are small and Ee = ε (or e = EHε),
equation M34
[A28]
Assuming physiological noise is slowing varying under width of sliding window, Pp. (In the limit of window width = 1 frame, P = p and exact cancellation occurs.)
equation M35
[A29]
Clearly, resulting noise ε[epsilon with overline] in the reconstructed image is reduced if thermal noise ε at a given frame is similar to the noise averaged over a number of consecutive frames [epsilon with overline].
Appendix B
Thermal noise in EPI raw data of a time series can be expressed by independent random variables O1i,E1i,O2i,E2iONi,ENi, where O and E denote odd and even k-space lines of time frames 1 through N of the ith coil.
The sensitivity map formed using time frames 1 and 2 is
equation M36
The bracket notation represents filling in odd and even k-space lines by its contents, and F denotes Fourier Transform. Since raw data is entirely made up with thermal noise, then Ŝi,w = 2 = e where e is thermal noise in the image domain, and [epsilon with overline] is thermal noise in k-space.
For sliding window width = 4 (using time frames 1–4), the sensitivity map is
equation M37
Thermal noise ε in equation [A29] of the first time frame is O1i or E1i depending upon which half of k-space is sampled. Problems arise in obtaining thermal noise e, since applying a Fourier Transform to O1i or E1i directly will result in aliasing. This leads to the reason for fully sampling k-space in this experiment: the Fourier Transform can be applied to data of each time frame without aliasing; e.g., first image im = F {[O1i, E1i]} with thermal noise e embedded. We then calculate correlation between im1 and Ŝi,w. With eight coils, sensitivity Ŝi,w has Rician distribution, and im1 is Rayleigh distributed.
For sliding window width = 2, the correlation coefficient is equation M38
For window width = 4, equation M39
In general, for window width = W, equation M40
1. Uurbil K, Adriany G, Andersen P, Chen W, Gruetter R, Hu X, Merkle H, Kim D-S, Kim S-G, Strupp J, Zhu XH, Ogawa S. Magnetic Resonance Studies of Brain Function and Neurochemistry. Annual Review of Biomedical Engineering. 2000;2(1):633–660. [PubMed]
2. Schmidt CF, Degonda N, Luechinger R, Henke K, Boesiger P. Sensitivity-encoded (SENSE) echo planar fMRI at 3T in the medial temporal lobe. Neuroimage. 2005;25(2):625–641. [PubMed]
3. Weiger M, Pruessmann KP, Osterbauer R, Bornert P, Boesiger P, Jezzard P. Sensitivity-encoded single-shot spiral imaging for reduced susceptibility artifacts in BOLD fMRI. Magn Reson Med. 2002;48(5):860–866. [PubMed]
4. Kellman P, Epstein FH, McVeigh ER. Adaptive sensitivity encoding incorporating temporal filtering (TSENSE) Magn Reson Med. 2001;45(5):846–852. [PubMed]
5. Preibisch C, Pilatus U, Bunke J, Hoogenraad F, Zanella F, Lanfermann H. Functional MRI using sensitivity-encoded echo planar imaging (SENSE-EPI) Neuroimage. 2003;19(2 Pt 1):412–421. [PubMed]
6. de Zwart JA, van Gelderen PPK, Duyn JH. Application of sensitivity-encoded echo-planar imaging for blood oxygen level-dependent functional brain imaging. Magnetic Resonance in Medicine. 2002;48(6):1011–1020. [PubMed]
7. Glover GH, Law CS. Spiral-in/out BOLD fMRI for increased SNR and reduced susceptibility artifacts. Magn Reson Med. 2001;46(3):515–522. [PubMed]
8. Hsu JJ, Glover GH. Rapid MRI method for mapping the longitudinal relaxation time. J Magn Reson. 2006;181(1):98–106. [PubMed]
9. Jackson JI, Meyer CH, Nishimura DG, Macovski A. Selection of a convolution function for Fourier inversion using gridding. IEEE Trans Med Imaging. 1991;10:473–478. [PubMed]
10. Pruessmann KP, Weiger M, Börnert P, Boesiger P. Advances in sensitivity encoding with arbitrary k-space trajectories. Magnetic Resonance in Medicine. 2001;46(4):638–651. [PubMed]
11. Liu C, Moseley ME, Bammer R. Simultaneous phase correction and SENSE reconstruction for navigated multi-shot DWI with non-cartesian k-space sampling. Magn Reson Med. 2005;54(6):1412–1422. [PubMed]
12. Glover GH, Lai S. Self-navigated spiral fMRI: interleaved versus single-shot. Magnetic Resonance in Medicine. 1998;39(3):361–368. [PubMed]
13. Gudbjartsson H, Patz S. The Rician distribution of noisy MRI data. Magn Reson Med. 1995;34(6):910–914. [PMC free article] [PubMed]
14. Kellman P, McVeigh ER. Image reconstruction in SNR units: a general method for SNR measurement. Magn Reson Med. 2005;54(6):1439–1447. [PMC free article] [PubMed]
15. Kruger G, Glover GH. Physiological noise in oxygenation-sensitive magnetic resonance imaging. Magn Reson Med. 2001;46(4):631–637. [PubMed]
16. Lee AT, Glover GH. Discrimination of large venous vessels in time-course spiral Blood Oxygen Dependent Magnetic resonance functional neuroimaging. Magn Reson Med. 1995;33:745–754. [PubMed]
17. Kruggel F, Pelegrini-Issac M, Benali H. Estimating the effective degrees of freedom in univariate multiple regression analysis. Medical Image Analysis. 2002;6(1):63–75. [PubMed]
18. Bandettini PA, Jesmanowicz A, Wong EC, Hyde JS. Processing strategies for time-course data sets in functional mri of the human brain. Magnetic Resonance in Medicine. 1993;30(2):161–173. [PubMed]
19. Miller I, Freund JE. Probability and Statistics for Engineers. Prentice-Hall; 1977.
20. Pratt W. Digital Image Processing. New York: John Wiley & Sons, Inc.; 1991.
21. de Zwart JA, van Gelderen P, Golay X, Ikonomidou VN, Duyn JH. Accelerated parallel imaging for functional imaging of the human brain. NMR Biomed. 2006;19(3):342–351. [PubMed]
22. Golay X, de Zwart JA, Ho YC, Sitoh YY. Parallel imaging techniques in functional MRI. Top Magn Reson Imaging. 2004;15(4):255–265. [PubMed]
23. Lutcke H, Merboldt KD, Frahm J. The cost of parallel imaging in functional MRI of the human brain. Magn Reson Imaging. 2006;24(1):1–5. [PubMed]
24. Worsley KJ, Friston KJ. Analysis of fMRI time-series revisited--again [comment] Neuroimage. 1995;2(3):173–181. [PubMed]
25. Dowdy SW. Statistics for Research. Wiley-Interscience; 1991.
26. Huettel SA, Song AW, McCarthy G. Functional Magnetic Resonance Imaging. Sunderland, MA: Sinauer Associates; 2004. p. 492.