Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Magn Reson Med. Author manuscript; available in PMC 2014 April 14.
Published in final edited form as:
PMCID: PMC3985775

Parallel Imaging Reconstruction for Arbitrary Trajectories Using k-Space Sparse Matrices (kSPA)


Although the concept of receiving MR signal using multiple coils simultaneously has been known for over two decades, the technique has only recently become clinically available as a result of the development of several effective parallel imaging reconstruction algorithms. Despite the success of these algorithms, it remains a challenge in many applications to rapidly and reliably reconstruct an image from partially-acquired general non-Cartesian k-space data. Such applications include, for example, three-dimensional (3D) imaging, functional MRI (fMRI), perfusion-weighted imaging, and diffusion tensor imaging (DTI), in which a large number of images have to be reconstructed. In this work, a systematic k-space–based reconstruction algorithm based on k-space sparse matrices (kSPA) is introduced. This algorithm formulates the image reconstruction problem as a system of sparse linear equations in k-space. The inversion of this system of equations is achieved by computing a sparse approximate inverse matrix. The algorithm is demonstrated using both simulated and in vivo data, and the resulting image quality is comparable to that of the iterative sensitivity encoding (SENSE) algorithm. The kSPA algorithm is noniterative and the computed sparse approximate inverse can be applied repetitively to reconstruct all subsequent images. This algorithm, therefore, is particularly suitable for the aforementioned applications.

Keywords: magnetic resonance imaging, parallel imaging, kSPA, sparse matrix, sparse approximate inverse, functional imaging, massive parallel imaging

Parallel magnetic resonance imaging (MRI) utilizes multiple coils to simultaneously receive radio frequency (RF) signals emitted from a scan subject (14). Aided by the spatial distribution of the RF coils’ reception sensitivity, parallel MRI has been widely used for improving imaging speed or reducing artifacts. A key problem in developing parallel MRI has been the computational difficulty in forming an image from data acquired by multiple coils. Over the past several years, two feasible classes of algorithm have been introduced for parallel imaging reconstruction: image-domain algorithms and k-space algorithms. Image-domain algorithms directly compute the image from the k-space data acquired by each coil, where as the k-space algorithms compute the spectrum of the image.

The most successful image-domain algorithm has been the iterative sensitivity encoding (SENSE) algorithm for arbitrary trajectories by Pruessmann et al. (5,6). Another image-based algorithm named “sensitivity profiles from an array of coils for encoding and reconstruction in parallel” (SPACE RIP) has also been introduced by Kyriakos et al (7) to reconstruct an image column-by-column along the phase encoding direction. k-Space algorithms include, for example, simultaneous acquisition of spatial harmonics (SMASH) by Sodickson and Manning (8) and generalized GRAPPA by Bydder et al. (9), generalized autocalibrating partially parallel acquisitions (GRAPPA) by Griswold et al. (10,11), and parallel imaging with adaptive radius in k-space (PARS) by Yeh et al. (12), among other methods (13,14). The iterative SENSE algorithm expressed the k-space data as a linear combination of the spatially-encoded magnetization in which the multiplicative spatial encoding function is the product of coil sensitivity and the Fourier encoding function. The resulting system of linear equations is solved with the conjugate gradient (CG) or related method in an iterative fashion since the sheer size of the design matrix precludes a direct solution of the inverse problem (15). To improve the reconstruction speed, the matrix-vector multiplication encountered in each iteration is replaced by a gridding and inverse gridding procedure, which is an important innovation that has made iterative SENSE feasible for arbitrary sampling trajectories with a reasonable reconstruction speed (6). With the perfect knowledge of coil sensitivity, the iterative SENSE algorithm is shown to provide an accurate estimation of the true underlying image.

GRAPPA was first introduced to reconstruct an image that is partially acquired on a Cartesian grid (10). The essential idea is to estimate the missing k-space data points by linearly combining its acquired neighboring points. The weights used in this estimation are first trained on some calibration lines that are typically acquired near the center of the k-space. In addition, there have been some recent efforts to extend this method beyond Cartesian trajectory (16,17). The PARS algorithm estimates the k-space data on a grid using its neighboring data points sampled on arbitrary trajectories similar to the gridding procedure. Contrary to GRAPPA, the combination weights of PARS are calculated using coil sensitivity maps rather than the calibration lines. Both GRAPPA and PARS compute each individual coil image first and combine the resulting coil images through a sum-of-squares reconstruction.

Despite the successes of these existing algorithms, it remains a challenge to rapidly and reliably reconstruct an image from undersampled k-space data. One major issue is the difficulty in processing the large amount of data generated by higher spatial resolution, more coil elements, and 3D imaging. In addition, functional studies, such as functional MRI (fMRI) (1820), perfusion imaging (21,22), and diffusion tensor imaging (DTI) (2327), pose an especially difficult challenge, in which thousands of images have to be reconstructed for a single study. For example, a typical whole-brain DTI study generates on the order of 1000 images. Without parallel computing implementation, it typically takes the iterative SENSE algorithm minutes to reconstruct a 256 × 256 image acquired with an eight-channel receiving coil. Even with a speed of one minute per image, it takes over 16 hours to reconstruct every 1000 images. Such long computational time makes parallel imaging essentially impractical for routine clinical functional studies. For dynamic studies or DTI studies, in principle, the design matrix only needs to be inverted once and then can be applied repetitively to reconstruct all subsequent images. However, iterative reconstruction has precluded such approach. Alternative means that allow the precomputing of the complex weights needed for parallel imaging reconstruction are, therefore, highly desirable for these types of applications.

