PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Neuroimage. Author manuscript; available in PMC 2010 August 1.
Published in final edited form as:
PMCID: PMC2782710
NIHMSID: NIHMS106675

Superresolution parallel magnetic resonance imaging: Application to functional and spectroscopic imaging

Abstract

Standard parallel magnetic resonance imaging (MRI) techniques suffer from residual aliasing artifacts when the coil sensitivities vary within the image voxel. In this work, a parallel MRI approach known as Superresolution SENSE (SURE-SENSE) is presented in which acceleration is performed by acquiring only the central region of k-space instead of increasing the sampling distance over the complete k-space matrix and reconstruction is explicitly based on intra-voxel coil sensitivity variation. In SURE-SENSE, parallel MRI reconstruction is formulated as a superresolution imaging problem where a collection of low resolution images acquired with multiple receiver coils are combined into a single image with higher spatial resolution using coil sensitivities acquired with high spatial resolution. The effective acceleration of conventional gradient encoding is given by the gain in spatial resolution, which is dictated by the degree of variation of the different coil sensitivity profiles within the low resolution image voxel. Since SURE-SENSE is an ill-posed inverse problem, Tikhonov regularization is employed to control noise amplification. Unlike standard SENSE, for which acceleration is constrained to the phase-encoding dimension/s, SURE-SENSE allows acceleration along all encoding directions — for example, two-dimensional acceleration of a 2D echo-planar acquisition. SURE-SENSE is particularly suitable for low spatial resolution imaging modalities such as spectroscopic imaging and functional imaging with high temporal resolution. Application to echo-planar functional and spectroscopic imaging in human brain is presented using two-dimensional acceleration with a 32-channel receiver coil.

Keywords: Parallel imaging, Superresolution, SENSE, fMRI, Spectroscopic imaging

Introduction

Magnetic resonance imaging (MRI) methods involve imaging objects with high spatial frequency content in a limited amount of time. However, information over only a limited k-space range is usually acquired in practice due to SNR and time constraints. For example, in functional MRI (fMRI) (Belliveau et al., 1991) k-space coverage is traded off for increased temporal resolution. In MR spectroscopic imaging (MRSI) (Brown et al., 1982), which is constrained by relatively low SNR, k-space coverage is sacrificed to achieve an adequate SNR within a feasible acquisition time. The lack of high k-space information leads to limited spatial resolution and Gibbs ringing when the Fourier transform is directly applied to reconstruct the image. Constrained image reconstruction techniques using prior information (Liang et al., 1992) have been proposed to achieve superresolution image reconstruction, i.e. to estimate high k-space values without actually measuring them. For example, the finite spatial support of an image can be used to perform extrapolation of k-space at expense of SNR loss. However, this method performs well only at positions close to the periphery of the object being imaged (Plevritis and Macovski, 1995). For experiments with temporal repetitions such as fMRI and MRSI; k-space substitution (Jones et al., 1993), also known as the key-hole method, was proposed to fill the missing high k-space values of the series of low resolution acquisitions using a high resolution reference. However, this method is vulnerable to artifacts due to inconsistencies between the reference and the actual acquisition. An improvement of this approach, known as generalized series reconstruction (Liang and Lauterbur, 1991), forms a parametric model using the high resolution reference to fit the series of low resolution acquisitions and thus reduce the effect of data replacement inconsistencies. Alternatively, superresolution reconstruction can be performed by combination of several low resolution images acquired with sub-pixel differences (Elad and Feuer, 1997). This method is well developed for picture and video applications and was employed before in MRI by applying a sub-pixel spatial shift to each of the low resolution acquisitions (Peled and Yeshurun, 2001). However, its application is very limited since a spatial shift is equivalent to a linear phase modulation in k-space, which does not represent new information to increase the k-space coverage of the acquisition.

Parallel MRI (Sodickson and Manning,1997; Pruessmann et al.,1999) has been introduced as a method to accelerate the sequential gradient-encoding process by reconstructing an image from fewer acquired k-space points using multiple receiver coils with different spatially varying sensitivities. The standard strategy for acceleration is to reduce the density of k-space sampling beyond the Nyquist limit while maintaining the k-space extension in order to preserve the spatial resolution of the fully-sampled acquisition. The rationale for this sub-sampling scheme is that the coil sensitivities are spatially smooth and retrieve k-space information only from the neighborhood of the actual gradient-encoding point. However, any arbitrary k-space sub-sampling pattern can in principle be employed at the expense of increasing the computational cost and decreasing the stability of the inverse reconstruction, i.e. increasing the nominal noise amplification in the reconstruction (g-factor) (Pruessmann et al., 1999; Sodickson and McKenzie, 2001). On the other hand, standard parallel MRI performed at a spatial resolution that presents intra-voxel coil sensitivity variation suffers from residual aliasing artifacts which are depicted as spurious side lobes around the aliasing positions in the reconstructed point spread function (PSF) (Zhao et al., 2005). The minimum-norm SENSE (MN-SENSE) technique was proposed to remove the residual aliasing artifact by performing an intra-voxel reconstruction of the PSF using coil sensitivities acquired with higher spatial resolution (Sanchez-Gonzalez et al., 2006). However, while this method improves the spatial distinctiveness of image voxels, it does not increase the number of voxels and hence the underlying spatial resolution.

Receiver arrays with a large number of small coils tend to have rapidly varying coil sensitivity profiles, and therefore offer the promise of high accelerations for parallel imaging. However, standard Sensitivity Encoding (SENSE) reconstruction with many-element arrays is exposed to residual aliasing artifacts due to potential intra-voxel coil sensitivity variations. On the other hand, many-element arrays open the door for other k-space sub-sampling patterns that might not be feasible with few-element arrays. For example, highly accelerated parallel imaging using only one gradient-encoding step was proposed in the Single Echo Acquisition (SEA) technique with a 64-channel planar array (McDougall and Wright, 2005) and in the Inverse Imaging (InI) technique with a 90-channel helmet array for human brain fMRI (Lin et al., 2006, 2008). However, the price to pay for these extreme levels of acceleration is reconstruction with low spatial resolution as dictated by the degree of variation of the coil sensitivity maps.

