Magnetic Resonance in Medicine 58:1257–1265 (2007) Non-Cartesian Data Reconstruction Using GRAPPA Operator Gridding (GROG) Nicole Seiberlich,1 Felix A. Breuer,2 Martin Blaimer,3 Kestutis Barkauskas,4 Peter M. Jakob,1,2 and Mark A. Griswold3 A novel approach that uses the concepts of parallel imaging to grid data sampled along a non-Cartesian trajectory using GRAPPA operator gridding (GROG) is described. GROG shifts any acquired data point to its nearest Cartesian location, thereby converting non-Cartesian to Cartesian data. Unlike other parallel imaging methods, GROG synthesizes the net weight for a shift in any direction from a single basis set of weights along the logical k-space directions. Given the vastly reduced size of the basis set, GROG calibration and reconstruction requires fewer operations and less calibration data than other parallel imaging methods for gridding. Instead of calculating and applying a density compensation function (DCF), GROG requires only local averaging, as the reconstructed points fall upon the Cartesian grid. Simulations are performed to demonstrate that the root mean square error (RMSE) values of images gridded with GROG are similar to those for images gridded using the gold-standard convolution gridding. Finally, GROG is compared to the convolution gridding technique using data sampled along radial, spiral, rosette, and BLADE (a.k.a. periodically rotated overlapping parallel lines with enhanced reconstruction [PROPELLER]) trajectories. Magn Reson Med 58:1257–1265, 2007. © 2007 Wiley-Liss, Inc. Key words: parallel imaging; GRAPPA; gridding; non-Cartesian; image reconstruction Non-Cartesian imaging has advantages over standard Cartesian imaging (1–7) due to, for example, efficient k-space coverage or suppression of off-resonance effects. However, all acquired points do not necessarily fall onto a grid and must be resampled onto a Cartesian matrix before a Fourier transform can be performed. The convolution gridding approach of Jackson et al. (8) is the gold-standard method to prepare the raw data prior to the fast Fourier transform (FFT), but it requires a density compensation function (DCF) to even out the sampling density throughout the trajectory. Unlike the radial trajectory, computing an analytic DCF is not trivial for trajectories with exotic geome- 1Department of Experimental Physics 5, University of Würzburg, Würzburg, Germany. 2Research Center Magnetic Resonance Bavaria (MRB), Würzburg, Germany. 3Department of Radiology, University Hospitals of Cleveland and Case Western Reserve University, Cleveland, Ohio, USA. 4Department of Biomedical Engineering, Case Western Reserve University, Cleveland, Ohio, USA. Grant sponsor: German Research Society (Deutsche Forschungsgemeinschaft, DFG); Grant number: JA 827/4⫺4; Grant sponsor: Siemens Medical Solutions. *Address reprint requests to: Nicole Seiberlich, Physikalisches Institut, EP 5, Universität Würzburg, Am Hubland, D-97074 Würzburg, Germany. E-mail: neseiber@physik.uni-wuerzburg.de Received 11 February 2007; revised 2 September 2007; accepted 9 September 2007. DOI 10.1002/mrm.21435 Published online 29 October 2007 in Wiley InterScience (www.interscience. wiley.com). © 2007 Wiley-Liss, Inc. tries (9), such as BLADE (a.k.a. periodically rotated overlapping parallel lines with enhanced reconstruction [PROPELLER] (5)), rosette (6), or stochastic trajectories (7). Other gridding methods (10 –13) can be used to grid data acquired along such trajectories without a DCF, although these methods have other disadvantages such as noise enhancement, strenuous parameterization, or iterative reconstructions. As an alternative to convolution gridding, parallel imaging concepts can be applied (14 –22). The basic idea behind parallel imaging is to use several coils at different spatial locations to acquire Nyquist undersampled data, and then to reconstruct the missing points by combining the coil information using a reconstruction algorithm such as generalized autocalibrating partially parallel acquisitions (GRAPPA) (20). In the case of using parallel imaging for gridding, the missing points are the Cartesian points, and the acquired data are the non-Cartesian points. By forming linear combinations of the multichannel non-Cartesian data points, the appropriate spatial harmonics of the missing Cartesian points can be synthesized. A well-suited parallel imaging method to shift non-Cartesian points to nearby Cartesian locations is the GRAPPA operator gridding (GROG) approach (23). As recently described by Griswold et al. (24), the GRAPPA operator formalism has several properties which make it useful for gridding. First, a GRAPPA operator can be used to shift a point in k-space by an arbitrary amount n⌬k: s共kx,ky ⫹ n⌬ky兲 ⬇ Gn 䡠 s共kx,ky兲, [1] where s(kx,ky) is the acquired point, Gn are the appropriate coil weighting factors (weights) for the desired shift, and s(kx,ky⫹n⌬ky) is the vector containing the signal from each receiver coil at the desired location. It is important to note that this formulation is equivalent to GRAPPA with a single source point and a single target point. The weight set is simply a square matrix of size NC ⫻ NC, where NC is the number of coil elements used for the acquisition. Because this weight set is similar to a ladder or propagator operator in quantum mechanics, the term “GRAPPA operator” is used to describe it. The GRAPPA operator can be derived in the same fashion as standard GRAPPA weights; namely, a fit of points with the appropriate relationship is performed with the use of the pseudoinverse: S ACS共kx,ky ⫹ n⌬ky兲 䡠 pinv共SACS共kx,ky兲兲 ⬇ Gn. [2] In these equations, the signal matrices SACS are made up of a collection of signal vectors from the Cartesian autocalibration signal with the appropriate distance relation- 1257 1258 Seiberlich et al. ship, analogous to the weight set determination in GRAPPA. A second important property of the GRAPPA operator is that a small shift operator G␦ can be derived from an operator for a larger shift Gn␦ by taking the nth root of the larger operator (25): G ␦ ⫽ Gn䡠␦1/n, [3] where n is not restricted to integer values. This implies that the calibration step for an arbitrary shift in one direction must be performed only once, as any smaller shift can be derived from a larger shift. The operation above can be performed by first diagonalizing the GROG weight set for the base shift G1, and using this diagonal form to calculate the smaller shift operator G␦ by employing the ␦th power of the eigenvalues: G 1 ⫽ E 䡠 V 䡠 E⫺1 and G␦ ⫽ E 䡠 V␦ 䡠 E⫺1 , [4] where E is a matrix containing the eigenvectors of the base shift matrix, and V is the diagonal matrix containing the eigenvalues. Because the weights are calculated using a Nyquist sampled calibration dataset, the ⌬k ⫽ 1 GROG weights are completely determined, and the principle roots of the appropriate weight set (i.e., those where the eigenvalue has a nonnegative real part) can always be used to calculate the weight sets for smaller shifts. Thus, no time-consuming root determination must be performed as would be the case if employing undersampled calibration data. Finally, GRAPPA operators for shifts in orthogonal directions can be derived and subsequently applied to generate a multi-dimensional shift: s共kx ⫹ ␦x,ky ⫹ ␦y,kz ⫹ ␦z兲 ⫽ G␦共xyz兲 䡠 s共kx,ky,kz兲 ⫽ G␦x 䡠 G␦y 䡠 G␦z 䡠 s共kx,ky,kz兲. [5] Only unit shifts along the logical k-space axes (read, phase, and partition for 3D imaging) need to be calibrated rather than every possible shift, and a shift in any direction by any amount can be calculated from this single group of orthogonal basis weights. GROG makes use of these properties of the GRAPPA operator to rapidly map non-Cartesian data points to the nearest Cartesian grid location. This is performed through the following steps: ● Gx and Gy weights (for a 2D image) are determined using the read lines (Gx) and the phase encoding lines (Gy) of a low-resolution Cartesian calibration signal by fitting each point to the point adjacent to it in the appropriate direction (Eq. [2]). In the case of Nyquistsampled Cartesian data, the base weights Gx and Gy are the GRAPPA operators for a shift of ⌬k ⫽ 1 in either the read or phase encoding direction. ● The distance needed to shift a non-Cartesian point (analogous to the GRAPPA source point) to the nearest Cartesian point (analogous to the GRAPPA target point) is calculated using the k-space trajectory. For FIG. 1. GROG gridding of non-Cartesian points. Cartesian destinations are at the intersections of straight, finely dotted lines, where consecutive samples of an arbitrary trajectory are represented by solid circles. GROG grids a non-Cartesian data point by shifting it to its nearest Cartesian location via an appropriate weight set. For example, the GROG weights Gx0.22 䡠 Gy0.29 are applied to the upper sampling location. When applied throughout the trajectory, GROG results in a purely Cartesian k-space. example, if a data point lies at [kx,ky] ⫽ [19.77,5.83], the nearest Cartesian location is [kx,ky] ⫽ [20,6], and the distance to be shifted is [⌬kx,⌬ky] ⫽ [⫹0.23,⫹0.17]. ● As shown in Eq. [3], the appropriate operators for the small shifts in the x- and y-direction are calculated using Gx and Gy as the base weights. For instance, if a point must be shifted by ⌬kx ⫽ 0.23 and ⌬ky ⫽ 0.17 and the base weights are calibrated for ⌬k ⫽ 1 shifts, the appropriate weight set would be calculated as shown in Eq. [5]: G ␦kx,␦ky ⫽ Gx0.23 䡠 Gy0.17 . [6] ● The shifted points are calculated by applying the appropriate operators to the non-Cartesian data points as shown in Eq. [1] and deposited in the proper Cartesian location in k-space ● In cases where multiple non-Cartesian points map to the same Cartesian grid point, a simple average is performed to arrive at the final value at that Cartesian location A schematic depicting the shifting of non-Cartesian kspace points to their nearest Cartesian locations, and the proper weight calculation, is shown in Fig. 1. The effective DCF for GROG, i.e., the simple averaging of the shifted points that are mapped to the same Cartesian location, can be understood as follows: each non-Cartesian point that is shifted to the same Cartesian point should ideally have the same value (given perfect GROG weights, no noise, and an absence of experimental factors such as relaxation). In other words, after shifting the non-Cartesian GRAPPA Operator Gridding points with GROG, the resulting dataset can be treated as a dataset of Cartesian points, some of which have been acquired multiple times. Thus, the shifted points belonging to the same Cartesian location generated from different non-Cartesian points can simply be averaged. To employ GROG to grid non-Cartesian points, a calibration signal must be used to determine the weights for steps in the Gx and Gy directions (for 2D imaging). For arbitrary k-space trajectories, this calibration signal can be a low-resolution Cartesian dataset where steps of ⌬k ⫽ 1 are performed (i.e., Gx is calculated from the points in the read direction, and Gy from the points along the phase encoding direction). For some trajectories, this low-resolution Cartesian dataset is automatically present, as in the BLADE trajectory, in which one blade can be used to calculate the base weights. It is important to note that for a unit-sampled base grid (i.e., equal and isotropic image and grid FOV), the maximum shift magnitude in one direction is 0.5 ⌬k. While a step of this size is possible to reconstruct with GROG using most commercially available coils, an array with a small number of detector elements or poorly arranged coils can result in artifacts due to insufficient variation of coil sensitivity along an accelerated dimension. Although the trajectory affects the distribution of shifts along a given direction, most shifts are smaller than this maximum. Since the shifts are variable in size and direction, artifacts that result from GROG are expected to be incoherent and appear as a noise enhancement in the image. METHODS Simulation To examine the effects of noise on the gridded images, simulated datasets were constructed to compare convolution gridding and GROG. A standard Shepp-Logan phantom was employed with an eight-element one-ring head coil array, in which sensitivities were derived using an analytic integration of the Biot-Savart equations. The Cartesian data were resampled as radial data (200 projections, 256 readout points, base matrix of 128 ⫻ 128) by sinc interpolation. Progressively larger levels of normal complex noise (real and imaginary standard deviations [␴N] of 0.1, 0.5, 1, and 2.5, corresponding to approximately 0.1%, 0.5%, 0.95%, and 2.36% of the maximum amplitude of the direct current [DC] signal) were applied to the radial kspace data. The data were then gridded with both GROG and the standard convolution gridding in order to compare the root mean square error (RMSE) of the resulting images as compared to the noiseless Shepp-Logan phantom. Note that the region of k-space support for the Cartesian phantom was reduced with a radial mask of the same extent as the radial trajectory to enforce isotropic image resolution; coincidentally, this will slightly affect the appearance of the standard noiseless phantom. Additionally, a 1D interpolation along the Nyquist sampled read direction (i.e., the radial arms) of each of the noisy radial datasets was performed in order to increase the data oversampling in the read direction to a factor of 4. This is expected to improve the GROG reconstruction, as the additional radial points are densely sampled and therefore more non-Cartesian 1259 points contribute to each individual Cartesian location (i.e., each Cartesian point is reconstructed from non-Cartesian points which must be shifted in different directions, which has an averaging effect on the resulting Cartesian points). To examine the effect of different numbers of array elements on the RMSE of the GROG-gridded radial images, several coils were simulated with similar geometries but different numbers of coil elements. The geometry chosen was that of a head coil, i.e., the elements were arranged in a circle around the object. Following the generation of a coil sensitivity map for each of the different coils, radial datasets for each coil were simulated as stated above, and these datasets were gridded using GROG. The RMSE of each image as compared with the noiseless Shepp-Logan phantom was then calculated. In addition, a comparison was made between the inherent DCF resulting from GROG gridding of fully- and undersampled radial datasets and the corresponding DCFs proposed by Pipe (26). To this end, a fully-sampled radial dataset with one Cartesian arm along the x-direction was simulated as described above (128 projections, 82 readout points, base matrix of 82 ⫻ 82, corresponding to a dataset that is sampled according to the Nyquist criterion in the azimuthal direction). The dataset was gridded using GROG, and the contribution of the Cartesian arm to the GROG k-space, i.e., the effective DCF along this radial arm, was examined. In this way, the effect of the averaging of multiple non-Cartesian points that end up at the same Cartesian location could be examined. The effective GROG DCF was then compared to the standard radial Ram-Lak DCF. To examine the effective GROG DCF of an undersampled dataset, three fourths of the radial arms were removed, and the remaining arms (32 projections, R ⫽ 4) were gridded using GROG. The contribution of the radial arm along the x-direction was then compared to the undersampled DCF shown to have the optimal signal-tonoise ratio (SNR) and lowest artifact energy by Pipe (26). Finally, to demonstrate that GROG can be used to grid data acquired along exotic trajectories without the use of a DCF, a rosette dataset was simulated using the eight-channel coil described above. The trajectory is made up of 21,600 points (divided into 180 read points for each of 120 phase encoding steps). Before gridding, the data were interpolated in the read direction to yield a four-fold oversampled dataset, as discussed in the paragraphs above. In Vivo In order to compare GROG with the gold standard of convolution gridding, in vivo radial, spiral, and BLADE datasets were acquired with parameters as given in Table 1. Informed consent from the volunteers was obtained before each study. Since radial and spiral are the most commonly applied trajectories in non-Cartesian imaging, their wellknown artifacts from convolution gridding can be readily compared to the GROG reconstruction results. The BLADE trajectory was selected for two reasons. The first is that at least one of the blades can be used to calculate the base weights, as described in the first section of this article. In addition, GROG completely avoids the computation of an analytic DCF for the BLADE trajectory, which is a nontrivial but necessary process prior to convolution gridding. 1260 Seiberlich et al. Table 1 Scanners and Parameters Used for the In Vivo Images Gridded Using GROG for the Radial, Spiral, and BLADE Trajectories Read points Phase encoding steps Number of coils Base matrix size Scanner Sequence TR (ms) TE (ms) DCF Base grid Kernel width Radial Spiral BLADE 512 256 12 256 1.5T Espree FLASH 20 5 Ram-Lak 2 5 7289a 4 8 192 3T Trio Gradient echo 2500 30 Modified Ram-Lak 2 3 256 410b 17 256 1.5T Espree FLASH 20 5 Iterativec 2 5 aA center-out constant-linear-velocity spiral trajectory was used. 410 phase-encoding steps for the BLADE trajectory are comprised of 41 parallel lines for each of 10 blades. cSee Ref. 28 for a detailed description of the iterative BLADE DCF. FLASH ⫽ fast low-angle shot. bThe Before gridding each dataset with GROG, as described above, the data were interpolated along the read direction to an oversampling factor of 4, which is expected to improve the reconstruction quality. No additional DCF was applied to the data before gridding with GROG. In addition, each dataset was also gridded with standard convolution gridding (parameters shown in Table 1). A Ram-Lak DCF was applied to the radial data before performing convolution gridding, and a modified Ram-Lak to the constant-angular-velocity spiral data (27). An iteratively-calculated BLADE DCF (28) was applied prior to convolution gridding since an analytical DCF has yet to be described. RESULTS The results from the simulations are shown in Figs. 2– 6. As can be seen in the profiles and the RMSE in Fig. 2, images gridded with GROG without oversampling (dashed gray line) contain elevated noise levels when compared to the convolution gridding images (solid gray line). At low noise levels (␴N ⬍ 1), the profiles of both reconstructions are nearly identical to that of the Shepp-Logan phantom (solid black line), implying that these two methods offer similar image quality for high SNR values as the additional noise due to GROG is insignificant. However, at higher FIG. 2. Representative profiles of the results of convolution gridding (solid gray line), GROG without oversampling (dashed gray line), and GROG interpolated to yield an oversampling factor of 4 (thick black line) with different noise levels (␴N), compared to the noiseless Shepp-Logan phantom (thick black line). As the noise level of the dataset increases, the GROG without oversampling shows a considerable increase in noise, as indicated by the RMSE values, whereas the GROG images with oversampling have approximately the same noise level as the convolution gridding images. GRAPPA Operator Gridding FIG. 3. Examples of radial datasets with a noise level of ␴N ⫽ 1 gridded with convolution gridding (left) and GROG (right). As can be seen in the RMSE values shown in Fig. 2, there is no appreciable difference between the convolution gridding and the GROG image. noise levels, the profile in the lower signal areas of the GROG images deviates more strongly from the standard profile than that of the convolution gridding image, which is also reflected in the higher RMSE of the GROG images. SNR losses can be traced to the application of a weight set to a noisy data point; the underlying noise of a given datum can be enhanced by the shifting process, especially for large shifts. However, the RMSE values of the images reconstructed from data interpolated to an oversampling factor of four (thin black line) are similar to those of the convolution gridding images. This effect is to be expected, as the higher oversampling factor in the read direction indicates that more non-Cartesian points are shifted to the same Cartesian point, thus effectively averaging out the errors that occur when the weights are applied to noisy data points. Thus, when gridding low SNR data with GROG, a 1D interpolation along the read direction must be used to increase the read oversampling to avoid SNR 1261 FIG. 5. The effective GROG DCFs (solid lines) and the standard DCFs (dotted lines) for fully-sampled (black) and R ⫽ 4 undersampled (gray) radial datasets. For the central portion of k-space, the GROG DCFs follow the standard radial DCFs closely. In addition, in the R ⫽ 4 case, the GROG DCF follows the optimal SNR DCF for undersampled radial data as proposed by Pipe (26); at the point where the Nyquist criterion is violated in the azimuthal direction, the DCF becomes constant. losses. Alternatively, the oversampling factor in the read direction can be increased, which would yield the same results. Figure 3 shows the ␴N ⫽ 1 images gridded with the standard convolution gridding (left) and with GROG (interpolation factor 4, right). As can be seen in this figure, no additional noise enhancement is observed in the image gridded with GROG. Figure 4 is a plot of the results of the RMSE calculation for radial datasets simulated with different numbers of coil elements and gridded using GROG. As can be seen in the figure, GROG is not effective when two coil elements are FIG. 4. The RMSE values images reconstructed using GROG from radial datasets simulated with different numbers of coil elements. As can be seen in the graph, the use of coils with too few elements leads to a high RMSE, indicating that the image contains artifacts due to the GROG weights. The inset image on the left shows such a case; performing GROG on radial data acquired with two coil elements leads to artifacts in the reconstructed image. As the number of coils increases, the RMSE also decreases, until a minimum is reached. The use of six-coil elements offers approximately the same image quality as coils with more elements (dotted line in chart). The middle and right inset images, reconstructed using six and 12 elements, respectively, are visually indistinguishable from each other and from the reference image. 1262 Seiberlich et al. FIG. 6. An example of an image simulated using the rosette trajectory and gridded with GROG. The trajectory used is shown in the image on the left, where the black path delineates one phase encoding step, and the reference Shepp-Logan phantom is shown in the center. The image resulting from the use of GROG to grid the simulated rosette data is on the right, and is visually indistinguishable from the reference image (RMSE ⫽ 0.17%). used, as only coil sensitivity variation can be employed to determine the GROG weights in one direction. This leads to the appearance of artifacts in the gridded image (inset, left image) and a high RMSE. Similarly, when a small number of coils is employed, the RMSE of the GROGgridded image as compared to the standard Shepp-Logan phantom indicates that additional noise enhancement is present. However, as the number of coil elements increases, the RMSE values of the images approach a constant value; with the simulated coil arrangement, there is no appreciable difference between the use of six and 20 elements. The images in the center and right of the inset were reconstructed using GROG with six and 12 channels, and appear virtually identical. Because most commercially available coil arrays have at least six elements (with the tendency toward larger number of elements increasing), GROG can be performed with these coils on clinical systems without concern for additional noise enhancement due to insufficient coil sensitivity variations. The GROG effective DCFs for fully- and undersampled radial data are compared with the Ram-Lak and undersampled DCF proposed by Pipe (26) in Fig. 5. The GROG effective DCFs are shown in as solid lines, and the standard DCFs as dotted lines. As can be seen in the figure, the effective DCFs of GROG closely follow the standard radial DCFs, especially in the center of k-space, where the data weighting is essential to account for correlations between data points. It is important to note that the effective GROG DCF does not change after the data have become more sparse than dictated by the Nyquist criterion, as evidenced by the R ⫽ 4 undersampled data (the arrow in Fig. 5 indicates the point at which the Nyquist criterion is no longer fulfilled). This is consistent with the idea of the DCF as a tool for normalizing data correlation; as soon as the data no longer fulfill the Nyquist criterion in the azimuthal direction, compensation for varying density is no longer necessary. The image resulting from the rosette trajectory data gridded using GROG is shown in Fig. 6. The RMSE of this image in comparison to the noiseless Shepp-Logan phantom was calculated to be 0.17%, and the images are visually indistinguishable. It is important to note that the rosette data were gridded without a DCF or any additional gridding parameters, which is not possible with most other gridding techniques. The results of the gridding of the in vivo radial, spiral, and BLADE trajectories are shown in Fig. 7. In this figure, the results of convolution gridding are shown on the top, and the GROG images on the bottom. Upon visual inspection of the results, GROG and convolution gridding yield equivalent results, where any contrast or resolution differences are not evident. DISCUSSION The GROG has been demonstrated for radial, spiral, and BLADE trajectories. One advantage of GROG is that no precalculated DCFs or other parameters are required for gridding, whereas convolution gridding requires a DCF in addition to other parameters. While this is not a difficulty for the radial or spiral trajectory, the ability to grid BLADE, rosette, or stochastic data without having to calculate a DCF is highly advantageous. It is important to note that the effective DCF used in GROG, i.e., the averaging of shifted points that map to the same Cartesian location, cannot be used for other gridding techniques, because GROG explicitly calculates the values of the Cartesian points. Thus, after applying the appropriate GROG weights to each nonCartesian point, the resulting dataset is made up of purely Cartesian points, which can simply be averaged. It would also be possible to weight the shifted points with scaling factors that depend on the distances of their GROG shifts, although this method of calculating the DCF has not been examined. In addition, for undersampled datasets, GROG automatically performs an approximation of the high SNR, low artifact energy DCF proposed by Pipe (26) for undersampled datasets. Thus, undersampled data are also correctly density-compensated without the need for considerations about the degree of undersampling present in the dataset. Despite the fact that no additional parameters are needed for GROG, the in vivo GROG images presented here are visually indistinguishable from the convolution gridding images in terms of contrast and resolution. GROG is also less computationally intensive than convolutionbased methods. For a kernel size of five, convolution gridding must process 25 gridding points for each actual nonCartesian data point acquired. GROG, however, uses in essence a kernel size of one, which means that the gridding GRAPPA Operator Gridding 1263 FIG. 7. In vivo examples of images regridded with the standard convolution gridding (top) and with GROG (bottom) for radial (left), spiral (center), and BLADE (right) trajectories. operation can be accomplished in much less time and with fewer computational requirements. There are other methods besides GROG that perform data gridding without the need for a DCF. For instance, Uniform Resampling/Block Uniform Resampling (URS/ BURS) (10) is a method that performs data resampling by transforming the gridding problem into a linear equation which can be solved using singular value decomposition (SVD). As in GROG, no subsampling is employed for the gridding process and no DCF is required. However, this family of methods has several drawbacks. In the URS method, the large number of data samples leads to an inconveniently large linear equation. The BURS method is somewhat more practical, although it is highly sensitive to noise due to the need for a matrix inversion in the SVD. The extension to these approaches, regularized block uniform resampling (rBURS) (11), addresses this noise sensitivity problem, although the results are strongly dependent on the parameterization of the matrix inversion problem, i.e., the regularization and the size of the region of support. In comparison with GROG, which requires no parameterization, the URS family is much more difficult to employ, and is more computationally intensive. Several iterative methods have also been proposed that operate without a DCF. One such approach is iterative next-neighbor regridding (INNG) (12), which differs from convolution gridding in that a multiplication in the frequency domain is substituted for a convolution in the time domain. In this way, an effective sinc convolution can be performed in the time domain (through the multiplication of a box-car function with the image in the frequency domain) instead of using an approximation (such as a Kaiser-Bessel window). However, a number of iterations of the INNG algorithm are required before an artifact-free image results from the non-Cartesian data. Each iteration involves one Fourier transformation and one inverse Fourier transformation, which becomes quite time-consuming when large oversampling factors are used. A second iterative gridding algorithm is deconvolution-interpolation gridding (DING) (13), which also formulates the gridding problem as a linear equation. Unlike the URS family, this linear equation is solved using a conjugate gradient optimization. Because no subsampling is used, this method is less computationally intensive than INNG, although each iteration requires the same number of Fourier transformations. However, the convolution window used in DING must be optimized, again leading to potential difficulties in parameterization. The advantage of GROG over such methods is that GROG is a direct method, and requires no iteration or parameterization (i.e., subsampling, iteration numbers, window widths, or window shapes). In addition, all of the methods described above assume that the nonCartesian data fulfill the Nyquist criterion, which is not an assumption of GROG. However, it is important to note that unlike convolution gridding, URS/BURS/rBURS, INNG, or DING, GROG can only be used to grid data acquired with a multielement coil; this is not a requirement for most other gridding methods. For datasets with low SNR, the use of GROG without interpolation may increase the noise in the final images as compared with convolution gridding, as shown in Fig. 2. This SNR loss can be attributed to the application of weight sets to the individual noisy points, which amplifies the noise in the shifted point. However, if the data are first interpolated along the read direction to increase the read oversampling, more non-Cartesian points can be used to reconstruct a single Cartesian point. The use of more nonCartesian points, which must be shifted in different directions, acts as an effective averaging, which combats the noise enhancement. As demonstrated in Fig. 2, interpolation of the noisy data to an oversampling factor of four before gridding with GROG leads to images with RMSE values which do not differ greatly from those of the images resulting from convolution gridding. Thus, noise enhancement is not observed when using GROG to grid non-Cartesian data collected with clinical multielement coils, even in datasets with lower SNR. 1264 To calculate the base weights used in GROG, a calibration dataset is required. For some non-Cartesian trajectories, such as BLADE, this calibration dataset is automatically acquired with the rest of the data. However, for many others, additional data must be acquired. A method of obtaining the weights from the non-Cartesian data themselves is currently a work-in-progress (29). In addition, the non-Cartesian data must be acquired with a coil array that offers sufficient sensitivity variation to allow the calculation of robust weights. As demonstrated in the examples shown here, standard coils, which are readily available in a clinical setting, possess the sensitivity variations needed to reconstruct images using GROG. While the 18-channel abdomen coil used to acquire the BLADE data may not be currently commonplace, the 12-channel and eight-channel head coils used to acquire the radial and spiral data, respectively, are both standard and commercially available. Thus, the need for multichannel data from appropriate coils is not a major obstacle to use of GROG. An important feature of GROG is that this method can be used to grid undersampled data as well as fully-sampled data, which is not the case in convolution gridding. Convolution gridding methods assume that the non-Cartesian data fulfill the Nyquist criterion throughout the trajectory. Despite acceptable streak artifacts for angularly undersampled radial trajectories, gridding with convolutionbased methods essentially “fills up” more k-space locations than were actually measured. In the extreme case of a single point, GROG can be used to shift this point to the nearest Cartesian location, yielding a single point in kspace surrounded by zeros, a feat that cannot be accomplished with convolution gridding. The ability of GROG to grid undersampled non-Cartesian datasets to yield Cartesian data points near the sampled locations and zeros in all other locations could even be advantageous for other types of non-Cartesian parallel imaging reconstructions. For instance, non-Cartesian GRAPPA algorithms (30 –32) first reconstruct missing data points, and then grid the reconstructed data. The undersampled data can instead first be gridded with GROG followed by a Cartesian GRAPPA reconstruction of the resultant Cartesian data points. Preliminary results (33,34) show that this approach could greatly simplify non-Cartesian GRAPPA reconstructions. Similarly, the ability of GROG to shift acquired points by small amounts could be advantageous for a number of other applications (35). Finally, the GRAPPA operator gridding approach is less burdensome than other possible parallel imaging-based gridding techniques. GROG requires only one weight set for each orthogonal direction in k-space, i.e., 2D data require two weight sets. Using the method described in the first section of this article, the specific weight sets needed to shift the non-Cartesian points to the nearest Cartesian locations are then calculated out of these two weight sets. Thus, both the calibration and reconstruction processes are fast due to the small number of weight sets. In contrast, one could imagine using a technique such as parallel magnetic resonance imaging with adaptive radius in k-space (PARS) (21,22) for data gridding. PARS examines the local vicinity of a missing point to find possible source points, and then using these source points and a sensitivity map, calculates a weight set that allows one to reconstruct the Seiberlich et al. missing point. PARS is clearly inferior to GROG in terms of computational complexity, due to the repeated search operations. In addition, because GROG requires such a small number of weight sets, these weights can be calculated from simple Cartesian calibration data. Thus, GROG does not require coil maps, weight interpolation, or differentlyspaced Cartesian calibration data, which simplifies the reconstruction with respect to other possible parallel imaging-based gridding schemes. CONCLUSIONS The GROG technique proposed here is an alternative method for the gridding of non-Cartesian datasets. Instead of employing a convolution window as in the gold-standard gridding, the GRAPPA operator is instead used to shift the non-Cartesian data points to their nearest Cartesian locations. The calculation and application of the GROG weights is computationally efficient, as only one weight set is needed per orthogonal direction in k-space. Gridding of data acquired using radial, spiral, rosette, and BLADE trajectories has been demonstrated, where visually equivalent contrast and resolution were obtained with GROG and standard convolution gridding. Simulations were conducted which show that noise enhancement resulting from the use of GROG is minimal when a 1D interpolation is performed to increase the oversampling in the read direction. Only the trajectory, a small Cartesian calibration dataset and multichannel data are needed to perform GROG. This simplicity is advantageous for trajectories, such as BLADE, rosette, or stochastic sampling, that would otherwise require the calculation of a nontrivial DCF prior to convolution gridding. Additionally, single points or undersampled data can also be gridded and properly density-compensated with GROG. In summary, GROG demonstrates a novel use of parallel imaging to perform tasks other than the reconstruction of undersampled data. REFERENCES 1. Ahn CB, Kim JH, Cho ZH. High-speed spiral-scan echo planar NMR imaging. IEEE Trans Med Imaging 1986;5:2–7. 2. Meyer CH, Hu BS, Nishimura DG, Macovski A. Fast spiral coronary artery imaging. Magn Reson Med 1992;28:202–213 3. Lauterbur PC. Image formation by induced local interactions: examples employing nuclear magnetic resonance. Nature 1973;242:190 –191. 4. Glover GH, Pauly JM. Projection reconstruction techniques for reduction of motion effects in MRI. Magn Reson Med 1992;28:275–289. 5. Pipe JG. Motion-correction with PROPELLER MRI: application to head motion and free-breathing cardiac imaging. Magn Reson Med 1999;42: 963–969. 6. Noll, DC. Multishot rosette trajectories for spectrally selective MR imaging. IEEE Trans Med Imaging 1997;16:372–377. 7. Scheffler K, Hennig J. Frequency resolved single-shot MR imaging using stochastic k-space trajectories. Magn Reson Med 1996;35:569 – 576. 8. 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. 9. Pipe JG, Menon P. Sampling density compensation in MRI: rationale and an iterative numerical solution. Magn Reson Med 1999;41:179 – 186. 10. Rosenfeld D. An optimal and efficient new gridding algorithm using singular value decomposition. Magn Reson Med 1998;40:14 –23. GRAPPA Operator Gridding 11. Rosenfeld D. New approach to gridding using regularization and estimation theory. Magn Reson Med 2002;48:193–202. 12. Moriguchi H, Duerk JL. Iterative next-neighbor regridding (INNG): improved reconstruction from nonuniformly sampled k-space data using rescaled matrices. Magn Reson Med 2004;51:343–352. 13. Gabr RE, Aksit P, Bottomley PA, Youssef AB, Kadah YM. Deconvolution-interpolation gridding (DING): accurate reconstruction for arbitrary k-space trajectories. Magn Reson Med 2006;56:1182–1191. 14. Sodickson DK, Manning WJ. Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays. Magn Reson Med 1997;38:591– 603. 15. Jakob PM, Griswold MA, Edelman RR, Sodickson DK. AUTO-SMASH: a self-calibrating technique for SMASH imaging. Simultaneous acquisition of spatial harmonics. MAGMA 1998;7:42–54. 16. Pruessmann KP, Weiger M, Scheidegger MB, and Boesiger P. SENSE: sensitivity encoding for fast MRI. Magn Reson Med 1999;42:952–962. 17. Kyriakos WE, Panych LP, Kacher DF, Westin CF, Bao SM, Mulkern RV, Jolesz FA. Sensitivity profiles from an array of coils for encoding and reconstruction in parallel (SPACE RIP). Magn Reson Med 2000;44:301– 308. 18. Griswold MA, Jakob PM, Nittka M, Goldfarb JW, Haase A. Partially parallel imaging with localized sensitivities (PILS). Magn Reson Med 2000;44:602– 609. 19. Heidemann RM, Griswold MA, Haase A, Jakob PM. VD-AUTO-SMASH imaging. Magn Reson Med 2001;45:1066 –1074. 20. 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:1202–1210. 21. Yeh EN, McKenzie CA, Ohliger MA, Sodickson DK. Parallel magnetic resonance imaging with adaptive radius in k-space (PARS): constrained image reconstruction using k-space locality in radiofrequency coil encoded data. Magn Reson Med 2005;53:1383–1392. 22. Samsonov AA, Block WF, Arunachalam A, Field AS. Advances in locally constrained k-space-based parallel MRI. Magn Reson Med 2006; 55:431– 438. 23. Seiberlich N, Blaimer M, Barkauskas KJ, Breuer FA, Jakob PM, Griswold MA. Non-Cartesian data gridding using GRAPPA operator gridding (GROG). In: Proceedings of the 23rd Annual Meeting of the ESMRMB, Warsaw, Poland, 2006, 204 –205. Proc ESMRMB; 2006. p 300. 1265 24. Griswold MA, Blaimer M, Breuer FA, Heidemann RM, Muller M, Jakob JM. Parallel magnetic resonance imaging using the GRAPPA operator formalism. Magn Reson Med 2005;54:1553–1556. 25. Blaimer M, Breuer FA, Heidemann RM, Mueller M, Griswold MA, Jakob PM. Is parallel MRI without a prior information possible? In: Proceedings of the 12th Annual Meeting of ISMRM, Kyoto, Japan, 2004 (Abstract 2417). 26. Pipe JG. Reconstructing MR images from undersampled data: dataweighting considerations. Magn Reson Med 2000;43:867– 875. 27. Hoge RD, Kwan RK, Pike GB. Density compensation functions for spiral MRI. Magn Reson Med 1997;38:117–128. 28. Barkauskas KJ, Griswold MA, Sunshine JL, Duerk JL. Density compensation for trueFISP BLADE imaging in the steady-state (TrueBLISS). In: Proceedings of the 14th Annual Meeting of ISMRM, Seattle, WA, USA, 2006 (Abstract 2949). 29. Seiberlich N, Breuer FA, Blaimer M, Jakob PM, Griswold MA. Selfcalibrated GRAPPA operator gridder (SC-GROG). In: Proceedings of the 15th Annual Meeting of ISMRM, Berlin, Germany, 2007 (Abstract 153). 30. Griswold M, Heidemann RM, Jakob PM. Direct parallel imaging reconstruction of radially sampled data using GRAPPA with relative shifts. In: Proceedings of the 11th Annual Meeting of ISMRM, Toronto, Ontario, Canada, 2003 (Abstract 2349). 31. Heidemann RM, Griswold MA, Seiberlich N, Kruger G, Kannengiesser SA, Kiefer B, Wiggins G, Wald LL, Jakob PM. Direction parallel imaging reconstructions for spiral trajectories using GRAPPA. Magn Reson Med 2006;56:317–326. 32. Heberlein K, Hu X. Auto-calibrated parallel spiral imaging. Magn Reson Med 2006;55:619 – 625. 33. Seiberlich N, Heidemann RM, Breuer FA, Blaimer M, Griswold MA, Jakob PM. Pseudo-Cartesian GRAPPA reconstruction of undersampled non-Cartesian data. In: Proceedings of the 14th Annual Meeting of ISMRM, Seattle, WA, USA, 2006 (Abstract 2436). 34. Seiberlich N, Breuer FA, Heidemann RM, Blaimer M, Griswold MA, Jakob PM. Reconstruction of arbitrary non-Cartesian trajectories using pseudo-Cartesian GRAPPA in conjunction with GRAPPA operator gridding (GROG). In: Proceedings of the 15th Annual Meeting of ISMRM, Berlin, Germany, 2007 (Abstract 334). 35. Seiberlich N, Breuer FA, Moriguchi H, Jakob PM, Griswold MA. Fully autocalibrated parallel imaging for arbitrary trajectories using a combination of GRAPPA operator gridding and conjugate-gradient optimization. In: Proceedings of the 15th Annual Meeting of ISMRM, Berlin, Germany, 2007 (Abstract 1742).