In this work, a systematic parallel imaging reconstruction algorithm in k-space, termed k-space sparse matrices (kSPA), is proposed that particularly suits such kind of repetitive image reconstruction process. The image reconstruction problem is formulated in the k-space as a linear algebra problem. The kSPA algorithm then solves the problem by taking advantage of the sparsity of the resulting matrices. With this new algorithm, a sparse approximate reconstruction matrix (i.e., the inverse of the design matrix) can be directly computed. Multiple images can then be reconstructed in a fraction of the time needed for iterative SENSE by repetitively applying the reconstruction matrix through a matrix and vector multiplication. This algorithm is demonstrated with both simulated and in vivo data.


Because of the discrete nature of MR image acquisition, all mathematical treatment of the reconstruction algorithm involves only discrete samples of continuous functions.

Parallel Imaging Reconstruction in k-Space

Assuming that the n-th coil has a receiving sensitivity of sn(rρ) with a Fourier transform of sn(kρ) on a Cartesian grid (ρ = 1, …, N2 for an N × N grid), then the signal received by the n-th coil in the k-space can be approximated as,


Here, m(kρ) represents the true magnetization to be imaged in the k-space, and the sign “*” represents the two-dimensional (2D) convolution.

Because the field of view (FOV) of an imaging experiment is always limited, MR data in any k-space location can be calculated through interpolation following the Nyquist theorem. Specifically, for an arbitrary location κμ in k-space, the data acquired by the n-th coil can be written as,


where c(kκ) is the interpolation kernel. Ideally, c(kκ) should be a sinc function of infinite length. In practice, c(kκ) is usually implemented as a finite-width kernel, for example, a Kaiser-Bessel window (28). To differentiate a Cartesian grid from an arbitrary sampling location, the Latin letter k is used to indicate a Cartesian grid point, while the Greek letter κ is used to indicate an arbitrary location.

Combining Eqs. [1] and [2], the data on an arbitrary k-space location can be written as,


where sn(κμ) is the spectrum of the coil sensitivity inter-polated at an arbitrary k-space location κμ according to:


With multiple receiving coils and a number of sampling locations, Eq. [3] forms a system of linear equations that can be denoted as,


Here, d is a column vector stacked with the k-space data acquired by all coils; m is also a column vector with the k-space value to be estimated; G is the coefficient matrix. Specifically,


Here, nκ denotes the total number sampling locations in k-space and nc denotes the total number of receiving coils.

In principle, Eq. [5] can be solved with various techniques for solving systems of linear equations, for example, through computing the pseudoinversion of matrix G. However, the typical size of matrix G is prohibitively large and prevents such kind of matrix operation on a typical personal computer. For example, with an eight-channel receiving coil, a reduction factor of 2, and an image size of 256 × 256, the size of the G matrix will be around 260,000 × 65,536. To store this matrix requires 128 GB of memory with double floating-point precision.

To reduce the storage requirement, we study the sparsity of the matrix G by exploring the physical properties of the coil sensitivity. According to the Biot-Savart law, the coil sensitivity decays roughly quadratically with the distance from the coil location. Although the actual decaying speed depends on the size and shape of the coil, in general, when the object is sufficiently far away from the coil, the sensitivity is a smooth function containing only low spatial frequency components. In other words, the convolution kernel defined by the coil sensitivity is very compact leading to a sparse matrix G (Fig. 1a).

FIG. 1
Illustration of sparse matrices: (a) an example of the sparse design matrix G resulting from an echo-planar imaging experiment with an image size of 32 × 32, an eight-channel coil, and a reduction factor of 4; (b) a schematic illustration of the ...

Let ws be the cutoff bandwidth (BW) of the sensitivity, i.e.,


Here, ||·||2 denotes the vector Euclidean norm. In this work, ws was chosen such that at the cutoff frequency sn(kρ) decreases to around 0.36% of its peak value. With ws = 8π/128 (i.e., eight pixels in a 256 × 256 matrix), it requires approximately 700 MB instead of 128 GB of memory to allocate the aforementioned matrix G, and the storage decreases linearly with increasing reduction factor. With this assumption, Eq. [5] results in a set of sparse linear equations, which can be solved , for example, by the LSQR algorithm introduced by Paige and Saunders (29). Nevertheless, the application of the LSQR algorithm is not the focus of this work and no further description of the implementation will be given.

The kSPA Algorithm

Although LSQR or any other suitable iterative method can provide a very accurate solution to Eq. [5], iterative methods are time consuming for repetitive image reconstruction, such as fMRI or DTI, since multiple iterations are generally required for each image. Ideally, one would like to find a generalized inverse matrix G+ applicable to all images such that


where I is the identity matrix. Note that, since G usually is an overdetermined system, it does not have a right inverse G+ such that


A least-squares solution to Eq. [10] is given by


Here, GH refers to the complex conjugate transpose of G. Unfortunately, G+ can not be computed directly for such large systems. Even if G+ can be computed in a reasonable amount of time, it is impractical to store G+ in the memory or even on the hard drive. Therefore, approximation methods are necessary in order to compute and store G+.

