Sensitivity Encoding Reconstruction for MRI with Gridding Algorithm

Lianjun Zhang, Gang Liu

Computer College Science and Technology, Shandong University of Technology, Zibo, China.

DOI: 10.4236/jcc.2021.92002                

Abstract

The Sensitivity Encoding (SENSE) parallel reconstruction scheme for magnetic resonance imaging (MRI) is studied and implemented with gridding algorithm in this paper. In this paper, the sensitivity map profile, field map information and the spiral k-space data collected from an array of receiver coils are used to reconstruct un-aliased images from under-sampled data. The gridding algorithm is implemented with SENSE due to its ability in evaluating forward and adjoins operators with non-Cartesian sampled data. This paper also analyzes the performance of SENSE with real data set and identifies the computational issues that need to be improved for further research.

Keywords

Parallel Imaging, SENSE, Gridding Algorithm

1. Introduction

MRI is a relatively young technology, it has reached a juncture where future advances are limited by the scanning and reconstruction time. Since current MRI scanners already operate at the physical limit of data acquisition speed, mainly due to technical difficulties in producing rapidly switching magnetic field gradients, people resort to parallel imaging techniques that apply phased array coils and parallel reconstruction methods for faster MR imaging. In the two decades that followed, the MRI community has also witnessed great progress in other parallel imaging schemes. Among these schemes perhaps the most well-known are SENSE presented by Pruessmann [1] and SMASH presented by Sodickson [2]. Compared with SMASH, SENSE has its own advantage in providing parallel imaging with arbitrary coil configuration, though at the expense of extra loss in SNR. Despite the fact that these two techniques can significantly compress scanning time, usually computationally intensive inverse procedures are required to unfold images from undersampled data, especially in the case of non-cartesian data sets. Aimed at this problem, Pruessmann proposed an interative reconstruction scheme [3] that combines fast Fourier transform (FFT) with forward and inverse gridding operation. The computation complexity of SENSE is then greatly reduced from O(N⁴) to O(N²logN). This method makes sensitivity-encoded imaging practical with arbitrary k-space acquisition pattern. In this manner, the purpose of this term project is to implement the SENSE reconstruction scheme with non-cartesian k-space trajectories, and identify computational issues that need to be improved for further generalizing this algorithm for the reconstruction of three dimensional MRI dataset. The performance of SENSE with real data set and identifies the computational issues to be improved for research.

2. SENSE Parallel Imaging Methodology in MRI for Gridding Algorithm

Parallel imaging in MRI is often performed by placing an array of receiver coils around the object to be imaged, with each independent receiver coil lending spatially distinct reception profiles to the collected data sets. The following equation indicates the formation of baseband image in a typical MRI experiment:

$\rho \left(t\right)={\displaystyle \int d\left(r\right)S\left(r\right){\text{e}}^{-i\omega \left(r\right)t}{\text{e}}^{-i2\pi k\left(t\right)r}\text{d}r}+\epsilon \left(t\right)$, (1)

where d(r) is a continuous function of the object’s traverse magnetization at location r; d(r) along with its trajectories can be obtained directly from MRI scanners in k-space, the frequency space; S(r) is the spatial sensitivity profile of the receiver coil; ω(r) is the B0 field inhomogeneity present at location r; k(t) is the input data trajectory at time t; ε(t) is the noise term at time t. Equation (1) shows that the k-space data, the coil sensitivity profile and the field inhomogeneity information are three essential components for parallel image formation. It also suggests that the data acquired from multiple independent coils and the features of those coils can be integrated into a larger system of equations. With each coil receiving its own data, weighted by S_l(r), where $l=1,\cdots ,L$, i.e. the complex spatial sensitivity profile of coil l, the parallel imaging model in an integrated matrix form can thus be represented as follows (take L = 3 as an example):

$\left[\begin{array}{c}{\rho }_1\\ {\rho }_2\\ {\rho }_3\end{array}\right]=\left[\begin{array}{c}F\cdot S_1\\ F\cdot S_2\\ F\cdot S_3\end{array}\right]d+\left[\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\\ {\epsilon }_3\end{array}\right]$, (2)

