Research of Sensitivity Encoding Reconstruction for MRI with Non-Cartesian K-Space Trajectories

Lianjun Zhang¹, Lifang Zhao², Gang Liu¹

¹School of Computer, Shandong University of Technology, Zibo, China.
²Library, Shandong University of Technology, Zibo, China.

DOI: 10.4236/jcc.2021.91001                

Abstract

The Sensitivity Encoding (SENSE) parallel reconstruction scheme for magnetic resonance imaging (MRI) is implemented with non-cartesian sampled k-space trajectories in this paper. SENSE has the special capability to reduce the scanning time for MRI experiments while maintaining the image resolution with under-sampling data sets. In this manner, it has become an increasingly popular technique for multiple MRI data acquisition and image reconstruction schemes. The gridding algorithm is also implemented with SENSE due to its ability in evaluating forward and adjoin operator with non-cartesian sampled data. 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 unaliased images from under-sampled data. The performance of SENSE with real data set identifies the computational issues to be improved for researched.

Keywords

Parallel Imaging, SENSE, Gridding, K-Space, MRI

1. Introduction

Although 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. Phased array coils were first exploited in 1987 to reduce phase-encoding lines [1], which is proportional to the amount of scanning time. 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 [2] and SMASH presented by Sodickson [3]. 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 under sampled data, especially in the case of non-cartesian data sets. Aimed at this problem, Pruessmann proposed an interative reconstruction scheme [4] that combines fast Fourier transform (FFT) with forward and inverse gridding operation. The computation complexity of SENSE is then greatly reduced from O(N4) to O(N2logN). 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, evaluate its performance with “Double Vision” data set from the “ISMRM Reconstruction Challenge”, and identify computational issues that need to be improved for further generalizing this algorithm for the reconstruction of three dimensional MRI dataset.

2. Sense Parallel Imaging in MRI

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.

However, the SENSE parallel imaging scheme [5] 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: 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)

This equation can be transformed into matrix notation and illustrated in Figure 1:

Figure 1. Basic SENSE imaging scheme relation.

$I=S\cdot \rho$ (4)

To simplify the formulation, in Equation (4) we do not consider the issue of noise correlation.

According to the matrix theory knowledge, we can easily obtain the full FOV image from taking generalized of the sensitivity matrix S:

$\rho ={\left(S^{H}S\right)}^{-1}{S}^H\cdot I$ (5)

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. In this report, 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}}$ , (6)

$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}}$ , (7)

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.

Finally, the result of the third step is the desired image from the adjoint operator evaluated via gridding algorithm. 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

The reconstructed data set in this report is the “Double Vision” data of the “ISMRM Reconstruction Challenge”. The “Double Vision” data originate from 12 axial images in the abdomen collected using a torso phase-array coil, collected at 320 × 320 matrix and resized to four different sizes, 64 × 64 matrix, 128 × 128 matrix, 256 × 256 matrix and 512 × 512 matrix, to test the effectiveness of reconstruction. Field maps were collected using gradient echo images at TE = 3, with a 96 × 96 collected matrix resized to the above mentioned four different sizes. Sensitivity maps are also provided with the dataset. The original data sets were synthesized over 8 spiral trajectories, each with 20,000 points. The k-space trajectories at the first and last interleaves are shown in Figure 2 and Figure 4. The gray region in both figures represents areas that past interleaves have covered in the sequence. The reconstructed images of four different sizes are shown by Figures 3-7.

Figure 2. First shot trajectory.

Figure 3. Last shot trajectory.

Figure 4. 64 × 64 reconstruction.

Figure 5. 128 × 128 reconstruction.

Figure 6. 256 × 256 reconstruction.

Figure 7. 512 × 512 reconstruction.

The above figures demonstrate that the SENSE parallel imaging scheme with gridding algorithm can effectively handle non-cartesian k-space trajectories. Compare Figure 5 with Figure 6 and Figure 7, image resolution increase corresponding to the data set matrix size. However, improvement in image resolution is not significant between Figure 6 and Figure 7. 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. Table 1 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 from arbitrary k-space trajectories. 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, together with FFT can reduce the computation complexity as opposed to brute force methods. Still, as can be seen from Figure 6 to Figure 7, the reconstructed images contain a small portion of the artificially added Gaussian white noise in the given data set. In future effort of this term project, this artifact can be eliminated by combining regularization methods with SENSE implementation [8]. Also, in order to optimize computation efficiency of the SENSE implementation [9], the Toeplitz-based iterative methods [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 reconstruction.

Table 1. Computation efficiency of SENSE implementation.

4. Conclusion

Sense parallel imaging based on gridding algorithm can improve the efficiency and speed of image reconstruction from arbitrary k-space trajectory. This term project is a wise investment of time and effort for my future research effort on three dimensional SENSE and it is exciting experience in applying critical thinking skills to the state of the art image reconstruction techniques in my research area. The reconstructed image contains a small number of Additive white Gaussian noise in a given dataset, which can be eliminated by a combination of regularization and SENSE implementations.

Conflicts of Interest

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

References

[1]	Carlson, J. (1987) An Algorithm for NMR Imaging Reconstruction Based on Multiple RF Receiver Coils. Journal of Magnetic Resonance, 74, 376-380. https://doi.org/10.1016/0022-2364(87)90348-9
[2]	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%3C952::AID-MRM16%3E3.0.CO;2-S
[3]	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
[4]	Pruessmann, K. (2001) Advances in Sensitivity Encoding with Arbitrary k-Space Trajectories. Magnetic Resonance in Medicine, 46, 638-651. https://doi.org/10.1002/mrm.1241
[5]	Blaimer, M. (2004) SMASH, SENSE, PILS, GRAPPA, How to Choose the Optimal Method. Topics in Magnetic Resonance Imaging, 15, 223-236. https://doi.org/10.1097/01.rmr.0000136558.09801.dd
[6]	Alessandrini, M., Basarab, A., Liebgott, H. and Bernard, O. (2013) Cardiac Motion Assessment from Echocardiographic Image Sequences by Means of the Structure Multivector. In: IEEE International Ultrasonics Symposium (IUS), Prague, July 2013, 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]	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. Information Processing in Medical Imaging, Springer, 723-732. https://doi.org/10.1007/978-3-319-19992-4_57
[9]	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