Our method is to approximate G+ with a sparse matrix. The concept of sparse approximate inverse was first introduced by Benson and further developed by other authors as a preconditioner for the CG method (3032). In Benson’s original work (32), the matrix to be inverted has been assumed to be real and square. In our case, G is complex-valued and overdetermined. We propose to first find a sparse approximate inverse of (GHG)+ that is a symmetric square matrix. Let M = (GHG), then M is a N2 × N2 symmetric matrix whose entries are given by


Note that M is a sparse matrix, since Mρρ′ is only nonzero when there exists at least one sampling location κμ such that both sn(κμ-kρ) and sn(κμkρ′) are nonzero. More specifically, given a kρ, Mρρ′, is only nonzero for those kρ′ such that ||kρkρ′||2 ≤ 2w, (Fig. 1b).

A sparse approximate inverse M+ can be found by solving the following minimization problem:


subject to constraints on the sparsity pattern of M+, i.e., the number and position of the nonzero entries of M+. Here, ||·||F denotes the Frobenius norm of a matrix (33). For practical implementation, this minimization problem can be divided into N2 independent least squares problems and M+ can be solved for row-by-row.

To solve Eq. [14], a sparsity pattern needs to be first prescribed for M+. The sparse pattern can be determined by studying the analytical expression of the inverse of a square matrix. Following the Cayley-Hamilton theorem (33), if a rank-n square matrix M is nonsingular, then M−1 can be expressed as a linear combination of the powers of M, with the order of powers ranging from 0 to n − 1. If M is scaled properly such that ||M||F < 1, then M−1 can be approximated with some low-order powers of M. The physical significance of this low-order approximation means that a k-space sample is only affected by its neighboring samples within a certain distance. Given such a distance w, for each row of M+, Eq. [14] can be written explicitly as,




This equation can be solved for each row of M+. For the ρth row, ρ′ is the set of k-space grid points such that ||kρkρ′||2w; μ is the set of sampling locations such that both sn(κμ-kρ) and sn(κμkσ) are nonzero for a given pair of kρ′ and kσ; δ(·) is the Kronecker delta function. For the given ρ-th row, the range of kσ is given by ||kρkσ||2wσ with wσw + 2ws (Fig. 1c). Note that although the ranges of sn(kρ) and M+ are both drawn circularly in Fig. 1b and c, in practice, a square range can be used for simplicity.

Equation [15] can be further modified to include the apodization correction needed to address the effect of the convolution kernel:


An additional advantage of using the convolution kernel instead of the delta function in the right hand side of the equation is to reduce the ringing artifact caused by the sharp transition of the delta function. Equation [16] forms a relatively small set of linear equations with (w + 1)2 unknowns and (wσ + 1)2 equations, where wσ defines the maximum distance between kρ and kσ. In a matrix format, this set of equations can be written as


The column vector mρ+ is the stack of Mρρ+; Mρ is the coefficient matrix, which basically is the submatrix of M corresponding to mρ+; and cρ is the column vector formed by c(kρkσ). For the above equation to be overdetermined, it is required that wσ > w. In this work, wσ is chosen to be w + ws.

The set of equations defined by Eq. [17] is solved for each row of M+ by computing the pseudoinverse of Mρ with the truncated singular value decomposition (SVD) method (33):


The threshold of the singular value was chosen to be the maximum singular value multiplied by the size of the Cartesian grid (i.e., N) and the floating point precision. This threshold was relatively conservative, representing only the level of quantization error. A larger threshold may be chosen to further suppress noise and accelerate the computation. For each ρ, mρ+ yields the ρ-th row of M+. After all rows of M+ are computed, G+ is calculated following Eq. [12]. Specifically, each element of G+ is computed as,


The summation is over those ρ′ such that ||kρkρ′||2w.

The above proposed algorithm for k-space parallel imaging reconstruction can be summarized as:

kSPA Algorithm

  1. For each coil, compute its sensitivity in the k-space on a ws × ws Cartesian grid.
  2. For the ρ-th row of M+, find all kρ′ such that ||kρkρ′||2w and all kσ such that ||kρkσ||2wσ with wσ = w + ws. Solve Eq. [16] for Mρρ+ with SVD.
  3. For the ρ-th row of G+, find all κμ such that ||κμkρ||2ws + w. Compute Gρ(μ,n)+ with Eq. [19].
  4. Repeat step 2 and step 3 until all rows of G+ are computed.
  5. Compute the k-space data points on the Cartesian grid with m = G+d.
  6. Multiply m with a Fermi window (etc.) and Fourier transform m to obtain an image.

Note that if there is sufficient memory, sn(κμkρ) can be precomputed and stored in the matrix G. This memory requirement for storing G usually can be easily satisfied. By directly accessing the matrix elements rather than computing it every time, steps 2 and 3 can be further sped up.

Fast Implementation of kSPA

The main computational burden of kSPA lies in the computation of M+. For an arbitrary sampling pattern, because the summation over μ is different, the set of equations defined by Eq. [16] is different for different row of M+. As a result, Eq. [16] needs to be solved independently for every row of M+. However, if the sampling pattern has a certain shift-invariant property, then the amount of computation can be dramatically reduced. For example, for an undersampled Cartesian trajectory where R phase encoding lines are skipped, the sampling pattern is periodic with a period of R grid points. Because of this self-resembling sampling pattern, there are only R sets of unique system of linear equations excluding some edge points (discussed at the end of this section). As a result, Eq. [16] only needs to be solved for R consecutive sampling points on a representative line along the phase encoding direction. These R sets of weights can then be selectively used by other grid points based on their relative locations. By solving only R sets of equations, the amount of computation is reduced by a factor of N2/R. For example, with N = 256 and R = 4, the speed-up factor is 4096. Of course, for Cartesian trajectory, aliasing can also simply be unfolded in the image domain pixel by pixel.