where ρ₁ denotes the signal vector received from coil 1 and the image vector ρ is formed by stacking ρ_l’s into a single column; S_l is the diagonal matrix holding the complex spatial sensitivity profiles on the diagonal entries from the lth coils; F denotes the imaging system matrix.

From the above mathematical model, parallel imaging makes use of signal processing methods to reduce the scanning time and the amount of acquired k-space data while still maintaining the same spatial resolution, i.e. the same area of k-space is still effectively covered in the Nyquist sense with the acquired data from phased array coils. However, according to common signal processing knowledge on sampling theory, if any parallel imaging scheme enables a reduced number of phase encoding lines, the interval between which is 1/Field of View (FOV), aliasing will consequently show up in the reconstructed image when conventional FFT reconstruction is employed. Figure 1 demonstrates this issue with comparison between two images reconstructed from fully sampled dataset and under-sampled dataset.

However, the SENSE parallel imaging scheme [4] combines the above mentioned aliasing effect and those aliased signals from multiple coils to reveal the gap between faster imaging speed and aliasing. Take cartesian sampled k-space as an example, the locations and distances between phase encoding lines are often known before MR experiments. From these locations and distances, it is possible to calculate the amount of aliasing accurately in the reconstructed image of each coil. Assume that the scan time is reduced by the acceleration factor of R, the FOV in the reconstructed image is reduced, or folded, by a same factor of R. Therefore, the acquired signal intensity at an aliased reconstructed pixel is indeed a superposition of signal intensity from the desired location in the object plus R-1 times of pixel values that are displaced by integer multiples of FOV/R.

The SENSE imaging scheme works in this way [5]: To obtain a desired image from the acquired full FOV image data, the image value of pixel (x, y) should be weighted with the coil sensitivity profile S at the corresponding locations (x, y_l), where $l=1,\cdots ,R$. And the accurate reconstructed pixel value at the location (x, y) with the k’th coil information is as follows:

$I_k\left(x,y\right)=S_k\left(x,y_1\right)\rho \left(x,y_1\right)+S_k\left(x,y_2\right)\rho \left(x,y_2\right)+\cdots +S_k\left(x,y_R\right)\rho \left(x,y_R\right)$. (3)

The gridding algorithm is usually employed to deal with data sets that do not fall on regular cartesian grids in the k-space. The gridding approach allows people to collect and process non-cartesian data sets in a way that previously developed cartesian reconstruction methods, such as FFT based schemes, can be applied without much modification. The gridding algorithm is performed in a way that the acquired spiral trajectory data is first re-sampled onto cartesian sampled grids before the FFTW library is called for image reconstruction. In the re-sampling process, the gridding algorithm provides an approximation of the adjoint and forward operators in the SENSE imaging scheme. The basic concept of adjoint and forward operators can be expressed as follows:

$\rho \left(x_n\right)={\displaystyle \underset{l=0}{\overset{l-1}{\sum }}{\text{e}}^{+i2\pi f\left(x_n\right)\tau l}}{\displaystyle \underset{m=0}{\overset{M}{\sum }}{a}_l\left(t\right)d\left(m\right){\text{e}}^{+i2\pi k_m{x}_n}}$, (4)

$d\left(m\right)={\displaystyle \underset{l=0}{\overset{l-1}{\sum }}{a}_l\left(t\right)}{\displaystyle \underset{n=0}{\overset{N}{\sum }}\rho \left(x_n\right){\text{e}}^{-i2\pi k_m{x}_n}{\text{e}}^{-i2\pi f\left(x_n\right)\tau l}}$, (5)

where d(m) demotes the k-space sampled data; a(t) denotes the interpolation window function for time segmentation; in this report, the hann window function is calculated and applied to time segmentation; f(x_n) denotes the B0 field inhomogeneity map; k_m denotes the k-space trajectories. From these two formulations, it is evident that the adjoint operator takes image data value and k-space trajectories as input and returns the k-space data value on the cartesian grid points as output. The forward operator, on the contrary, takes k-space data value and k-space trajectories as input and outputs the image value at each corresponding pixel locations. Although these two operators are distinct from each other, the evaluation of both operators is indispensable for the conjugate gradient method. The following paragraphs will explain the implementation of these two operators in detail.

