Progressive Magnetic Resonance Image Reconstruction Based on Iterative Solution of a Sparse Linear System

Image reconstruction from nonuniformly sampled spatial frequency domain data is an important problem that arises in computed imaging. Current reconstruction techniques suffer from limitations in their model and implementation. In this paper, we present a new reconstruction method that is based on solving a system of linear equations using an efficient iterative approach. Image pixel intensities are related to the measured frequency domain data through a set of linear equations. Although the system matrix is too dense and large to solve by direct inversion in practice, a simple orthogonal transformation to the rows of this matrix is applied to convert the matrix into a sparse one up to a certain chosen level of energy preservation. The transformed system is subsequently solved using the conjugate gradient method. This method is applied to reconstruct images of a numerical phantom as well as magnetic resonance images from experimental spiral imaging data. The results support the theory and demonstrate that the computational load of this method is similar to that of standard gridding, illustrating its practical utility.


INTRODUCTION
In this paper, we consider the problem of image reconstruction in magnetic resonance imaging (MRI) without the loss of generality to other tomographic modalities. In the past decade, there has been a great deal of interest and advances in ultrafast MRI. Echo-planar imaging (EPI) is a technique that is widely used for dynamic studies. However, it is limited in speed by how fast the gradients can be switched [1]. To lessen this limitation, nonrectilinear k-space trajectories such as spiral trajectory have been used to allow smoother gradient waveforms. This led to a more efficient acquisition process besides the reduction of inhomogeneity and motion artifacts, but image reconstruction from the nonuniformly sampled k-space data remains less than optimal in speed and accuracy. Given that conventional MRI relies on discrete Fourier transformation as its reconstruction tool, most existing methods reconstruct the image by gridding, that is, estimating a uniformly sampled k-space from the nonuniformly sampled version.
In his seminal work published in 1985, O'Sullivan [2] showed that the ideal resampling kernel is theoretically an infinite Sinc function and that the convolution with such a kernel followed by sampling on a rectilinear grid results in accurate reconstruction. Furthermore, the optimization of an objective function based on a truncated kernel determined that the prolate-spheroidal function is the optimal kernel, which can be well approximated by the Kaiser-Bessel window. This procedure has been utilized in practice as the classical gridding technique for the past two decades. Several papers addressed various aspects of this method, including the optimal selection of the width and form of the gridding kernel for spiral imaging [3], the practical implementation and applications to clinical studies [4,5], the compensation of sampling density variations [6], the efficient use of computer arithmetic for faster computation [7], and the optimized selection of sampling density compensation factors [8]. The performance of nonrectilinear sampling in the presence of inhomogeneities was also shown to be more robust [9]. Therefore, it is fair to say that gridding is well established in MRI reconstruction. Nevertheless, several issues remain as open questions in using such a method. These issues include the choice of kernel size and shape for 2 International Journal of Biomedical Imaging different k-space trajectories, the possible consequences of undersampling some areas in the k-space, and the reduction of the characteristic artifacts associated with gridding.
The above issues led to the development of methods for estimating the rectilinear points from the nonrectilinear ones as a solution to a linear system of equations. This system of equations has a rather large size that is equal to the total number of points. That is, to obtain N × N rectilinear k-space samples, a linear system of size N 2 must be solved. Even though the solution to this system would in fact be the optimal solution to the gridding problem, it is not feasible in practice. As a result, a number of attempts have been directed towards obtaining practical approximate solutions to this system of equations [10][11][12][13]. Moreover, other authors attempted to estimate the image from the sampled k-space using similar techniques [14][15][16][17][18]. Even though these methods offer some advantages over the conventional gridding method, they still have some limitations. First, there is no direct measure of the solution within their formulation and the methods are inflexible in trading off computation time and reconstruction accuracy. The latter is particularly important given the general consensus that gridding techniques are perceived as slow and inflexible in reconstruction time. Second, the computational requirements of the image domain solution techniques [14][15][16][17][18] are still very high for routine use. Therefore, another solution strategy that offers to overcome these limitations is still desired.
In this work, we describe a novel solution to the problem of reconstruction from nonuniformly sampled spatial frequency data. The problem is formulated as a linear system of equations that maps the acquired k-space points into the pixel values that represent the image. Even though the linear system is large in size and dense, a simple transformation and thresholding are utilized to convert the system into a sparse form that requires much smaller computational and storage efforts. The sparse system is then solved using the conjugate gradient iterative technique, which enables further control of accuracy versus the computation time. The theory of the new method is described and the results of applying the technique to reconstruct images of a numerical phantom as well as data from a spiral imaging sequence are presented.