The current work presents a novel method for parallel MRI in which acceleration is performed by acquiring only the central region of k-space instead of increasing the sampling distance over the complete k-space matrix. The proposed method, termed Super-resolution SENSE (SURE-SENSE), increases the spatial resolution of the fully-sampled low resolution acquisition using coil sensitivities acquired with higher resolution. Regularization of the ill-conditioned inverse reconstruction is performed to control noise amplification due to the relatively large weights required to reconstruct high k-space values from low resolution data. The attainable increase in spatial resolution is determined by the degree of variation of the coil sensitivities within the acquired image voxel. Application of SURE-SENSE to increase the spatio-temporal resolution of fMRI is presented. Unlike standard SENSE, which is constrained to acceleration of phase-encoding dimensions, SURE-SENSE in this case allows two-dimensional acceleration of the echo-planar acquisition. Application to MRSI of human brain is presented as a method to reduce lipid contamination and to enhance the spatial resolution of the metabolite maps in two dimensions.

Methods

Superresolution SENSE (SURE-SENSE)

The goal of superresolution SENSE is to reconstruct a single image with higher resolution from fully-sampled low resolution images acquired with multiple receiver coils using high resolution coil sensitivity maps (Fig. 1). Since the image acquired by each coil is weighted by the corresponding spatial sensitivity of the coil, superresolution reconstruction is feasible if the different sensitivities are varying within the low resolution image voxel.

Fig. 1
Conceptual illustration of the superresolution parallel MRI technique. The low resolution multi-coil data are combined into a single image with higher resolution by performing an intra-voxel reconstruction with high spatial resolution coil sensitivity ...

SURE-SENSE is formulated in the image domain following the generalized parallel imaging model for arbitrary sub-sampling trajectories (Sodickson and McKenzie, 2001) considering that image data are acquired from a central k-space region and coil sensitivity data are acquired from an extended k-space region. Both data sets are acquired on a grid given by the Nyquist sampling distance (Δk=1/FOV, where FOV is the field of view). The signal acquired by each coil can be represented as:

Sl(k)=rρ(r)Cl(r)ej2πk·rdr,l=1,2,,Nc,
(1)

where r is the position vector, k=γ2π0tG(τ)dτ is the k-space vector determined by the gradient vector G(t), ρ(r) is the object function, cl (r) is the complex-valued spatially varying coil sensitivity and Nc is the number of coils. Considering the acquisition of Nk image data points and Ns sensitivity points (Ns = RNk, where R is the overall sampling reduction factor), a discretized version of Eq. (1) in matrix form is given by:

FNksl=ΠFNsClρ,
(2)

where Fn (n × n) is the spatial discrete Fourier transform (DFT) matrix, sl (Nk ×1) is the low resolution image vector for the l-th coil, Cl (Ns × Ns) is a diagonal matrix containing the l-th coil sensitivity values along the diagonal, and ρ(Ns ×1) is the target object function at high spatial resolution. Π(Nk × Ns) is the low-pass k-space filter operator, where the element π(i,j) is equal to 1 if the k-space position with index j is sampled and equal to 0 otherwise. The encoding equation for the l-th coil in the image domain can be expressed as:

sl=FNk1ΠFNsClρ=Elρ.
(3)

Note that FNk1ΠFNs represents the low-pass k-space filter in the image domain. The complete encoding equation is obtained by concatenating the individual encoding equations:

s=Eρ,s=[s1sNc],E=[E1ENc],
(4)

where s is the multi-coil image vector at low resolution (NkNc ×1) and E is the encoding matrix (NkNc × Ns). Noise correlation between coils is removed by pre-whitening the image vector and the encoding matrix using the noise covariance matrix estimated from a noise-only acquisition (Pruessmann et al., 2001; Lin et al., 2004). Pre-whitening is performed by xw=ψ12x, where Ψ (Nc ×Nc) is the noise covariance matrix for the array coil and xw (Nc ×1) is the multi-coil data. Ψ was estimated using a sample average estimate from a noise-only acquisition (nt) switching off the RF excitation:

ψ^=1Ntt=1Nt(ntn¯)(ntn¯)H,
(5)

where Nt is the number of time points and [n with macron] (Nc ×1) is the average over time of nt.

The proposed k-space sub-sampling pattern, which reduces the extent of sampled k-space while maintaining the Nyquist sampling distance, allows for acceleration along the readout dimension. Two-dimensional acceleration of imaging sequences with two spatial dimensions will therefore be feasible with SURE-SENSE, unlike with standard SENSE where the acceleration is limited to the phase-encoding dimension. Two-dimensional acceleration provides better conditioning of the inverse problem and thus allows for higher acceleration factors than one-dimensional acceleration for the same overall acceleration factor when an array with two-dimensional coil sensitivity encoding is employed (Weiger et al., 2002).

The inverse of the encoding equation will provide an image reconstructed onto the superresolution grid, where the acquired low resolution data are fitted to delta functions in the high resolution grid using the high resolution coil sensitivities. Direct inversion will result in large noise amplification since the encoding matrix for SURE-SENSE is intrinsically ill-conditioned as compared with standard SENSE due to the lower coil sensitivity variation within the low resolution voxel than across aliased voxels. Tikhonov regularization is employed to control the noise amplification in the reconstruction (g-factor) (Lin et al., 2004, 2005). The least-squares solution using Tikhonov regularization with diagonal weighting can be represented as:

ρ^=argminρ{||Eρs||22+λ2||ρ||22},
(6)

where λ is the regularization parameter. Tikhonov regularization constrains the power of the solution ( ||ρ||22) thus controlling noise amplification while attenuating solution components with low singular values compared to λ (Hansen,1998). The Tikhonov weighting function for the i-th singular value σi is given by wi=σi2σi2+λ2, which presents a smooth roll-off behavior along the singular value spectrum instead of the sharp cut-off imposed by other techniques such as the truncated singular value decomposition (TSVD). The regularization parameter (λ2) was set to the average power of the high resolution reference used for sensitivity calibration to attenuate components with a low squared singular value compared to the average power of the reference. For the SURE-SENSE encoding matrix, low singular values represent high spatial frequency components, and therefore the regularization approach trades off SNR and spatial resolution. The nominal gain in spatial resolution is given by inversion of the encoding matrix without regularization. The effective gain in spatial resolution is dictated by the degree of regularization necessary to achieve an adequate SNR.