For an arbitrary non-Cartesian trajectory, the sampling pattern is not periodic; therefore the above implementation based on the periodicity of the sampling pattern is no longer feasible. For arbitrary trajectories, we use a block-by-block implementation of kSPA that significantly reduces the total amount of computation required. We first observe that, although the matrix Mρ defined in Eqs. [16] and [17] is different at different k-space grid point, this matrix does not differ significantly for two adjacent grid points. Specifically, given two neighboring points kρ1 and kρ2, the set of kρ′ satisfying ||kρ1kρ′||2w and the set of kρ′ satisfying ||kρ2kρ′||2w consist of mostly the same grid points except for a few points near the edge of the circle. Based on this observation, we can divide the k-space into a number of blocks (Fig. 2). Each block has a half-width of wb with wb < w. For each block, the matrix Mρ is computed based on the center point of the block as indicated by the small white circle in Fig. 2. The pseudo-inverse of Mρ is computed and the result is applied for every grid point inside the block to compute mρ+ following Eq. [18]. By computing the pseudoinverse only once for the whole block rather than computing one for each grid point, the amount of computation required for inverting Mρ is reduced by a factor of (2wb + 1)2. Because the matrix-vector multiplication in Eq. [18] is relatively inexpensive, the total speed-up factor is approximately (2wb + 1)2, which equals to 529 with wb = 11.

FIG. 2
An illustration of the block-wise implementation of kSPA. The solid straight lines indicate the Cartesian grid, which is divided into small blocks shown in dotted squares. Each block has a half-width of wb. For each block, Mρ is constructed once ...

It should be pointed out that, for grid points within a distance of w from the edge of the matrix, mρ+ cannot be computed accurately using the block-wise implementation. This limitation is a result of the fact that there are no sampling points outside the matrix. For example, a grid point on the left edge of the matrix has only sampling points on its right hand side, while a grid point in the middle of the matrix has sampling points on both sides. The resulting equations defined by Eq. [16] are vastly different. Therefore, in the block-wise implementation, the edge points need to be computed separately. Alternatively, the block-wise implementation can be applied to all grid points, while the synthesized data towards the edge has to be filtered out using a low-pass filter, for example using a Fermi filter. For simplicity, the latter approach is implemented in this work.

Sensitivity Estimation

The sensitivity information required for kSPA reconstruction can be obtained with any established sensitivity estimation method (5,3436). For example, the sensitivity can be obtained with an additional body coil image that is assumed to have homogeneous sensitivity. In the absence of a body coil image, the sum of squares image can also be used for sensitivity estimation (Fig. 3a). Another common method is to use a low resolution image that is either acquired with an extra navigator or with a self-navigated trajectory, for example a variable density spiral (37). The self-navigation approach avoids the need for a separate calibration scan and errors from misregistration.

FIG. 3
The estimation of coil sensitivity map using 2D thin-plate spline fitting. (a) An example of raw coil sensitivity map calculated by dividing the coil image with the sum-of-squares image. This raw sensitivity map is typically noisy. (b) The raw sensitivity ...

Because the kSPA algorithm utilizes only the low-frequency component of the coil-sensitivity, it is important to preprocess the raw coil sensitivity and remove high-frequency noise and residual tissue contrast (Fig. 3a). Any signal void or noise in the sensitivity map will contribute errors in the estimation of the low-frequency component of the coil-sensitivity. Here, we use the thin-plate spline function for smoothing the coil sensitivity because it not only provides sensitivity estimation for regions of signal voids inside the object, but also provides some sensitivity extrapolation outside the object. Given a set of n control points (xi,yi) on a 2D plane, the thin-plate spline function is given by,


Given a set of measured points on the sensitivity map, the smoothed real and imaginary part of the full sensitivity map can be estimated independently by solving the following minimization problem:


Here s(xi, yi) indicates either the real or the imaginary part of the measured sensitivity at control point (xi, yi); p is a smoothing parameter between 0 and 1; p = 0 corresponds to linear fitting and p = 1 corresponds to thin-plate spline interpolation. The optimal value of p can be determined according to the algorithm by Reinsch (38). For simplicity, we use a “good” empirical number, p = 0.0002, in this work. Matlab (The MathWorks, Inc.) provides an implementation of the described thin-plate spline smoothing method in its function call “tpaps.”

To remove any potential fast oscillating phase in the raw coil sensitivity that is common for all coils, the phase of the sensitivity map of a randomly chosen coil is subtracted from all coils. Each raw coil sensitivity map is masked based on a threshold set on the order of the noise level of the sum-of-squares image (Fig. 3b). A small set of support points is then randomly selected within the masked region and used to solve Eq. [22] (Fig. 3c). This set of control points consists of only 2.5% of the total points inside the masked region. As can be seen, the thin-plate spline function not only provides a good fitting to the support points within the object, but also provides a reasonable extrapolation beyond the boundary of the object (Fig. 3d–f).