THEORY
Consider f (x, y) as a continuous-space spatial domain intensity distribution and let the available k-space (frequency domain) samples be F(kx i , ky i ), where i = 0, 1, . . . , L − 1, and L is the number of samples. The spatial domain can be modeled as piecewise constant function consisting of the sum of shifted gate-like functions representing the pixels of the image at the desired resolution. That is, Consequently, the continuous Fourier transform of f (x, y) can be obtained as which reduces to Here w x and w y correspond to the half of the pixel width in x and y directions, respectively. Given the available arbitrarily located samples in the k-space, the above equation leads to a linear system whose equations are obtained by substituting the set of L values for k x and k y . Hence, the linear system matrix has dimensions of L rows and M ×N columns. The righthand side vector of this linear system consists of the values of the measured data indexed by i. Given this index and from the known locations of the samples in the k-space, the system matrix can be computed for a given arrangement of the pixels. The unknowns in this linear system are the pixel intensity values. Specifically, the matrix equation takes the form in (4) and using normalized uniformly sampled spatial points (i.e., x n = n · Δx and y m = m · Δy), (4) takes the form in (5) Yasser M. Kadah et al.
In order to simplify the solution, we multiply the rows of the system matrix by the N · M-point discrete Fourier transform matrix H in the following form: Note that the inverse of this matrix is simply its Hermitian. This multiplication compacts most of the energy within each row into only a small number of elements in a similar fashion to the transform coding-based compression techniques. As a result, instead of the original dense format of the system matrix, we obtain an equivalent system with a sparse matrix, allowing us to tap into the rich literature of matrix storage and solution techniques in this area. In order to keep the linear system unchanged while taking advantage of this property, we multiply by H and its inverse such that where V is the 1D M · N-point discrete Fourier transform of the 1D listing of all pixel values in the 2D image and the matrix M has most of the energy in each row compacted into only a few elements. In order to make the matrix M really sparse, in each row, the element values are sorted by absolute value and only the first few row elements above a threshold are retained. Once this process is completed, the matrix is stored in a sparse matrix format and the conjugate gradient iteration is applied to solve the linear system. A block diagram of the new method is illustrated in Figure 1.

Energy compacting transformation
The energy compacting transformation is implemented by Fourier-transforming the rows of the linear system matrix A. The outcome of this transformation is sorted by magnitude and only the largest few elements are retained such that their energy is above a certain fraction of the total energy in each row. Subsequently, such elements are stored by their row and column positions and values as a part of the new sparse linear system matrix. Given that this process is done in a row-by-row fashion, the storage demands for this method are not very high since the huge system matrix need not be constructed or stored.

Sparse matrix manipulation
The storage of the above-derived sparse system can be performed using several techniques. The row-indexed storage method [19], which requires only twice the size of the nonzero elements for their storage, is used in our implementation. With this representation, matrix-vector multiplication operations are only equal in complexity to the number of nonzero elements and always computationally feasible. Also, Hermitian operations are rather simple and do not pose a substantial computational burden.

Conjugate gradient iterative solution
The method of conjugate gradient is used to solve the sparse system. It optimizes the solution of a linear system by removing the error components in a number of directions that span the space of the solution [20]. The original formulation of this iteration requires the system to be real, square, symmetric, and positive definite for the algorithm to work and provide a unique solution to the system [21]. Here, a variant of this method is applied to solve the normal equations of the system where the system matrix (termed the Grammian matrix) is complex Hermitian, positive definite, and square. The approach based on the delta-criterion [22] has a computational advantage over other methods [9] in the way the stopping condition is checked every step where no further matrix-vector multiplication is involved. In particular, the conjugate gradient algorithm for solving the normal equation A H A x = A H b is described in Algorithm 1.