SURE-SENSE reconstruction requires the inversion of the complete encoding matrix. 1D-SURE is implemented using a line-by-line matrix inversion approach. 2D-SURE is implemented using a conjugate gradient (CG) solution with pre-conditioning as in the case of non-Cartesian SENSE (Pruessmann et al., 2001). For SURE-SENSE, density correction and regridding are not performed since the reconstruction lies on a Cartesian grid. Considering the original normal equations from Eq. (6) (EHE2I)ρ =EHs, the CG iterations are applied to the following transformed system:

P(EHE+λ2I)Pbi=PEHs,
(7)

where the elements of the diagonal pre-conditioning matrix P are given by pi,i=(l=1Nccl(ri)2+λ)12 and bi is the partial result for the i-th iteration. The final result after Ni iterations is then given by [rho with circumflex] = PbNi.

Spatial resolution analysis

The spatial resolution of the reconstruction was analyzed using the full-width at half-maximum (FWHM) of the point spread function (PSF). The effective gain in spatial resolution is defined as K=FWHM-DFT/FWHM-SURE, where FWHM-DFT is the FWHM of the conventional DFT reconstruction of the low resolution data with k-space zero-filling and FWHM-SURE is the FWHM of the SURE-SENSE reconstruction. The nominal gain in spatial resolution is given by K using the FWHM of SURE-SENSE without regularization. The PSF for each spatial position r was obtained by reconstructing the low resolution representation of a single source point located at r. The source point is modeled as a delta function at the corresponding voxel position, which is multiplied by the high resolution coil sensitivities and passed trough the low-pass filter. The PSF is given by the resulting SURE-SENSE reconstruction of the low resolution source point.

Noise amplification analysis

Noise amplification in the inverse reconstruction was assessed using the g-factor formalism (Pruessmann et al., 1999). For 1D-SURE, the analytical g-factor was computed using the matrix E. For the conjugate gradient reconstruction, the g-factor was computed by reconstructing a time-series of simulated noise-only images (Gaussian distribution, mean: 0, standard deviation: 1) with and without SURE acceleration. The corresponding g-factor is given by the ratioσSURE(r)R×σFULL(r), where σSURE(r) and σFULL(r) are the standard deviations along the time dimension for each of the noise-only reconstructions and R is the overall sampling reduction factor (Eggers et al., 2005).

Signal to noise ratio (SNR)

Using the property that the SNR in MRI is proportional to the square root of the acquisition time and the voxel volume (Macovski, 1996); the SNR of SURE-SENSE reconstruction (SNRSURE) with respect to the SNR of the low resolution DFT reconstruction (SNRlow) is given by:

SNRSURE=SNRlowKg,
(8)

where K is the effective spatial resolution gain and g is the noise amplification factor as defined above. Note that the SNR is spatially varying since all the quantities involved are spatially varying. Using the same property, the relationship between SNRSURE and the SNR of the fully-sampled DFT reconstruction (SNRfull) is given by:

SNRSURE=RKgSNRfull.
(9)

Note that if the effective gain in spatial resolution approaches the theoretical limit (K=R), Eq. (9) is similar to the SNR relationship in standard SENSE.

Simulation experiments

One-dimensional simulation

Simulation with one spatial dimension was performed assuming an 8-channel planar array with coil sensitivities computed according to the Biot–Savart law. The array was simulated using a field of view of 256×256 mm2 and rectangular coils of size 40×256 mm2 located along the first spatial dimension with a 10% overlap between adjacent coils. The coil sensitivities were computed for a region of interest located at 80 mm from the array using two spatial grids: 64×64 and 32×32. The central vertical line was used for the one-dimensional simulation. Multi-coil data were generated multiplying an object given by the central vertical line of the Shepp–Logan phantom (Fig. 2) by the simulated coil sensitivities. The central vertical line is represented by column 32 of the 64×64 phantom and by column 16 of the 32×32 phantom. A superresolution factor of 2 was tested on the simulated low resolution data given by the central k-space region, i.e. discarding the outer k-space values of the full k-space data. For each reconstruction, the corresponding regularization parameter was set to have an average g-factor of 2.5.

Fig. 2
1D simulation with two-fold superresolution factor for two different target spatial resolutions: (a) 64-point grid and (b) 32-point grid. W is the FOV length (256 mm). Accelerated acquisition was simulated by discarding the outer k-space values of the ...

Two-dimensional simulation

Simulation with two spatial dimensions were performed using the Biot–Savart model of the 32-element head array coil with soccer-ball geometry (Wiggins et al., 2006) which is used in the in vivo experiments and the Shepp–Logan phantom as object function. Coil sensitivity maps were simulated using a field of view of 220×220 mm2 and an image matrix of 128×128. Noise-free multi-coil data were generated by multiplying the numerical phantom with the sensitivity maps. Gaussian noise corresponding to SNR=100 was added to simulate the SNR that might be measured in an fMRI experiment. 2D superresolution factors of 2×2 and 4×4 were tested on the low resolution data given by the central 64×64 and 32×32 k-space region of the full k-space data respectively.

In vivo experiments

Human brain data were acquired using a 3 Tesla MR scanner (Tim Trio, Siemens Medical Solutions, Erlangen, Germany) equipped with Sonata gradients (maximum amplitude: 40 mT/m, slew rate: 200 mT/m/ms). A 32-element head array coil with soccer-ball geometry which provides sensitivity encoding along all the spatial dimensions was used for RF reception (Wiggins et al., 2006), while RF transmission was performed with a quadrature body coil.

Coil sensitivity calibration was performed using unprocessed in vivo sensitivity references (Sodickson and McKenzie, 2001). In this approach, multi-coil reference images are employed directly as coil sensitivities to solve the inverse problem, followed by post multiplication by the sum-of-squares combination of the reference images to remove additional magnetization density information introduced by the use of the unprocessed reference images. In other words, the image reconstructed by the inversion of an encoding matrix constructed from raw coil reference data is the pixelwise quotient of the true image and the reference combination; therefore the true image can be recovered by post-multiplying the result by the reference combination. This approach is preferred, since the spatial smoothing inherent in explicit coil sensitivity estimation methods such as polynomial fitting (Pruessmann et al., 1999) may limit the performance of SURE-SENSE.

Functional MRI

