A note on the iterative MRI reconstruction from nonuniform k-space data

In magnetic resonance imaging (MRI) the raw data is measured in k-space, the domain of spatial frequencies. Methods that use a non-Cartesian (e.g. spiral) sampling grid in k-space are becoming increasingly important. Reconstruction is usually performed by resampling the data onto a Cartesian grid and use the standard Fast Fourier Transformation (FFT). This model is called gridding. Another approach, the inverse model, is based on an implicit discretisation. The gridding model is the explicit computation of the picture with given Fourier samples. The inverse model is the implicit computation.
We solve both by using the NFFT 3.0 library of S. Kunis and D. Potts. The NFFT 3.0 which includes the NFFT as well as the inverse NFFT can simply be used for both approaches.
For more information please read related paper.
From the MR Research Center at the University of Aarhus (Denmark) and from Philips Research (Hamburg, Germany) we got real-life data from a MRT scanner.

The above picture shows the reconstructed data of Philips Research. Please click on the images to see all slices.

For numerical results we used a 3-dimensional shepp-logan phantom as example. The left picture shows the original phantom of size 256*256*36. The right picture shows the image we reconstructed with SPIRAL data and Voronoi weights after ten iteration. Please click on the images to see all 36 slices.

18th slice of the phantom (left) and the profile of the 128th row of this slice (right) after 1-10 iterations with INTERLEAVED SPIRAL data and without weights. Please click on the images to see the sequence.

64th slice of the 128x128x128 phantom (left) and the profile of the 64th row of this slice (right) after 1-30 iterations with 3d-RADIAL data and without weights. Please click on the images to see the sequence. The RMS after 30 iteration is 0.3309. Note that the RMS after 30 iteration with weights is 0.2104. The result is not improved within further iterations (see Table 4.1 of the related paper).

18th slice of the phantom (top) and the profile of the 128th row of this slice (bottom) after one (left) , five (middle) and ten iterations (right) with SPIRAL data, 30% noise and Voronoi weights.

64th slice of the phantom (top) and the profile of the 64th row of this slice (bottom) with 3d-RADIAL data and without weights : from left to right after one, five and ten iterations.

64th slice of the phantom (top) and the profile of the 64th row of this slice (bottom) with 3d-RADIAL data and analytic weights : from left to right after one, five and ten iterations.

18th slice of the phantom (top) and the profile of the 128th row of this slice (bottom) after one iteration with SPIRAL data: from left to right, without weights, with approximative weights and with Voronoi weights.

18th slice of the phantom (top) and the profile of the 128th row of this slice (bottom) after ten iteration with SPIRAL data: from left to right, without weights, with approximative weights and with Voronoi weights.

18th slice of the phantom (top) and the profile of the 128th row of this slice (bottom) after one iteration with RADIAL data: from left to right, without weights, with approximative weights and with Voronoi weights.

18th slice of the phantom (top) and the profile of the 128th row of this slice (bottom) after ten iteration with RADIAL data: from left to right, without weights, with approximative weights and with Voronoi weights.

18th slice of the phantom (top) and the profile of the 128th row of this slice (bottom) after one iteration with INTERLEAVED SPIRAL data: from left to right, without weights, with approximative weights and with Voronoi weights.

Our programs are written in C and MATLAB. Please download the NFFT3 and see the folder ~/examples/lowlevel/mri . First one has to install FFTW and NFFT 3.0. Read the containing README file. Different nodes, different sizes of the phantom and different reconstruction methods can be chosen. The software can also handle other phantoms or real-life MRT data. An example for the usage is following:

% N is the number of points per row /column, Z is the number of slices.
% First we construct the Phantom. Then we construct the nodes and compute
% the k-Space data of the phantom. Next we compute the weights. The reconstruction follows.
% Then we visualize the phantom and compute the signal to noise ratio.

N = 256; Z = 36;

construct_phantom(N,Z);

M = construct_nodes_radial(N);

system(['./construct_data ' 'output_phantom_nfft.dat ' int2str(N) ' ' int2str(M) ' ' int2str(Z)]);

precompute_weights('output_phantom_nfft.dat',M);

system(['./reconstruct_data_2+1d ' 'output_data_phantom_nfft.dat ' ...
         int2str(N) ' ' int2str(M) ' ' int2str(Z)  ' 1 1']);

visualize_phantom('pics_2+1d/pic', N,Z);

snr('pics_2+1d/snr.txt');