The kSPA algorithm is applied for various k-space trajectories including a Cartesian trajectory, a spiral trajectory, and a random trajectory. All sampling trajectories correspond to an image matrix size of 128 × 128. The spiral trajectory was designed based on the analytic method by Glover (39). It contains eight interleaves. Each interleaf has 1861 sampling points. A Shepp-Logan phantom and an eight-channel receiving coil were used to simulate the k-space data via inverse gridding (6,27,40). The image seen by each coil was simulated by multiplying the phantom image with the corresponding coil sensitivity map. To reduce the effect of k-space circular convolution resulting from this image-domain multiplication and subsequent discrete Fourier transform, both the phantom image and coil sensitivity maps were upsampled by a factor of 2 prior to the multiplication of the phantom image with coil sensitivity maps. To improve the accuracy of inverse gridding, each coil image was further zero-padded by a factor of 2 before performing Fourier transform. In addition, the edges of the coil sensitivity maps that were outside the imaging object are tapered with a Fermi filter in order to further suppress the energy leakage caused by Gibbs ringing. Images were reconstructed using reduction factors ranging from 1 to 4.


In vivo brain images of a healthy volunteer were acquired using a spiral readout trajectory on a 1.5T whole-body system (GE Signa; GE Healthcare, Waukesha, WI, USA) equipped with a maximum gradient of 50 mT/m and a slew rate of 150 mT/m/s. An eight-channel head coil (MRI Devices Corporation, Pewaukee, WI, USA) was used for image acquisition. The scan parameters were: FOV = 24 cm, TR = 4 s, TE = 90 ms, BW = 125 kHz, and matrix size = 256 × 256. The spiral readout trajectory consists of 32 interleaves. Each interleaf has 1744 samples. Under-sampling in k-space was achieved by skipping a certain number of interleaves. Coil sensitivities were measured using a spiral-in navigator that samples the k-space on a 64 × 64 grid and precedes each spiral-out interleaf.

All image reconstructions were performed using Matlab (Version 7, Release 14; The MathWorks Inc., USA), running on a LINUX PC equipped with a 3.20 GHz Intel Xeon CPU and 5 GB random access memory (RAM). For the simulation, the kSPA reconstruction parameters were: ws = 12, w = 20, and block size wb = 13; for the in vivo study, ws = 8, w = 7, and wb = 5.



Figure 4 compares images reconstructed with gridding and kSPA for all three types of trajectories. Reduction factors range from 1 to 4 with one meaning no reduction. For each trajectory, the first row shows images reconstructed with a gridding algorithm (28) and the second row shows images reconstructed with kSPA. The kSPA parameters were ws = 12, w = 20, and block size wb = 13. These parameters were empirically optimized by varying their values as illustrated in the following paragraph. The total number of blocks is 22. As a result, Eq. [17] needs to be solved 22 times. As seen in Fig. 4, the kSPA algorithm results in excellent image quality for all three sampling trajectories. The corresponding relative computational time is listed in Table 1. Because the computational time varies with software implementation and hardware performance, the listed time is normalized by the iterative SENSE reconstruction time (1.3 min) for R = 1, hence, is for reference purpose only.

FIG. 4
kSPA reconstruction for three sampling trajectories with undersampling factors up to 4: (a) Cartesian, (b) spiral, and (c) random trajectory. The first row of each group shows images reconstructed with a simple gridding procedure. With k-space under-sampling ...
Table 1
Relative Reconstruction Time Required for Images Shown in Fig. 4*

To further assess the dependence of the kSPA reconstruction on the three window sizes, ws, w, and wb, the randomly sampled data with R = 4 were reconstructed using various combinations of the three parameters. A complete analysis of this relationship requires the variation of three independent variables. For the convenience of data presentation, three types of quantitative analysis were performed. The first analysis was performed by varying w while maintaining ws constant and wb = w/2; the second analysis was performed by varying ws while keeping both w and wb constant; the third analysis was performed by varying wb while keeping both ws and w constant. The resulting normalized image mean square errors (MSE) were plotted as a function of the corresponding variables (Fig. 5). Three sets of typical kSPA images were also shown in Fig. 6 to illustrate the image quality under different reconstruction parameters. In general, the image quality improves when increasing both ws and w. Severe residual artifacts remain in the image when insufficient window sizes are used.

FIG. 5
The mean square error (MSE) of kSPA reconstruction as a function of various parameters. MSE was computed using the simulated randomly sampled data with R = 4. a: MSE decreases with increasing w; (b) MSE decreases with increasing ws; and (c) MSE increases ...
FIG. 6
Some examples of images reconstructed at the marked positions on the curves in Fig. 5. a: From left to right: w = 5, 10, and 15; (b) from left to right: ws = 8, 12, and 18; (c) from left to right: wb = 3, 9, and 16. In general, larger w and ws improve ...


Figure 7 shows the in vivo results. The first row shows a typical coil image reconstructed with gridding for each reduction factor, while the second row shows the kSPA images. The kSPA parameters were ws = 8, w = 7, and block size wb = 5. The resulting sparse matrix M+ requires 208 MB storage space. As expected, the gridding-reconstructed images exhibit severe aliasing artifacts. However, such severe aliasing artifacts are not visible in the kSPA images. For comparison, images reconstructed with iterative SENSE (6) and the difference images between iterative SENSE and kSPA are also shown in Fig. 7. The same coil sensitivity maps were used for both kSPA and iterative SENSE reconstruction, and they were estimated using the described thin-plate spline fitting method. From Fig. 7, the image qualities are comparable for these two techniques. For reference, the relative reconstruction times are also listed in Table 1. Notice that all computations were implemented in Matlab as scripts, and no precompiled executable functions were used. These times are not intended for strict measures of the reconstruction speed of these two algorithms. Rather, they are shown here simply to illustrate the rough order of the computational time resulting from this particular Matlab implementation.

