An efficient iterative gridding reconstruction method with correction of off-resonance artifacts was developed, which is especially tailored for multiple-shot non-Cartesian imaging. The novelty of the method lies in that the transformation matrix for gridding (T) was constructed as the convolution of two sparse matrices, among which the former is determined by the sampling interval and the spatial distribution of the offresonance frequencies and the latter by the sampling trajectory and the target grid in the Cartesian space. The resulting T matrix is also sparse and can be solved efficiently with the iterative conjugate gradient algorithm. It was shown that, with the proposed method, the reconstruction speed in multipleshot non-Cartesian imaging can be improved significantly while retaining high reconstruction fidelity. More important, the method proposed allows tradeoff between the accuracy and the computation time of reconstruction, making customization of the use of such a method in different applications possible. The performance of the proposed method was demonstrated by numerical simulation and multiple-shot spiral imaging on rat brain at 4. Key words: non-Cartesian imaging; spiral imaging; gridding; reconstruction; off-resonance artifacts Ultrafast MRI methods employing non-Cartesian sampling schemes, such as spiral (1,2) and rosette scans (3,4), have been used in various applications. Image reconstruction from non-Cartesian k-space data, however, can be problematic in terms of speed and accuracy. The former problem arises from the fact that fast Fourier transformation (FFT) is not directly applicable (5-8), and the latter is because ultrafast imaging methods employing lengthy sampling durations are much more susceptible to the so-called off-resonance artifacts (ORA), relative to the conventional imaging methods (9).If FFT is to be used for image reconstruction, the nonCartesian k-space data need first be mapped onto a uniformly sampled k-space (i.e., gridding or resampling) (5-8). Gridding is commonly done by convoluting the acquired k-space data with a kernel function, followed by resampling in the Cartesian k-space (5). With this approach, correction of ORA can be carried out while gridding with the conjugate phase reconstruction methods (10-13). However, in the convolution-based gridding methods, the choices of the kernel size and the pre-and postcompensation functions may affect the outcome of reconstruction (5). In addition, the conjugate phase reconstruction methods for correcting ORA require field inhomogeneities varying smoothly across the image space (14) but this does not always hold in practical applications.Alternatively, gridding can be done by formulating a transformation matrix between the acquired k-space data and the target k-space data in the Cartesian space and solving it with single value decomposition (6) or iterative algorithms (7,8). With such approaches, subjective selection of the kernel and the pre-and postcompensation functions is no longer required. In practice, however, thes...