12. Self-calibrated CS-pMR image reconstruction from undersampled non-Cartesian data

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}$ .

We use the toy datasets available in pysap, more specifically a 2D brain slice and under-sampled Cartesian acquisition over 32 channels. We compare zero-order image reconstruction with self-calibrated multi-coil Compressed sensing reconstructions (analysis vs synthesis formulation) using the FISTA algorithm for the synthesis formulation and the Condat-Vu algorithm for the analysis formulation. The multicoil data $(y_\ell)_\ell$ is collected across multiple, say $L$ , channels. The sensitivity maps $(S_\ell)_\ell$ are automically calibrated from the central portion of k-space (e.g. 5%) for all channels $\ell=1, \ldots, L$ .

We remind that the synthesis formulation of the non-Cartesian CS-PMRI problem reads (minimization in the sparsifying domain):

\widehat{z} = \text{arg}\,\min_{z\in C^n_\Psi} \frac{1}{2} \sum_{\ell=1}^L\|y_\ell - F_\Omega S_\ell \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} \sum_{\ell=1}^L \|y_\ell - F_\Omega S_\ell 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

# Package import
from mri.operators import NonCartesianFFT, WaveletN, WaveletUD2
from mri.operators.utils import convert_locations_to_mask, \
    gridded_inverse_fourier_transform_nd
from mri.reconstructors import SelfCalibrationReconstructor
from mri.reconstructors.utils.extract_sensitivity_maps import get_Smaps

import pysap
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¶

# Loading input data
cartesian_ref_image = get_sample_data('2d-pmri')
#image = pysap.Image(data=np.sqrt(np.sum(cartesian_ref_image.data**2, axis=0)))
image = np.linalg.norm(cartesian_ref_image, axis=0)

# Obtain MRI non-Cartesian radial mask
radial_mask = get_sample_data("mri-radial-samples")
kspace_loc = radial_mask.data
mask = pysap.Image(data=convert_locations_to_mask(kspace_loc, image.shape))

plt.subplot(1, 2, 1)
plt.imshow(np.abs(image), cmap='gray')
plt.title("MRI Data")
plt.subplot(1, 2, 2)
plt.imshow(mask, 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 non cartesian acquisition mask. We then grid the kspace to get the gridded solution as a baseline.

# Get the kspace observation values for the kspace locations
fourier_op = NonCartesianFFT(
    samples=kspace_loc,
    shape=image.shape,
    n_coils=cartesian_ref_image.shape[0],
    implementation='finufft' # Use finufft for speed
)

kspace_obs = fourier_op.op(cartesian_ref_image.data)

Gridded solution

# Gridded solution
grid_space = np.linspace(-0.5, 0.5, num=image.shape[0])
grid2D = np.meshgrid(grid_space, grid_space)
grid_soln = np.asarray([
    gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs_ch,
                                         tuple(grid2D), 'linear')
    for kspace_obs_ch in kspace_obs
])
image_rec0 = pysap.Image(data=np.sqrt(np.sum(np.abs(grid_soln)**2, axis=0)))

plt.imshow(image_rec0, cmap='gray')
base_ssim = ssim(image_rec0, image)
print('The Base SSIM is : {}'.format(base_ssim))

## Estimate Sensitivity maps (Smaps)

# Obtain SMaps
Smaps, SOS = get_Smaps(
    k_space=kspace_obs,
    img_shape=fourier_op.shape,
    samples=kspace_loc,
    thresh=(0.01, 0.01),    # The cutoff threshold in each kspace direction
                            # between 0 and kspace_max (0.5)
    min_samples=kspace_loc.min(axis=0),
    max_samples=kspace_loc.max(axis=0),
    mode='gridding',
    method='linear',
    n_cpu=-1,
)

h=3;w=5;
f, axs = plt.subplots(h,w)
for i in range(h):
    for j in range(w):
        axs[i, j].imshow(np.abs(Smaps[3 * i + j]))
        axs[i, j].axis('off')
plt.subplots_adjust(wspace=0,hspace=0)
plt.show()

!pip install pynfft2

Setup Fourier operators with SENSE¶

# Setup Fourier Operator with SENSE. This would initialize the
# fourier operators in the GPU.
# and also pass the Smaps calculated above

fourier_op_sense = NonCartesianFFT(
    samples=kspace_loc,
    shape=image.shape,
    n_coils=cartesian_ref_image.shape[0],
    smaps=Smaps.astype(np.complex64),
    implementation="finufft",
)

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

# Setup the operators
linear_op = WaveletN(wavelet_name='sym8', nb_scale=4)
regularizer_op = SparseThreshold(Identity(), 1e-8, thresh_type="soft")

# Setup Reconstructor
reconstructor = SelfCalibrationReconstructor(
    fourier_op=fourier_op_sense,
    linear_op=linear_op,
    regularizer_op=regularizer_op,
    gradient_formulation='synthesis',
    lipschitz_cst=fourier_op_sense.impl.get_lipschitz_cst(),
    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

x_final, costs, metrics = reconstructor.reconstruct(
    kspace_data=kspace_obs.astype(np.complex64),
    optimization_alg='fista',
    num_iterations=100,
)

image_rec = pysap.Image(data=x_final)
recon_ssim = ssim(image_rec, image)

plt.imshow(np.abs(image_rec), cmap='gray')
plt.title('FISTA Reconstruction\nSSIM = {}'.format(recon_ssim))
plt.show()

POGM reconstruction¶

x_final, costs, metrics = reconstructor.reconstruct(
    kspace_data=kspace_obs,
    optimization_alg='pogm',
    num_iterations=200,
)

image_rec = pysap.Image(data=np.abs(x_final))
recon_ssim = ssim(image_rec, image)

plt.imshow(np.abs(image_rec), cmap='gray')
recon_ssim = ssim(image_rec, image)
plt.title('POGM Reconstruction\nSSIM = {}'.format(recon_ssim))
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,
)
regularizer_op = SparseThreshold(Identity(), 1e-9, thresh_type="soft")

# Setup Reconstructor
reconstructor = SelfCalibrationReconstructor(
    fourier_op=fourier_op_sense,
    linear_op=linear_op,
    regularizer_op=regularizer_op,
    gradient_formulation='analysis',
    verbose=1,
)

x_final, costs, metrics = reconstructor.reconstruct(
    kspace_data=kspace_obs,
    optimization_alg='condatvu',
    num_iterations=200,
)
image_rec = pysap.Image(data=np.abs(x_final))
plt.imshow(np.abs(image_rec), cmap='gray')
recon_ssim = ssim(image_rec, image)
plt.title('Condat-Vu Reconstruction\nSSIM = {}'.format(recon_ssim))
plt.show()