FIG. 7
In vivo kSPA reconstruction for spiral sampling with under-sampling factors up to 4. The first row shows a typical image from one coil element reconstructed with gridding. In the cases of un-dersampling, images are severely distorted. The second row shows ...


We have shown that kSPA is a k-space-based parallel imaging reconstruction algorithm that can be applied to arbitrary k-space sampling trajectories. While the iterative SENSE algorithm relates the image-domain data with the acquired k-space data through a sensitivity-encoded Fourier transform, kSPA expresses the acquired k-space data as a convolution between the spectrum of the coil sensitivity and the spectrum of the image. The kSPA reconstruction, therefore, is a deconvolution process based on undersampled data, which is accomplished by approximating the inverse of the design matrix with a sparse matrix. Its feasibility and accuracy is verified in the simulation study using various undersampling ratios and various trajectories including a Cartesian, a spiral, and a random trajectory. In vivo studies also showed that the image quality of kSPA is comparable to that of iterative SENSE with a slightly increased computational time for a single image. For repetitive image reconstruction, an important strength of kSPA is its significant increase in speed compared to the iterative SENSE approach, which becomes increasingly important for time-series data, DTI, and spectroscopic imaging.

Gridding Interpretation of kSPA

The image reconstruction process of kSPA offers an analogy to the traditional gridding procedure (28). Based on Eq. [19], the k-space data on a k-space grid point can be computed as


This equation can be interpreted as a double convolution. Each coil’s data are first gridded onto a Cartesian grid with a gridding kernel equaling to the complex conjugate of the sensitivity of that coil, i.e., sn(k). The gridded data from all coils are then summed together. Finally, the sum of all gridded data is convolved with each row of M+ to estimate each point of m(kρ). In general, M+ provides a spatially varying convolution kernel. This gridding interpretation is illustrated in Fig. 8.

FIG. 8
A schematic illustration of the kSPA algorithm shows a gridding analogy. From left to right, k-space data sampled on arbitrary trajectories (shown in small black circles) are first gridded onto Cartesian grids using a convolution kernel equaling to the ...

In the case of fully sampled single-coil acquisition, the gridding kernel becomes the conventional gridding kernel c(k) instead of sn(k), and M+ is simply the postsampling density compensation. More specifically, assuming the single-channel coil has homogeneous reception sensitivity, the Fourier transform of the sensitivity is a delta function. As a result, sn(k) equals to c(k) following Eq. [4]. If the coil has inhomogeneous reception sensitivity, in order to correct the sensitivity modulation, the gridding kernel has to be s*(k), that is, the spectrum of the coil sensitivity. When multiple coils are used, such gridding procedure can be carried out independently for each coil. However, if the k-space is undersampled, the resulting image after this gridding step is still aliased. The aliasing artifacts have to be removed by multiplying the gridded k-space data with the unaliasing matrix M+.

Effect of Kernel Width and Block Size

The quality of kSPA reconstruction depends on a number of parameters, namely the cutoff BW of coil sensitivity ws, the reconstruction kernel width w, and the block size wb. The choice of the sensitivity’s cutoff BW is determined by the physical design of the coil. The main reason of using a limited cutoff BW is to reduce the total amount of computation required in calculating the approximate inverse. As shown in Eq. [16], larger ws results in more summation operations when computing the convolution of the sensitivity, and, therefore, longer computational time. On the other hand, a sufficient width is necessary for accurately approximating the coil sensitivity and eventually removing the aliasing artifacts (Fig. 6b). In general, wider BW is required for larger sampling reduction factor, while the requirement can be relaxed for low reduction factors.

Given a sufficiently wide cutoff BW, the kSPA image quality is mainly determined by the reconstruction kernel width w. The reconstruction kernel width w limits the total number of nonzero entries in each row of the sparse approximate inverse M+, or, in other words, it decides how many neighboring points are used to synthesize each k-space data point. In principle, the larger w is, the more accurate the approximate inverse will be (Fig. 6a). In the extreme case, when both w and ws are equal to one-half of the imaging matrix size, the sparse approximate inverse M+ becomes the exact inverse of M in the least-squares sense (see Eq. [14]), and the kSPA image will be equivalent to the iterative SENSE image. However, the larger w is, the longer it takes to compute M+. Therefore, for fast image reconstruction, it is beneficial to use a minimum width for the reconstruction kernel without suffering severe image degradation. Mathematically, this minimum width is determined by the structure of the matrix M, which in turn is determined by coil design and k-space sampling pattern. Consequently, such a minimum width cannot be determined a priori. Nevertheless, a general guideline is provided by the Cayley-Hamilton theorem. This theorem expresses M−1 as a linear combination of powers of M. If only low-order powers are kept in the linear combination, then the reconstruction kernel should be w ≥ 2ws.

The block size wb, together with the reconstruction kernel width w are the two major parameters that affect the reconstruction speed. Clearly, the larger wb is, the faster the reconstruction is. However, to avoid artifacts caused by this block-wise implementation of kSPA, it is recommended that wb < w (Fig. 6c). As we have shown, the parameters used to reconstruct the in vivo data are significantly different from those used in simulation. This difference is mainly caused by differences in coil sensitivity used in simulating the k-space data.