A simple version of the adjoint operator can be concluded in three steps. In the first step, each data point in k-space is sampled with an arbitrary sampling function. The sampled data point is then convolved with a gridding kernel. in this report, the Kaiser Bessel function is used as the gridding kernel. The first step can be expressed in the following equation:

$\hat{d}\left(k_x,k_y\right)=\left[d\left(k_x,k_y\right)S\left(k_x,k_y\right)\times K\left(k_x,k_y\right)\right]$, (6)

where d(k_x, k_y) denotes the original k-space data point; S(k_x, k_y) denotes the k-space sampling function; K(k_x, k_y) denotes the Kaiser Bessel gridding kernel. In this step, the ideal image d(x, y) is actually blurred by the gridding convolution with the Fourier transform of the Kaiser Bessel gridding kernel. In the second step, the result of gridding convolution is re-sampled onto the cartesian grid point before performing the Fourier transform. It is just this step that takes the original spiral k-space data points and re-sample their contribution onto the adjacent cartesian grid points. And the corresponding formulation in the image domain is as follows:

$\hat{d}\left(x,y\right)=\left\{\left[d\left(x,y\right)\times s\left(x,y\right)\right]k\left(x,y\right)\right\}\times \psi \left(\frac{x}{FOVx},\frac{y}{FOVy}\right)$, (7)

where d(x, y) denotes the ideal image data; s(x, y) denotes the sampling function in the image domain; k(x, y) denotes the Kaiser Bessel gridding kernel in the image domain; ψ(x, y) denotes the re-sampling function in the image domain. It should be noted that the second step brings in some unwelcome effects, i.e. spiral shaped side lobes are usually created corresponding to the point spread function of the sampling pattern. To suppress the side lobes introduced in the second step, the third step deapodizes the image by the Fourier transform of the gridding kernel,

$\hat{d}\left(x,y\right)=\frac{1}{k\left(x,y\right)+const.}\left\{\left[d\left(x,y\right)\times s\left(x,y\right)\right]k\left(x,y\right)\right\}\times \psi \left(\frac{x}{FOVx},\frac{y}{FOVy}\right)$. (8)

While the third step introduces shading in image, it suppresses the side lobes created from the second step. Finally, the result of the third step is the desired image from the adjoint operator evaluated via gridding algorithm.

A simple implementation of inverse gridding is to reverse each steps in gridding. Assume that start with an image update from the conjugate gradient solver, and need to estimate the corresponding spiral data in k-space. From the above introduction of the gridding process, it is evident that the first step of inverse gridding should be apodization, i.e. to cancel the effect of the deapodization in the last step of gridding algorithm:

$d\left(x,y\right)=\frac{\hat{d}\left(x,y\right)}{k\left(x,y\right)+const.}$. (9)

The remaining steps of inverse gridding are simply the deconvolution and re-sampling in k-space. The only difference of these steps with those corresponding steps in the adjoint operator is that, this time we are estimating the values of the spiral trajectories in the k-space from the cartesian grid points, not the other way around. And the expression in k-space for the remaining steps is as follows:

$d\left(k_x,k_y\right)=\left[\frac{\hat{d}\left(x,y\right)}{k\left(x,y\right)+const.}\times K\left(k_x,k_y\right)\right]S\left(k_x,k_y\right)$. (10)

Although some literatures [6] [7] claims that density compensation is necessary either before gridding or after gridding. In actual implementation on the specific data set used in this report, however, it turns out that that these operations are redundant since there is no significant improvement in both image resolution and SNR. Thus, it is not necessary to discuss density compensation issues in detail.

3. Experimental Results and Discussion

The above figures demonstrate that the SENSE parallel imaging scheme with gridding algorithm can effectively handle non-cartesian k-space trajectories. Image resolution is increased corresponding to the data set matrix size. However, improvement in image resolution is not significant between Figure 1 and Figure 4. This is due to the fact that the original collected data set has merely 320 × 320 matrix size. Therefore, resizing matrix size by interpolation in k-space does not lead to significant improvement in image resolution. Figure 4 evaluates the computational efficiency of SENSE implementation. The estimated computation time is around 90% of actual computation timing. The difference in computation time is mainly due to memory issues.