A visual stimulation experiment with 8 blocks of 16 s of visual fixation and 16 s of flashing checkerboard was performed. Single slice data were acquired using an interleaved echo-planar imaging (EPI) sequence (repetition time (TR): 4 s, echo time (TE): 30 ms, spatial matrix: 256×256, field of view (FOV): 220×220 mm2, slice thickness: 3.4 mm, 64 scan repetitions). The high resolution reference was obtained from the first scan using the full k-space matrix and SURE-SENSE reconstruction was applied to the following down-sampled scan repetitions. The down-sampled data are given by the central 32×32 k-space data discarding the outer k-space data. In order to have a target spatial resolution with sufficient intra-voxel coil sensitivity variation, SURE-SENSE reconstruction was performed using the 32×32 central k-space matrix with a 64×64 target k-space matrix, which represents a two-dimensional sampling reduction factor of R=2×2. A 64×64 k-space matrix is usually employed in fMRI with high temporal resolution (Lin et al., 2006). For comparison, the original fully-sampled 64×64 data and the down-sampled 32×32 data were conventionally reconstructed by applying a discrete Fourier transform (DFT) to each channel and the composite image was computed using a sensitivity-weighted combination. Additionally, a standard SENSE reconstruction was applied to a regularly under-sampled data set with reduction factor of R=4×1, i.e. the fully-sampled 64×64 k-space data matrix was decimated by keeping the first row of every consecutive four rows. Note that standard SENSE does not allow for acceleration of the readout dimension therefore a higher one-dimensional acceleration was used to match the SURE-SENSE acceleration. Correlation and region of interest (ROI) analyses were performed using the TurboFire software package (Posse et al., 2001) with a correlation coefficient threshold of 0.6, i.e. both positive correlation values above 0.6 and negative correlation values below −0.6 were included in the activation map. The reference vector defined by the stimulation paradigm was convolved with the canonical hemodynamic response function defined in SPM99 (Friston et al., 1995). Motion correction and spatial filters were not employed. Average SNR was computed as the ratio of the mean value and standard deviation of the reconstructed time-series data along the temporal dimension.

Spectroscopic imaging

Human brain MRSI data with two spatial dimensions were acquired with Proton Echo Planar Spectroscopic Imaging (PEPSI) (Posse et al., 1995) in axial orientation using a 64×64×512 spatial-spectral matrix (x,y,v) where x and y are the spatial dimensions and v is the spectral dimension. The FOV was 256×256 mm2 and the slice thickness was 20 mm resulting in a nominal voxel size of 320 mm3 (in-plane nominal pixel size was 4×4 mm2). Data were filtered in k-space using a Hamming window which increased the voxel size to 820 mm3 (in-plane effective pixel size was 6.4×6.4 mm2). The spectral width was set to 1087 Hz. The 2D-PEPSI sequence consisted of water-suppression (WS), outer-volume-suppression (OVS), spin-echo RF excitation, phase-encoding for y and the echo-planar readout module for simultaneous encoding of x and t. Data acquisition included water-suppressed (WS) and non-water-suppressed (NWS) scans. The NWS scan was performed without the WS and OVS modules and it is used for spectral phase correction, eddy current correction and absolute metabolite concentration. The high resolution NWS and WS PEPSI data sets were acquired in 2 min each using single signal average, TR=2 s and TE=15 ms. Data were collected with 2-fold oversampling for each readout gradient separately to improve regridding performance and using a ramp sampling delay of 8 μs to limit chemical shift artifacts. Regridding was applied to correct for ramp sampling distortion of the kx-t trajectory. Spectral water images from the high resolution NWS data were employed for coil sensitivity calibration as described before in our SENSE-PEPSI implementations (Lin et al. 2007, Otazo et al. 2007). SURE-SENSE was applied to the central 32×32 k-space matrix using the 64×64 coil sensitivities, which represents a two-dimensional sampling reduction factor of R=2×2. Water images were obtained by spectral integration of the reconstructed NWS data. Lipid images were computed by spectral integration of the reconstructed WS data from 0 to 2.0 ppm. Metabolite images were obtained by spectral fitting using LCModel (Provencher, 1993) with analytically modeled basis sets. Spectral fitting errors in LCModel were computed using the Cramer–Rao Lower Bound (CRLB, the lowest bound of the standard deviation of the estimated metabolite concentration expressed as percentage of this concentration), which when multiplied by 2.0 represent 95% confidence intervals of the estimated concentration values. A threshold of 30% was imposed on the CRLB to accept voxels in the metabolite concentration maps. Average SNR in the WS data was computed using the SNR value from LCModel which is given by the ratio of the maximum in the spectrum-minus-baseline over the analysis window to twice the standard deviation of the residuals. Error maps were computed as the difference with respect to the concentration map from the fully-sampled DFT reconstruction.

Results

One-dimensional simulation

SURE-SENSE reconstruction increased significantly the spatial resolution of conventional Fourier reconstruction of the low resolution data (Fig. 2). The increase in spatial resolution is more significant at lower target spatial resolution (b) than at higher target spatial resolution (a) as depicted by the PSF reconstruction due to the higher intra-voxel variation of the coil sensitivity maps. The nominal mean gain in spatial resolution given by the SURE-SENSE reconstruction without regularization was 1.89 for 64-point target resolution (a) and 1.98 for 32-point target resolution (b). For both cases, the noise-free SURE-SENSE reconstruction is similar to the full DFT reconstruction. However, the noise amplification in the reconstruction (g-factor) presented high mean values of 1.2×104 and 5.2×102 respectively. Using a regularization parameter adjusted to have a mean g-factor of 2.5 for both cases, the mean gain in spatial resolution was reduced to 1.64 and 1.88 respectively. Note that the effective gain is closer to the nominal gain for the smaller 32-point target (b). Constraining the target spatial resolution, which is equivalent to having a smaller extrapolation region in k-space provides a better conditioning of the inverse problem which decreases the SNR penalty, and thus the effective superresolution gain converges to the nominal gain.

Two-dimensional simulation