In summary, the three parameters ws, w, and wb provide a tradeoff between image quality and reconstruction speed. Because kSPA is an approximate reconstruction method, residual artifacts may exist in images with high reduction factors (e.g., in the case of R = 4 in Fig. 4). However, such residual artifacts can be minimized by careful selecting the three reconstruction parameters (Figs. 5 and and6).6). When proper parameters are used, no significant difference between kSPA and SENSE is observed in the sensitivity to R.

Noniterative Reconstruction

The kSPA algorithm synthesizes k-space data directly. Similar to other k-space based algorithms, this process is noniterative. Such noniterative reconstruction offers some advantages in certain applications, such as repetitive image reconstruction.

Similar to GRAPPA and PARS, kSPA synthesizes k-space data using limited data support in the k-space. Particularly, kSPA shares some unique characteristics with PARS in that both compute the reconstruction weights purely based on the coil sensitivity maps instead of the measured raw data. In addition, by recognizing the compact structure of coil sensitivity’s spectrum, PARS and kSPA enable the computation of reconstruction weights by solving a set of small inverse subproblems and permit the storage and repetitive application of the reconstruction weights. However, kSPA differs fundamentally from PARS in the way the weights are computed and in the meaning of the weights. In PARS, the weights are the weighting factors that combine all coil sensitivity maps to synthesize a given coil sensitivity map. Accordingly, PARS reconstructs each coil image one-by-one. In kSPA, on the other hand, the weights are the inverse matrix of the autocorrelation function of the sampled spectrum of the coil sensitivity. kSPA reconstructs the image itself instead of the coil images. Because the weights only need to be computed once in kSPA, whereas a different set of weights have to be computed for different coils in PARS, kSPA offers a significant gain in speed, especially for situations in which there is a large number of coils.

Although kSPA is slightly slower than the iterative SENSE algorithm for a single image reconstruction in the example shown in Table 1, the result of kSPA reconstruction can be applied directly to reconstruct subsequent images through a simple matrix and vector multiplication, which can be accomplished in one second for this particular example. As a result, kSPA is particularly useful, for example, for fMRI and perfusion imaging, in which a large number of images are acquired in the same study. It is estimated that, for 1000 images, the image reconstruction time can be reduced by a factor of 100.


A general parallel imaging reconstruction algorithm termed kSPA has been introduced that is k-space–based and applies for arbitrary k-space sampling trajectories. This algorithm approximates the reconstruction in k-space by approximating the inverse of the design matrix with a sparse matrix. Such computation is noniterative and the results can be applied to reconstruct all subsequent images with the same prescription (e.g., slice location, FOV, resolution, and coil placement) acquired in the same imaging session. Both simulated and in vivo studies have shown that kSPA provides excellent image quality that is comparable to that of iterative SENSE. It is expected that kSPA will find great utility in various applications including fMRI, DTI, perfusion-weighted imaging, as well as spectroscopic imaging.


Grant sponsor: National Institutes of Health (NIH); Grant numbers: 1K99NS057943-01; 1R01NS35959; 1R01EB002771; Grant sponsor: Center of Advanced MR Technology of Stanford; Grant number: NCRR P41 RR 09784; Grant sponsor: Lucas Foundation.

We thank Michael Saunders, Ph.D., for discussions on LSQR, Murat Aksoy, M.S., for assisting the image acquisition, Huanzhou Yu, Ph.D., for discussions on sparse matrix, and Karen Chen, M.S., for her editorial assistance.


