10. Non-Cartesian MR image reconstruction - MRI acquisition and Reconstruction

In this tutorial we will reconstruct an MRI image from radial undersampled kspace measurements. Let us denote $\Omega$ the undersampling mask, the under-sampled Fourier transform now reads $F_{\Omega}$ .

Import neuroimaging data¶

We use the toy datasets available in pysap, more specifically a 2D brain slice and the radial under-sampling scheme. We compare zero-order image reconstruction with Compressed sensing reconstructions (analysis vs synthesis formulation) using the FISTA algorithm for the synthesis formulation and the Condat-Vu algorithm for the analysis formulation.

We remind that the synthesis formulation reads (minimization in the sparsifying domain):

\widehat{z} = \text{arg}\,\min_{z\in C^n_\Psi} \frac{1}{2} \|y - F_\Omega \Psi^*z \|_2^2 + \lambda \|z\|_1

(1)

and the image solution is given by $\widehat{x} = \Psi^*\widehat{z}$ . For an orthonormal wavelet transform, we have $n_\Psi=n$ while for a frame we may have $n_\Psi > n$ .

while the analysis formulation consists in minimizing the following cost function (min. in the image domain):

\widehat{x} = \text{arg}\,\min_{x\in C^n} \frac{1}{2} \|y - F_\Omega x\|_2^2 + \lambda \|\Psi x\|_1 \,.

(2)

Author: Chaithya G R & Philippe Ciuciu
Date: 01/06/2021
Target: ATSI MSc students, Paris-Saclay University

from mri.operators import NonCartesianFFT, WaveletUD2, WaveletN
from mri.operators.utils import convert_locations_to_mask, \
    gridded_inverse_fourier_transform_nd
from mri.reconstructors import SingleChannelReconstructor
from pysap.data import get_sample_data

# Third party import
from modopt.math.metrics import ssim
from modopt.opt.linear import Identity
from modopt.opt.proximity import SparseThreshold
import numpy as np
import matplotlib.pyplot as plt

Loading input data¶

image = get_sample_data('2d-mri').data.astype(np.complex64)

# Obtain MRI non-cartesian mask
radial_mask = get_sample_data("mri-radial-samples")
kspace_loc = radial_mask.data

View Input data¶

plt.subplot(1, 2, 1)
plt.imshow(np.abs(image), cmap='gray')
plt.title("MRI Data")
plt.subplot(1, 2, 2)
plt.imshow(convert_locations_to_mask(kspace_loc, image.shape), cmap='gray')
plt.title("K-space Sampling Mask")
plt.show()

Generate the kspace¶

From the 2D brain slice and the acquisition mask, we retrospectively undersample the k-space using a radial mask. We then reconstruct the zero-order solution as a baseline

Get the locations of the kspace samples

# Get the locations of the kspace samples and the associated observations
fourier_op = NonCartesianFFT(samples=kspace_loc, shape=image.shape, implementation='finufft')
kspace_obs = fourier_op.op(image)

Gridded solution

grid_space = np.linspace(-0.5, 0.5, num=image.shape[0])
grid2D = np.meshgrid(grid_space, grid_space)
grid_soln = gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs,
                                                 tuple(grid2D), 'linear')
base_ssim = ssim(grid_soln, image)
plt.imshow(np.abs(grid_soln), cmap='gray')
plt.title('Gridded solution : SSIM = ' + str(np.around(base_ssim, 3)))
plt.show()

FISTA optimization¶

We now want to refine the zero order solution using a FISTA optimization. The cost function is set to Proximity Cost + Gradient Cost

linear_op = WaveletN(wavelet_name="sym8", nb_scales=4)
regularizer_op = SparseThreshold(Identity(), 6 * 1e-7, thresh_type="soft")

Generate operators¶

reconstructor = SingleChannelReconstructor(
    fourier_op=fourier_op,
    linear_op=linear_op,
    regularizer_op=regularizer_op,
    gradient_formulation='synthesis',
    verbose=1,
)

Synthesis formulation: FISTA optimization¶

We now want to refine the zero order solution using a FISTA optimization. The cost function is set to Proximity Cost + Gradient Cost

image_rec, costs, metrics = reconstructor.reconstruct(
    kspace_data=kspace_obs,
    optimization_alg='fista',
    num_iterations=200,
)

recon_ssim = ssim(image_rec, image)
plt.imshow(np.abs(image_rec), cmap='gray')
plt.title('FISTA Reconstruction : SSIM = ' + str(np.around(recon_ssim, 3)))
plt.show()

POGM reconstruction¶

image_rec2, costs2, metrics2 = reconstructor.reconstruct(
    kspace_data=kspace_obs,
    optimization_alg='pogm',
    num_iterations=200,
)

recon2_ssim = ssim(image_rec2, image)
plt.imshow(np.abs(image_rec2), cmap='gray')
plt.title('POGM Reconstruction : SSIM = ' + str(np.around(recon2_ssim, 3)))
plt.show()

Analysis formulation: Condat-Vu reconstruction¶

#linear_op = WaveletN(wavelet_name="sym8", nb_scales=4)
linear_op = WaveletUD2(
    wavelet_id=24,
    nb_scale=4,
)

reconstructor = SingleChannelReconstructor(
    fourier_op=fourier_op,
    linear_op=linear_op,
    regularizer_op=regularizer_op,
    gradient_formulation='analysis',
    verbose=1,
)

image_rec, costs, metrics = reconstructor.reconstruct(
    kspace_data=kspace_obs,
    optimization_alg='condatvu',
    num_iterations=100,
)
plt.imshow(np.abs(image_rec), cmap='gray')
recon_ssim = ssim(image_rec, image)
plt.title('Condat-Vu Reconstruction\nSSIM = ' + str(recon_ssim))
plt.show()