^{1}

^{*}

^{1}

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.

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 [^{4}) to O(N^{2}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.

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:

ρ ( t ) = ∫ d ( r ) S ( r ) e − i ω ( r ) t e − i 2 π k ( t ) r d r + ε ( t ) , (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 , ⋯ , 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):

[ ρ 1 ρ 2 ρ 3 ] = [ F ⋅ S 1 F ⋅ S 2 F ⋅ S 3 ] d + [ ε 1 ε 2 ε 3 ] , (2)

where ρ_{1} 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 [

The SENSE imaging scheme works in this way [_{l}), where l = 1 , ⋯ , R . And the accurate reconstructed pixel value at the location (x, y) with the k’th coil information is as follows:

I k ( x , y ) = S k ( x , y 1 ) ρ ( x , y 1 ) + S k ( x , y 2 ) ρ ( x , y 2 ) + ⋯ + S k ( x , y R ) ρ ( x , y R ) . (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:

ρ ( x n ) = ∑ l = 0 l − 1 e + i 2 π f ( x n ) τ l ∑ m = 0 M a l ( t ) d ( m ) e + i 2 π k m x n , (4)

d ( m ) = ∑ l = 0 l − 1 a l ( t ) ∑ n = 0 N ρ ( x n ) e − i 2 π k m x n e − i 2 π f ( x n ) τ 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:

d ^ ( k x , k y ) = [ d ( k x , k y ) S ( k x , k y ) × K ( k x , k y ) ] , (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:

d ^ ( x , y ) = { [ d ( x , y ) × s ( x , y ) ] k ( x , y ) } × ψ ( x F O V x , y F O V y ) , (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,

d ^ ( x , y ) = 1 k ( x , y ) + c o n s t . { [ d ( x , y ) × s ( x , y ) ] k ( x , y ) } × ψ ( x F O V x , y F O V y ) . (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 ( x , y ) = d ^ ( x , y ) k ( x , y ) + c o n s t . . (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 ( k x , k y ) = [ d ^ ( x , y ) k ( x , y ) + c o n s t . × K ( k x , k y ) ] S ( k x , k y ) . (10)

Although some literatures [

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

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 [

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.

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.

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

Zhang, L.J. and Liu, G. (2021) Sensitivity Encoding Reconstruction for MRI with Gridding Algorithm. Journal of Computer and Communications, 9, 22-28. https://doi.org/10.4236/jcc.2021.92002