The noise-free SURE-SENSE reconstruction with R=2×2 was similar to the target object presenting an effective gain close to the theoretical 2×2 increase in spatial resolution (Fig. 3a). Two-dimensional acceleration presents lower and more uniform noise amplification than using a high one-dimensional acceleration (Weiger et al., 2002). For R=4×4, the noise-free SURE-SENSE reconstruction is less similar to the target object presenting slight blurring at the center of the image and residual Gibbs ringing due to the more stringent requirements on the intra-voxel coil sensitivity variations to allow for 16-fold increase in spatial resolution. Nevertheless, the SURE-SENSE reconstruction presents a significant increase in spatial resolution when comparing to the DFT-ZF reconstruction. However, this ideal increase in spatial resolution imposed an SNR penalty as it is shown in the SURE-SENSE reconstruction without regularization in the case of noisy data. SURE-SENSE with Tikhonov regularization controlled the noise amplification in the inverse reconstruction at the expense of reducing the gain in spatial resolution. Note that the SNR penalty is higher for R=4×4 as expected from the larger k-space extrapolation region to recover the full k-space data. Even though the target spatial resolution was not feasible due to the SNR penalty, SURE-SENSE reconstruction with regularization presented a significant gain in spatial resolution and reduction of Gibbs ringing when compared to the conventional DFT with zero-filling reconstruction of the low resolution data.

Fig. 3
Shepp–Logan phantom simulation for superresolution factors of (a) 2×2 and (b) 4×4 with target image matrix of 128×128 using a Biot–Savart model of the 32-channel soccer-ball array. Left: noise-free reconstruction ...

Functional MRI

SURE-SENSE reconstruction of the low spatial resolution time-series fMRI data yielded a significant increase in spatial resolution when compared to the DFT reconstruction with zero-filling (Fig. 4). The gain in spatial resolution is spatially varying, with higher values at positions where the coil sensitivity variation is stronger, e.g. in the periphery of the brain for the soccer-ball array coil. The average gain in spatial resolution was a factor of 2.51 with a maximum value of 3.87 (very close to the theoretical gain of 4) in peripheral regions and a minimum value of 1.55 in central regions where the coil sensitivities are broad and overlapped. This behavior can be explained considering that the regularization approach is broadening the PSF at positions with high nominal g-factor in order to obtain an adequate SNR. The resulting g-factor map after regularization presented a homogeneous distribution with a mean value of 1.54±0.44. Note that without regularization the g-factor map will look inhomogeneous with high values at regions with less intra-voxel coil sensitivity variation, e.g. central region for the soccer-ball array. The average SNR computed from the reconstructed time-series data was 14.4 for the fully-sampled DFT reconstruction, 23.9 for the zero-filled DFT reconstruction and 8.6 for SURE-SENSE. The average SNR for SURE-SENSE is higher than the value predicted by Eq. (8) (6.2), and by Eq. (9) (7.5), which is in part due to the relatively small number of temporal points used to estimate the SNR (64) and the effect of the shape of the PSF which is not completely taken into account since the factor K in Eq. (8) and Eq. (9) is just a FWHM ratio. Nevertheless, they are approximately in the range predicted by Eq. (8) and Eq. (9).

Fig. 4
(a) Reconstruction of fMRI data: fully-sampled data (64×64 matrix) using DFT reconstruction (DFT-FULL), low spatial resolution data (32×32 matrix, R=2×2) using DFT with zero-filling (DFT-ZF) and SURE-SENSE, and regularly undersampled ...

Standard SENSE reconstruction of data regularly undersampled by a factor of R=4×1 to match the SURE-SENSE reduction factor resulted in residual aliasing artifacts and localized noise enhancement whereas the spatial resolution was homogeneous and similar to the fully-sampled DFT reconstruction (Fig. 4a). Note that standard SENSE only removes the aliasing at the center of the voxel and residual aliasing artifacts due to intra-voxel coil sensitivity variations are present in the reconstructed image (Zhao et al., 2005).

The Tikhonov regularization along with pre-conditioning presented a fast and stable reconstruction for SURE-SENSE using the conjugate gradient algorithm with 12 iterations. Note that without the Tikhonov term the g-factor continuously increases with the number of iterations and the stopping point should be selected before convergence to maintain adequate SNR. The total reconstruction time was approximately 2 min (2 s per temporal repetition) using Matlab (The MathWorks, Natick, MA) on a 64-bit quad core workstation.

Visual activation maps obtained from the time-series data reconstructed with SURE-SENSE displayed a spatial pattern closer to the maps from the fully-sampled data when comparing to the zero-filled DFT reconstruction (Fig. 5a). The activation pattern with zero-filled DFT reconstruction is smooth due to partial volume effects and areas with small signal decrease (blue) become visible due to increased SNR in the low spatial resolution data. The activation pattern of SURE-SENSE is spatially more localized than in the zero-filled DFT reconstruction as expected from the increase in spatial resolution and the areas with negative signal changes are less visible due to increased noise level, similar to the high resolution data. SURE-SENSE presented an activation map with spatially varying spatial resolution following the pattern given by the K-map (Fig. 4b). Fig. 5b shows the activation maps for conventional DFT with zero-filling reconstruction of k-space data acquired with different spatial resolutions. The activation map for SURE-SENSE is very close to the DFT-64×64 (full k-space data) at the edges of the brain where the pattern is very discrete. As we move towards the center, the SURE-SENSE map falls in between the DFT-40×40 and DFT-48×48 which is consistent with the 2.5-fold increase in spatial resolution predicted by the K-map for that region. The activation map of standard SENSE reconstruction also suffered from some blurring and additionally presented smaller spatial extent of activation than the DFT-64×64 reconstruction (Fig. 5a) due to increased noise level. This decrease in spatial extent of activation is also more localized at the upper right part of the activation region, which is due in part to the residual aliasing in that region.

Fig. 5
(a) Activation maps for fully-sampled data using DFT reconstruction (DFT-FULL), for low spatial resolution data (R=2×2) using DFT with zero-filling (DFT-ZF) and SURE-SENSE reconstruction and for regularly undersampled data (R=4×1) using ...

Spectroscopic imaging

SURE-SENSE reconstruction reduced the strong effect of k-space truncation in the low spatial resolution PEPSI data set, resulting in metabolite maps with improved spatial resolution and spectra with reduced lipid contamination as compared to the DFT with zero-filling reconstruction. The spatial resolution gain and g-f actor maps displayed a similar pattern as in the fMRI experiment, since the same coil array, acceleration factor and image matrix were employed. The average gain in spatial resolution was 2.37 with a maximum value of 3.72 in peripheral regions and a minimum value of 1.48 in central regions. The average g-factor was 1.49±0.36. The average SNR of the WS reconstructed data was 11.4 for the fully-sampled DFT reconstruction, 19.8 for the zero-filled DFT reconstruction and 7.9 for SURE-SENSE. These values are approximately in the range predicted by Eq. (8) (5.7) and Eq. (9) (6.5).