1. Carlson JW. An algorithm for NMR imaging reconstruction based on multiple RF receiver coils. J Magn Reson. 1987;74:376–380.
2. Hutchinson M, Raff U. Fast MRI data acquisition using multiple detectors. Magn Reson Med. 1988;6:87–91. [PubMed]
3. Roemer PB, Edelstein WA, Hayes CE, Souza SP, Mueller OM. The NMR phased array. Magn Reson Med. 1990;16:192–225. [PubMed]
4. Carlson JW, Minemura T. Imaging time reduction through multiple receiver coil data acquisition and image reconstruction. Magn Reson Med. 1993;29:681–687. [PubMed]
5. Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: sensitivity encoding for fast MRI. Magn Reson Med. 1999;42:952–962. [PubMed]
6. Pruessmann KP, Weiger M, Bornert P, Boesiger P. Advances in sensitivity encoding with arbitrary k-space trajectories. Magn Reson Med. 2001;46:638–651. [PubMed]
7. 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. [PubMed]
8. Sodickson DK, Manning WJ. Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays. Magn Reson Med. 1997;38:591–603. [PubMed]
9. Bydder M, Larkman DJ, Hajnal JV. Generalized SMASH imaging. Magn Reson Med. 2002;47:160–170. [PubMed]
10. 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. [PubMed]
11. Griswold MA, Blaimer M, Breuer F, Heidemann RM, Mueller M, Jakob PM. Parallel magnetic resonance imaging using the GRAPPA operator formalism. Magn Reson Med. 2005;54:1553–1556. [PubMed]
12. 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. [PubMed]
13. Griswold MA, Jakob PM, Nittka M, Goldfarb JW, Haase A. Partially parallel imaging with localized sensitivities (PILS) Magn Reson Med. 2000;44:602–609. [PubMed]
14. Samsonov AA, Block WF, Arunachalam A, Field AS. Advances in locally constrained k-space-based parallel MRI. Magn Reson Med. 2006;55:431–438. [PubMed]
15. Hestenes MR, Stiefel EL. Method of conjugate gradients for solving linear systems. J Res Natl Bur Standards. 1952:409–435.
16. Heberlein K, Hu X. Auto-calibrated parallel spiral imaging. Magn Reson Med. 2006;55:619–625. [PubMed]
17. Heidemann RM, Griswold MA, Seiberlich N, Kruger G, Kannengiesser SA, Kiefer B, Wiggins G, Wald LL, Jakob PM. Direct parallel image reconstructions for spiral trajectories using GRAPPA. Magn Reson Med. 2006 [PubMed]
18. Ogawa S, Menon RS, Tank DW, Kim SG, Merkle H, Ellermann JM, Ugurbil K. Functional brain mapping by blood oxygenation level-dependent contrast magnetic resonance imaging. A comparison of signal characteristics with a biophysical model. Biophys J. 1993;64:803–812. [PubMed]
19. Ogawa S, Lee TM, Kay AR, Tank DW. Brain magnetic resonance imaging with contrast dependent on blood oxygenation. Proc Natl Acad Sci USA. 1990;87:9868–9872. [PubMed]
20. Golay X, Pruessmann KP, Weiger M, Crelier GR, Folkers PJ, Kollias SS, Boesiger P. PRESTO-SENSE: an ultrafast whole-brain fMRI technique. Magn Reson Med. 2000;43:779–786. [PubMed]
21. Warach S, Levin JM, Schomer DL, Holman BL, Edelman RR. Hyperperfusion of ictal seizure focus demonstrated by MR perfusion imaging. AJNR Am J Neuroradiol. 1994;15:965–968. [PubMed]
22. Warach S, Dashe JF, Edelman RR. Clinical outcome in ischemic stroke predicted by early diffusion-weighted and perfusion magnetic resonance imaging: a preliminary analysis. J Cereb Blood Flow Metab. 1996;16:53–59. [PubMed]
23. Bammer R, Auer M, Keeling SL, Augustin M, Stables LA, Prokesch RW, Stollberger R, Moseley ME, Fazekas F. Diffusion tensor imaging using single-shot SENSE-EPI. Magn Reson Med. 2002;48:128–136. [PubMed]
24. Basser PJ, Mattiello J, LeBihan D. MR diffusion tensor spectroscopy and imaging. Biophys J. 1994;66:259–267. [PubMed]
25. Jaermann T, Crelier G, Pruessmann KP, Golay X, Netsch T, van Muiswinkel AM, Mori S, van Zijl PC, Valavanis A, Kollias S, Boesiger P. SENSE-DTI at 3 T. Magn Reson Med. 2004;51:230–236. [PubMed]
26. Liu C, Bammer R, Kim DH, Moseley ME. Self-navigated interleaved spiral (SNAILS): application to high-resolution diffusion tensor imaging. Magn Reson Med. 2004;52:1388–1396. [PubMed]
27. Liu C, Moseley ME, Bammer R. Simultaneous phase correction and SENSE reconstruction for navigated multi-shot DWI with non-Cartesian k-space sampling. Magn Reson Med. 2005;54:1412–1422. [PubMed]
28. Jackson JI, Meyer CH, Nishimura DG, AM Selection of a convolution function for Fourier inversion using gridding. IEEE Trans Med Imaging. 1991;10:473–478. [PubMed]
29. Paige CC, Saunders MA. LSQR: An algorithm for sparse linear equations and sparse least squares. ACM Trans of Math Software. 1982;8(1):43–71.
30. Benzi M, Meyer C, Tuma M. A sparse approximate inverse preconditioner for the conjugate gradient method. Siam J Sci Comput. 1996;17:1135.
31. Benson M, Frederickson P. Iterative solution of large sparse linear systems arising in certain multidimensional approximation problems. Utilitas Math. 1982;22:127–140.
32. Benson M. Master’s Thesis. Thunder Bay, Canada: Lakehead University; 1973. Iterative solution of large scale linear systems.
33. Golub GH, Van Loan CF. Matrix computations. Baltimore, MD: Johns Hopkins; 1996. p. 590.
34. Bydder M, Larkman DJ, Hajnal JV. Combination of signals from array coils using image-based estimation of coil sensitivity profiles. Magn Reson Med. 2002;47:539–548. [PubMed]
35. Griswold MA, Breuer F, Blaimer M, Kannengiesser S, Heidemann RM, Mueller M, Nittka M, Jellus V, Kiefer B, Jakob PM. Autocalibrated coil sensitivity estimation for parallel imaging. NMR Biomed. 2006;19:316–324. [PubMed]
36. Lin FH, Chen YJ, Belliveau JW, Wald LL. A wavelet-based approximation of surface coil sensitivity profiles for correction of image intensity inhomogeneity and parallel imaging reconstruction. Hum Brain Mapp. 2003;19:96–111. [PubMed]
37. Kim DH, Adalsteinsson E, Spielman DM. Simple analytic variable density spiral design. Magn Reson Med. 2003;50:214–219. [PubMed]
38. Reinsch C. Smoothing by spline functions. Numer Math. 1967;10:177–183.
39. Glover GH. Simple analytic spiral k-space algorithm. Magn Reson Med. 1999;42:412–415. [PubMed]
40. Rasche V, Proksa R, Sinkus R, Bornert P, Eggers H. Resampling of data between arbitrary grids using convolution interpolation. IEEE Trans Med Imaging. 1999;18:385–392. [PubMed]