Simulations and experimental data acquisition
The new method was implemented to construct images from numerical phantoms as well as from experimental magnetic resonance imaging data acquired using a spiral imaging sequence. The numerical phantom data were constructed from polar samples of the analytical expression of the Shepp-Logan phantom [14]. The experimental data were acquired on a Siemens Magnetom Trio 3T MRI scanner using a spiral sequence. The spiral sequence had 16 interleaves with 2596 points each to sample a 256 × 256 k-space matrix using an eight-channel phased array heal coil. The sequence parameters were spin echo, TR / TE = 1000/15 ms, and FOV: 25.6 cm × 25.6 cm. The spiral gradient waveform was designed with a maximum gradient of 40 mT/m and a minimum rise time of 0.5 milliseconds. The spiral gradient readout had a duration of 13 milliseconds per segment and was sampled at 200 kHz. Figure 2 illustrates the energy compacting transformation based on simulation. A system matrix was constructed for a 64 × 64 gridding problem under polar sampling as shown in Figure 2(a). The resultant system matrix, which has a size of 64 2 × 64 2 , is shown in Figure 2(b). As can be seen, very few elements are nonzero as opposed to the original dense matrix. The average number of elements per row needed to retain 92% of the kernel energy for this particular case was 3, a significant reduction of elements. The result of reconstructing the image of the experimental data using conventional gridding is shown in Figure 3(a). This was the best result obtained from this method, using a Kaiser-Bessel window of width 3 and parameter 9. The results obtained using the new method are illustrated in Figures 3(b)-3(d) for energy levels of 90%, 80%, and 70%, respectively. The average numbers of sparse matrix elements/row for these three cases were 3.9, 1.49, and 1.14 elements, respectively. The number of iterations used for all images was 6. As can be observed, the quality of all of them is comparable to that of gridding while having a computational  advantage especially for the last image (nearly half the reconstruction complexity of conventional gridding). In Figures 3(e)-3(h), specially-windowed (window center set to 10% of the average intensity) versions of all images are shown to illustrate the residual artifacts. As can be seen, the artifacts differ between conventional gridding and the new method where the latter exhibits vertical and horizontal aliases of the same shape as the object. Note also that the level of residual artifacts increases as the energy retained decreases as expected. The selection of the energy level to be retained controls the average number of elements per row in the sparse matrix conversion process and hence the reconstruction accuracy and complexity. Figures 6 and 4 illustrate the effect of energy retained on the average number of elements per row in the sparse matrix and also the relative error in the reconstructed image. These figures suggest an optimal range between 85% and 95% of the energy to be retained for substantial reduction in complexity while maintaining a sufficiently low error. It is worth noting that the performance of conventional gridding was 8.5 for the equivalent average number of elements per row with an error of 0.37%. Even though the conventional gridding with parameters shown above needs only a single iteration as opposed to the new method, the complexity of multiple iterations of the new method can still be lower than that. Given the significantly higher error, it is evident that the new method has an advantage over the conventional gridding.

RESULTS AND DISCUSSION
To illustrate the rapid convergence of the conjugate gradient iteration, the results of reconstructing the image above with 90% of the energy retained for different numbers of iterations are shown in Figure 5. As can be seen, diagnostic quality reconstruction can be achieved with as low as 2-3 iterations. The convergence of this experiment is shown in Figure 7 where the delta-criterion was computed and used as the error measure. The error converges very fast to reach a level of approximately 10 −6 in the 10th iteration. Note that the oscillatory form of the convergence is mainly due to the fact that the problem solved is an approximation to the original one.
It should be noted that the size and location of the large matrix elements that are selected to achieve the kernel energy percentage varies from one row to another. This means that the kernel used to perform the mapping is spatially varying. We note also that the truncation using this method is optimal in the least-squares sense. The selection using this method avoids having the user choose more specific parameters of reconstruction like window size or neighborhood definition like other techniques [3,10]. Instead, a single parameter is used as a quality measure (the energy percentage) while another is used as a computation time control (number of conjugate gradient iterations).
The computational complexity of the new method is O(N 2 ) to obtain the solution to the linear system in addition to one N 2 -point 1D FFT required afterwards. This computational complexity is well within the range of that for conventional gridding methods as well as gridding techniques based on localized SVD solution. This complexity is several orders of magnitude below other iterative reconstruction techniques [14][15][16][17][18] given its linear dependence on the number of points in the image in the matrix inversion part. The ability of the user to control the reconstruction time is a unique feature in this method. This allows a quick, almost real-time computation of images for fast viewing by the radiologist. Furthermore, it allows the radiologist to increase the accuracy of the reconstruction of a selected image at will simply by allowing additional iterations to run. Given the low complexity of iterations, such process can be performed during image viewing using console control just like zooming or gamma curve selection with virtually no noticeable delay. The method of conjugate gradient has a number of advantages. First only a few iterations are usually required to reach a good accuracy for the solution. Another advantage is that the solution accuracy can be traded off with computation time rather flexibly. This allows the method to be customized for the particular application at hand by selecting a predetermined number of iterations that correspond to the desired computation time. If a high accuracy is desired, an efficient implementation of this method may rely on a measure of the solution update in such a way that the stopping criterion is to have an update that is insignificant compared to the present solution.
The potential applications of the new method are several. The most obvious is its use to generate real-time reconstruction of spiral imaging sequences. This can be useful in applications such as functional MRI (fMRI) where real-time activation detection and display are increasingly more available on commercial scanners. Another application is to generate fast spiral navigator reconstructions for use with segmented acquisition sequences. The new method can potentially be modified to include special reconstruction procedures such as field inhomogeneity correction where a distortion in the spatial domain needs to be taken into account. In this case, the reconstruction and artifact suppression can be achieved in a single step. Applications also extend to other modalities such as computerized tomography (CT). Further investigation of such applications is needed to assess their value in practice.

CONCLUSIONS
A new method for image reconstruction from nonuniformly sampled spatial frequency domain data is presented. This method provides the means for solving the reconstruction