The maps of NAA and creatine for SURE-SENSE were similar to ones from the DFT reconstruction of the full k-space data. Average error with respect to map derived from the full k-space data were 8.9% for NAA and 8.2% for creatine (Fig. 6). The average errors of the corresponding zero-filled DFT reconstruction NAA and creatine maps were 24.1% and 22.9%. Table 1 shows the average absolute concentration and Cramer–Rao lower bound (CRLB) values from LCModel fitting. The accuracy of spectral quantification as indicated by the CRLB decreased for SURE-SENSE as compared to the fully-sampled DFT reconstruction due to the lower SNR and the presence of residual lipid contamination. However, the strong lipid contamination in the zero-filled DFT reconstruction was highly reduced by SURE-SENSE (Fig. 7) resulting in CRLB values decreased by a factor of 1.46 for NAA and 1.38 for creatine on average.

Fig. 6
Water and lipids images using spectral integration, and NAA and creatine concentration maps using spectral LCModel fitting for fully-sampled data (64×64 matrix) using DFT reconstruction (DFT-FULL) and for low spatial resolution data (32×32 ...
Fig. 7
Absorption mode spectra (DATA) and corresponding LCModel fit (FIT) for a central gray matter voxel and a white matter voxel from the DFT reconstruction of the full k-space data (DFT-FULL), DFT with zero-filling (DFT-ZF) and SURE-SENSE reconstruction of ...
Table 1
Mean value and standard deviation of the absolute concentration and Cramer–Rao lower bound for DFT reconstruction of the fully-sampled data (DFT-FULL) and for DFT with zero-filling (DFT-ZF) and SURE-SENSE reconstruction of the low spatial resolution ...

The total reconstruction time was approximately 34 min (2 s per spectral point of the positive and negative echo data respectively) using Matlab (The MathWorks, Natick, MA) on a 64-bit quad core workstation.

Discussion

The idea of superresolution parallel imaging represents a transition from standard parallel imaging techniques such as SENSE (Pruessmann et al.,1999) and GRAPPA (Griswold et al., 2002), which employ k-space undersampling, to techniques with minimal gradient encoding such as inverse MRI (Lin et al., 2006) which are based on the acquisition of a single k-space point at the center. Acquiring a central k-space region instead of a single point reduces the severe ill-conditioning of the encoding equation in the inverse MRI method and reconstruction with higher spatial resolution is feasible with SURE-SENSE at the expense of decreasing the high acceleration factor. On the other hand, while relatively stronger intra-voxel coil sensitivity variation leads to residual aliasing artifacts in standard parallel imaging, it improves the performance of SURE-SENSE. Moreover, the reconstruction error for standard SENSE reconstruction is distributed over the entire object while for SURE-SENSE reconstruction errors are limited to the extent of the low resolution voxel. This reflects the difference in k-space sampling: localized k-space errors in standard SENSE lead to distributed artifacts in the image space, whereas errors in reconstructing extrapolated k-space in SURE-SENSE lead to localized errors in the image space. Application of SURE-SENSE to spectroscopic imaging is therefore particularly advantageous since one of the major challenges of standard SENSE is the unfolding of positions which include lipid resonances from peripheral regions. Imaging cortical regions with SURE-SENSE is also advantageous due to the increased intra-voxel coil sensitivity variation for positions close to the coil elements.

SURE-SENSE presents a different reconstruction paradigm when compared to the minimum-norm SENSE method (Sanchez-Gonzalez et al., 2006), even though both methods employ high resolution coil sensitivity maps. SURE-SENSE and minimum-norm SENSE are both implementations of strong SENSE (Pruessmann et al., 1999) and the equivalent voxel shape fitting formalism (Sodickson and McKenzie, 2001) using different target voxel shapes. For SURE-SENSE the target voxel shape is given by delta functions in the high resolution grid whereas for minimum-norm SENSE the target voxel shape is given by delta functions in the low resolution grid. The minimum-norm SENSE method aims to reconstruct a sharp PSF (at the intra-voxel level) but without increasing the underlying image spatial resolution, i.e. one reconstructed voxel is created for each original low resolution voxel. In SURE-SENSE, the goal is to reconstruct multiple voxels for each original low resolution voxel.

The maximum superresolution factor for SURE-SENSE depends on a number of factors, including the spatial nonuniformity of the coil sensitivity patterns and efficiency of regularization in the reconstruction process. SURE-SENSE is most suitable for reconstruction of intrinsic low spatial resolution MRI modalities where the acquisition of a reference image with high spatial resolution does not impose a time penalty, e.g. functional MRI and spectroscopic imaging. The spatial resolution in SURE-SENSE is spatially varying as is the case for inverse MRI. Regions with broad and overlapping coil sensitivity distributions are reconstructed with lower spatial resolution in order to maintain an adequate SNR. In practice, for interpretation of SURE-SENSE functional images, the K-map may be used to judge how much confidence to place in the focality of activation in any particular region. The proposed method is expected to be able to reconstruct images with a spatial resolution up to 3×3×3 mm3 and with 4- to 6-fold two-dimensional superresolution factors using an array coil similar to the 32-channel soccer-ball array. This target resolution is appropriate for functional MRI with high temporal resolution which is usually performed with spatial resolution of about 3 mm (Lin et al., 2006. For example, using a FOV around 200×200 mm2 and coil sensitivities acquired with a 64×64 matrix, it should be feasible to perform SURE-SENSE reconstruction using the 28×28 or 32×32 acquisition matrix. Higher superresolution acceleration factors should be feasible for lower target resolution, e.g. 11×11 raw data matrix to achieve a 32×32 target matrix (R=3×3) or 4×4 raw data matrix to achieve a 16×16 target matrix (R=4×4). We expect that the performance of SURE-SENSE reconstruction for human brain imaging will improve with the number of receiver coils using array geometries similar to the soccer-ball geometry of the 32-element array used in this work. Moreover, implementation at high magnetic field strength (7 T) will provide stronger spatial modulation of the coil sensitivity maps (Ohliger et al., 2003, Wiesinger et al., 2004). We are in the process of implementing the technique using a 96-channel array at 3 T (Wiggins et al., 2007) and a 32-channel array at 7 T with a 128×128 target image matrix.

SURE-SENSE is a parallel imaging technique that can be applied across a wide range of pulse sequence types like standard parallel imaging techniques. The principal advantage of SURE-SENSE in traditional applications such as functional MRI and MR spectroscopic imaging is not to get very high spatial resolution per se, but rather to accelerate the acquisition of an image with a target resolution that presents intra-voxel coil sensitivity variation. SURE-SENSE allows for acceleration of the readout dimension in echo-planar trajectories and therefore the overall SURE acceleration can be divided along the two spatial dimensions with the advantage of having lower and more uniform noise amplification than is possible with a high one-dimensional acceleration (Weiger et al., 2002). Readout acceleration in SURE-SENSE would require more extensive ramp sampling which only introduces a small SNR loss (Otazo et al., 2006) and the only extra step required is the regridding of the non-equidistant data in k-space to a Cartesian grid which does not introduce any significant extra blurring, in particular for SURE-SENSE acquisition that is performed at low spatial resolution. SURE-SENSE would be applicable to spiral trajectories as well reducing the extent of the spiral in k-space, thereby reducing gradient slew rate constraints and possible gradient induced peripheral nerve stimulation. An alternative way to accelerate the readout dimension is to increase the gradient strength/bandwidth and pay a well-defined square root of bandwidth price in SNR as opposed to pay only slightly smaller price with SURE-SENSE that results in spatially-varying spatial resolution. However, readout gradients are already operating at the limits imposed by physiological and technical constraints and further increases in gradient strength/bandwidth might not be feasible. Moreover, the capability of accelerating the readout dimension in echo-planar trajectories can be used to increase the spectral bandwidth in PEPSI, especially at high magnetic strength where the gradient performance is currently limited. In general, SURE-SENSE acceleration will reduce the sensitivity to motion and it is particularly interesting for faster 3D encoding where the available SNR is higher. Moreover, the high temporal resolution provided for fMRI can be used to correlate between hemodynamic onset timing and neuronal events (recorded by MEG or EEG).

High resolution coil sensitivity calibration is required for optimal performance of the SURE technique. The conventional way to estimate coil sensitivity maps for SENSE using spatial smoothing may limit the performance of the method. The “in vivo” sensitivity calibration method employed in this work allows for the use of coil sensitivity information with full spatial resolution. In future work we will explore k-space superresolution reconstruction with autocalibration, such as in the GRAPPA method (Griswold et al., 2002), which would represent an alternative to SENSE based reconstruction that includes all the high resolution information available from the coil sensitivity profiles.

The current implementation of SURE-SENSE with two-dimensional acceleration is performed using an iterative conjugate gradient algorithm, which might not be optimal in terms of reconstruction time, but which represents a practical solution to the computationally intensive inverse problem. Faster reconstruction can be accomplished using direct inversion methods similar to the non-iterative approach for non-Cartesian parallel MRI as described in Ying et al. (2005). Experiments that consist of series of images such as fMRI and MRSI represent the real bottleneck in reconstruction time, since several SURE-SENSE reconstructions need to be performed sequentially. Parallel computing can be used to substantially reduce the total reconstruction time since each individual reconstruction is parallelizable.

In conclusion, SURE-SENSE represents a novel alternative to standard SENSE particularly for accelerating imaging applications which use intrinsically low spatial resolution and echo-planar trajectories. Future work will characterize the optimal operation regimes of the method with respect to target spatial resolution, array geometry and number of receiver elements, and will explore implementation at 7 T to improve the tradeoff between spatial and temporal resolution.

Acknowledgments

This study was supported by National Institutes of Health Grants R01 HD040712, R01 NS037462, R01 EB000790-04, P41 RR14075, R01 EB002618-01, R21 EB007298, R01 EB002468, the Mental Illness and Neuroscience Discovery Institute (MIND), the Ibero-American Science and Technology Education Consortium (ISTEC), the Electrical and Computer Engineering Department of the University of New Mexico, NHRI-EX97-9715EC (National Health Research Institute, Taiwan), NSC 97-2320-B-002-058 (National Science Council, Taiwan), and NSC 97-2221-E-002-005 (National Science Council, Taiwan). We thank Lawrence Wald (Massachusetts General Hospital) for helpful discussions and technical support and Pierre-Gilles Henry (University of Minnesota) for the simulation of basis sets.

References

  • Belliveau JW, Kennedy DN, McKinstry RC, Buchbinder BR, Weisskoff RM, Cohen MS, Vevea JM, Brady TJ, Rosen BR. Functional mapping of the human visual cortex by magnetic resonance imaging. Science. 1991;254 (5032):716–719. [PubMed]
  • Brown TB, Kincaid BM, Ugurbil K. NMR chemical shift imaging in three dimensions. Proc Natl Acad Sci U S A. 1982;79 (11):3523–3526. [PubMed]
  • Eggers H, Mazurkewitz P, Boesiger P. Nose amplification in non-Cartesian parallel imaging. Proceedings of the 13th ISMRM Annual Meeting; Miami, USA. 2005. p. 2429.
  • Elad M, Feuer A. Restoration of a single superresolution image from several blurred, noisy, and undersampled measured images. IEEE Trans Image Process. 1997;6 (12):1646–1658. [PubMed]
  • Friston KJ, Holmes AP, Worsley KJ, Poline JB, Frith CD, Frackowiak RSJ. Statistical parametric maps in functional neuroimaging: a general linear approach. Hum Brain Map. 1995;2:189–210.
  • Griswold MA, Jakob PM, Heidemann RM, Nittka M, Jellus V, Wang J, Kiefer B, Haase A. Generalized autocalibrating partially parallel acquisitions (GRAPPA) Magn Reson Med. 2002;47 (6):1202–1210. [PubMed]
  • Hansen PC. Rank-deficient and Discrete Ill-posed Problems: Numerical Aspects of Linear Inversion. 2. SIAM; Philadelphia: 1998.
  • Jones RA, Haraldseth O, Muller TB, Rinck PA, Oksendal AN. K-space substitution: a novel dynamic imaging technique. Magn Reson Med. 1993;29 (6):830–834. [PubMed]
  • Liang ZP, Lauterbur PC. A generalized series approach to MR spectroscopic imaging. IEEE Trans Med Imaging. 1991;10 (2):132–137. [PubMed]
  • Liang ZP, Boada F, Constable T, Haacke EM, Lauterbur PC, Smith MR. Constrained reconstruction methods in MR imaging. Rev Magn Reson Med. 1992;4:67–185.
  • Lin FH, Kwong KK, Belliveau JW, Wald LL. Parallel imaging reconstruction using automatic regularization. Magn Reson Med. 2004;51 (3):559–567. [PubMed]
  • Lin FH, Huang TY, Chen NK, Wang FN, Stufflebeam SM, Belliveau JW, Wald LL, Kwong KK. Functional MRI using regularized parallel imaging acquisition. Magn Reson Med. 2005;54 (2):343–353. [PMC free article] [PubMed]
  • Lin FH, Wald LL, Ahlfors SP, Hamalainen MS, Kwong KK, Belliveau JW. Dynamic magnetic resonance inverse imaging of human brain function. Magn Reson Med. 2006;56 (4):787–802. [PubMed]
  • Lin FH, Tsai SY, Otazo R, Caprihan A, Wald LL, Belliveau JW, Posse S. Sensitivity-encoded (SENSE) proton echo-planar spectroscopic imaging (PEPSI) in human brain. Magn Reson Med. 2007;57 (2):249–257. [PubMed]
  • Lin FH, Witzel T, Mandeville JB, Polimeni JR, Zeffiro TA, Greve DN, Wiggins G, Wald LL, Belliveau JW. Event-related single-shot volumetric functional magnetic resonance inverse imaging of visual processing. Neuroimage. 2008;42 (1):230–247. [PMC free article] [PubMed]
  • Macovski A. Noise in MRI. Magn Reson Med. 1996;36 (3):494–497. [PubMed]
  • McDougall MP, Wright SM. 64-channel array coil for single echo acquisition magnetic resonance imaging. Magn Reson Med. 2005;54 (2):386–392. [PubMed]
  • Ohliger MA, Grant AK, Sodickson DK. Ultimate intrinsic signal-to-noise ratio for parallel MRI: electromagnetic field considerations. Magn Reson Med. 2003;50 (5):1018–1030. [PubMed]
  • Otazo R, Mueller B, Ugurbil K, Wald LL, Posse S. Signal-to-noise ratio and spectral linewidth improvements between 1.5 and 7 Tesla in proton echo-planar spectroscopic imaging. Magn Reson Med. 2006;56 (6):1200–1210. [PubMed]
  • Otazo R, Tsai SY, Lin FH, Posse S. Accelerated short-TE 3D proton echo-planar spectroscopic imaging using 2D-SENSE with a 32-channel array coil. Magn Reson Med. 2007;58 (6):1107–1116. [PubMed]
  • Peled S, Yeshurun Y. Superresolution in MRI: application to human white matter fiber tract visualization by diffusion tensor imaging. Magn Reson Med. 2001;45 (1):29–35. [PubMed]
  • Plevritis S, Macovski A. Spectral extrapolation of spatially bounded images. IEEE Trans Med Imaging. 1995;14 (3):487–497. [PubMed]
  • Posse S, Tedeschi G, Risinger R, Ogg R, Le Bihan D. High speed 1H spectroscopic imaging in human brain by echo planar spatial–spectral encoding. Magn Reson Med. 1995;33 (1):34–40. [PubMed]
  • Posse S, Binkofski F, Schneider F, Gembris D, Frings W, Habel U, Salloum JB, Mathiak K, Wiese S, Kiselev V, Graf T, Elghahwagi B, Grosse-Ruyken ML, Eickermann T. A new approach to measure single-event related brain activity using real-time fMRI: feasibility of sensory, motor, and higher cognitive tasks. Hum Brain Mapp. 2001;12 (1):25–41. [PubMed]
  • Provencher SW. Estimation of metabolite concentrations from localized in vivo proton NMR spectra. Magn Reson Med. 1993;30 (6):672–679. [PubMed]
  • Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: sensitivity encoding for fast MRI. Magn Reson Med. 1999;42 (5):952–962. [PubMed]
  • Pruessmann KP, Weiger M, Bornert P, Boesiger P. Advances in sensitivity encoding with arbitrary k-space trajectories. Magn Reson Med. 2001;46 (4):638–651. [PubMed]
  • Sodickson DK, Manning WJ. Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays. Magn Reson Med. 1997;38 (4):591–603. [PubMed]
  • Sodickson DK, McKenzie CA. A generalized approach to parallel magnetic resonance imaging. Med Phys. 2001;28 (8):1629–1643. [PubMed]
  • Sanchez-Gonzalez J, Tsao J, Dydak U, Desco M, Boesiger P, Pruessmann KP. Minimum-norm reconstruction for sensitivity-encoded magnetic resonance spectroscopic imaging. Magn Reson Med. 2006;55 (2):287–295. [PubMed]
  • Weiger M, Pruessmann KP, Boesiger P. 2D SENSE for faster 3D MRI. MAGMA. 2002;14 (1):10–19. [PubMed]
  • Wiesinger F, Van de Moortele PF, Adriany G, De Zanche N, Ugurbil K, Pruessmann KP. Parallel imaging performance as a function of field strength — an experimental investigation using electrodynamic scaling. Magn Reson Med. 2004;52 (5):953–964. [PubMed]
  • Wiggins GC, Triantafyllou C, Potthast A, Reykowski A, Nittka M, Wald LL. 32-channel 3 Tesla receive-only phased-array head coil with soccer-ball element geometry. Magn Reson Med. 2006;56 (1):216–223. [PubMed]
  • Wiggins GC, Alagappan V, Potthast A, Schmitt M, Wiggins CJ, Fischer H, Jahns K, Benner T, Polimeni J, Wald LL. Design optimization and SNR performance of 3 T 96 channel phased array head coils. Proceedings of the 15th ISMRM Annual Meeting; Berlin, Germany. 2007. p. 243.
  • Ying L, Haldar J, Liang ZP. An efficient non-iterative reconstruction algorithm for parallel MRI with arbitrary k-space trajectories. Proceedings of the 27th Conf IEEE Eng Med Biol Soc; 2005. pp. 1344–1347. [PubMed]
  • Zhao X, Prost RW, Li Z, Li SJ. Reduction of artifacts by optimization of the sensitivity map in sensitivity-encoded spectroscopic imaging. Magn Reson Med. 2005;53 (1):30–34. [PubMed]