The implementation of SENSE parallel imaging scheme has shown the effectiveness of reconstructing image. By incorporating data from multiple receiver coils, the under sampling effect in each coil’s data can be effectively compensated. The combination of conjugate gradient algorithm, forward and inverse gridding, with FFT can reduce the computation complexity as opposed to brute force methods. as can be seen from Figures 1-4, the reconstructed images contain a small portion of the artificially added Gaussian white noise in the given data set. Also, in order to optimize computation efficiency of the SENSE In future effort of this term project, this artifact can be eliminated by combining regularization methods, the Toeplitz-based iterative methods [7] [8] serves as a promising approach since it replace the repetitive evaluation of both operators with only one step forward and adjoint gridding. This approach will have great impact on the overall computation performance of the SENSE imaging scheme. These approaches will continue to improve the current SENSE implementation and lay a solid foundation for implementing three dimensional SENSE.

Figure 1. 64 × 64 reconstruction.

Figure 2. 128 × 128 reconstruction.

Figure 3. 256 × 256 reconstruction.

Figure 4. 512 × 512 reconstruction.

This approach will have great impact on the overall computation performance of the SENSE imaging scheme. These approaches will continue to improve the current SENSE implementation and lay a solid foundation for implementing three dimensional SENSE.

4. Conclusion

A wise investment of time and effort on three-dimensional SENSE and it is an exciting experience in applying critical skills to the state of the art image reconstruction techniques. The gridding algorithm is also implemented with SENSE due to its ability in evaluating forward and adjoint operator with non-cartesian sampled data. That fast template matching algorithm is capable of generating high-quality spatiotemporal atlas from sparse-sampled data in a short computation time.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Pruessmann, K. (1999) SENSE: Sensitivity Encoding for Fast MRI. Magnetic Resonance in Medicine, 42, 952-962. https://doi.org/10.1002/(SICI)1522-2594(199911)42:5<952::AID-MRM16>3.0.CO;2-S
[2]	Sodickson, D. (1997) Simultaneous Acquisition of Spatial Harmonics (SMASH): Fast Imaging with Radiofrequency Coil Arrays. Magnetic Resonance in Medicine, 38, 591-603. https://doi.org/10.1002/mrm.1910380414
[3]	Pruessmann, K. (2001) Advances in Sensitivity Encoding with Arbitary K-Space Trajectories. Magnetic Resonance in Medicine, 46, 638-651. https://doi.org/10.1002/mrm.1241
[4]	Blaimer, M. (2004) SMASH, SENSE, PILS, GRAPPA, How to Choose the Optimal Method. Topics in Magnetic Resonance Imaging, 15, 223-236.
[5]	Woo, J., Xing, F., Lee, J., Stone, M. and Prince, J.L. (2015) Construction of an Unbiased Spatio-Temporal Atlas of the Tongue during Speech. International Conference on Information Processing in Medical Imaging, 28 June-3 July, Sabhal Mor Ostaig, 723-732. https://doi.org/10.1007/978-3-319-19992-4_57
[6]	Alessandrini, M., Basarab, A., Liebgott, H. and Bernard, O. (2013) Cardiac Motion Assessment from Echocardiographic Image Sequences by Means of the Structure Multivector. 2013 IEEE International Ultrasonics Symposium (IUS), 21-25 July 2013, Prague, 1541-1544. https://doi.org/10.1109/ULTSYM.2013.0392
[7]	Woo, J., Stone, M. and Prince, J. (2015) Multimodal Registration via Mutual Information Incorporating Geometric and Spatial Context. IEEE Transactions on Image Processing, 24, 757-769. https://doi.org/10.1109/TIP.2014.2387019
[8]	Fessler, J. (2017) Toeplitz-Based Iterative Image Reconstruction for MRI with Correction for Magnetic Field Inhomogeneity. IEEE Transactions on Signal Processing, 53, 3393-3402. https://doi.org/10.1109/TSP.2005.853152