Cartesian Calibrationless Reconstruction using GroupLASSO Regularizer

Author: Chaithya G R

In this tutorial we will reconstruct an MRI image from cartesian kspace measurements.

Import neuroimaging data

We use the toy datasets available in pysap, more specifically a 2D parallel MRI brain slice on 32 channels and the acquisition cartesian scheme.

Package import

from mri.operators import FFT, WaveletN
from mri.reconstructors import CalibrationlessReconstructor
from pysap.data import get_sample_data

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

Loading input data

cartesian_ref_image = get_sample_data('2d-pmri').data
image = np.linalg.norm(cartesian_ref_image, axis=0)
# Obtain MRI cartesian mask
mask = get_sample_data("cartesian-mri-mask").data

View Input

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()

Out:

/volatile/Chaithya/actions-runner/_work/pysap/pysap/examples/pysap-mri/CalibrationlessReconstruction_GL_Cartesian.py:43: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()

Generate the kspace

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

# Get the locations of the kspace samples and the associated observations
fourier_op = FFT(mask=mask, shape=image.shape,
                 n_coils=cartesian_ref_image.shape[0])
kspace_obs = fourier_op.op(cartesian_ref_image)

Zero order solution

zero_soln = np.linalg.norm(fourier_op.adj_op(kspace_obs), axis=0)
base_ssim = ssim(zero_soln, image)
plt.imshow(np.abs(zero_soln), cmap='gray')
plt.title('Zero Order Solution : SSIM = ' + str(np.around(base_ssim, 2)))
plt.show()

Out:

/volatile/Chaithya/actions-runner/_work/pysap/pysap/examples/pysap-mri/CalibrationlessReconstruction_GL_Cartesian.py:64: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  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

# Setup the operators
linear_op = WaveletN(
    wavelet_name='sym8',
    nb_scale=4,
    n_coils=cartesian_ref_image.shape[0],
)
coeffs = linear_op.op(cartesian_ref_image)
regularizer_op = GroupLASSO(weights=6e-8)
# Setup Reconstructor
reconstructor = CalibrationlessReconstructor(
    fourier_op=fourier_op,
    linear_op=linear_op,
    regularizer_op=regularizer_op,
    gradient_formulation='synthesis',
    verbose=1,
)

Out:

WARNING: Making input data immutable.
Lipschitz constant is 1.0999999344348907
The lipschitz constraint is satisfied

Run the FISTA reconstruction and view results

image_rec, costs, metrics = reconstructor.reconstruct(
    kspace_data=kspace_obs,
    optimization_alg='fista',
    num_iterations=100,
)
image_rec = np.linalg.norm(image_rec, axis=0)
recon_ssim = ssim(image_rec, image)
plt.imshow(np.abs(image_rec), cmap='gray')
plt.title('Iterative Reconstruction : SSIM = ' + str(np.around(recon_ssim, 2)))
plt.show()

Out:

WARNING: Making input data immutable.
 - mu:  6e-08
 - lipschitz constant:  1.0999999344348907
 - data:  (512, 512)
 - wavelet:  <mri.operators.linear.wavelet.WaveletN object at 0x7f3a681c9cf0> - 4
 - max iterations:  100
 - image variable shape:  (512, 512)
 - alpha variable shape:  (32, 291721)
----------------------------------------
Starting optimization...

  0%|          | 0/100 [00:00<?, ?it/s]
  1%|1         | 1/100 [00:02<03:33,  2.15s/it]
  2%|2         | 2/100 [00:04<03:30,  2.15s/it]
  3%|3         | 3/100 [00:06<03:28,  2.15s/it]
  4%|4         | 4/100 [00:08<03:25,  2.14s/it]
  5%|5         | 5/100 [00:10<03:23,  2.15s/it]
  6%|6         | 6/100 [00:12<03:21,  2.15s/it]
  7%|7         | 7/100 [00:15<03:19,  2.14s/it]
  8%|8         | 8/100 [00:17<03:17,  2.14s/it]
  9%|9         | 9/100 [00:19<03:14,  2.14s/it]
 10%|#         | 10/100 [00:21<03:12,  2.14s/it]
 11%|#1        | 11/100 [00:23<03:09,  2.13s/it]
 12%|#2        | 12/100 [00:25<03:07,  2.13s/it]
 13%|#3        | 13/100 [00:27<03:05,  2.13s/it]
 14%|#4        | 14/100 [00:29<03:03,  2.13s/it]
 15%|#5        | 15/100 [00:32<03:01,  2.14s/it]
 16%|#6        | 16/100 [00:34<02:59,  2.13s/it]
 17%|#7        | 17/100 [00:36<02:57,  2.13s/it]
 18%|#8        | 18/100 [00:38<02:54,  2.13s/it]
 19%|#9        | 19/100 [00:40<02:52,  2.13s/it]
 20%|##        | 20/100 [00:42<02:50,  2.13s/it]
 21%|##1       | 21/100 [00:44<02:48,  2.13s/it]
 22%|##2       | 22/100 [00:46<02:46,  2.13s/it]
 23%|##3       | 23/100 [00:49<02:43,  2.13s/it]
 24%|##4       | 24/100 [00:51<02:41,  2.13s/it]
 25%|##5       | 25/100 [00:53<02:39,  2.13s/it]
 26%|##6       | 26/100 [00:55<02:37,  2.13s/it]
 27%|##7       | 27/100 [00:57<02:35,  2.13s/it]
 28%|##8       | 28/100 [00:59<02:33,  2.13s/it]
 29%|##9       | 29/100 [01:01<02:31,  2.14s/it]
 30%|###       | 30/100 [01:04<02:29,  2.14s/it]
 31%|###1      | 31/100 [01:06<02:27,  2.14s/it]
 32%|###2      | 32/100 [01:08<02:25,  2.14s/it]
 33%|###3      | 33/100 [01:10<02:22,  2.13s/it]
 34%|###4      | 34/100 [01:12<02:20,  2.13s/it]
 35%|###5      | 35/100 [01:14<02:18,  2.13s/it]
 36%|###6      | 36/100 [01:16<02:16,  2.13s/it]
 37%|###7      | 37/100 [01:18<02:14,  2.13s/it]
 38%|###8      | 38/100 [01:21<02:12,  2.13s/it]
 39%|###9      | 39/100 [01:23<02:10,  2.13s/it]
 40%|####      | 40/100 [01:25<02:07,  2.13s/it]
 41%|####1     | 41/100 [01:27<02:05,  2.13s/it]
 42%|####2     | 42/100 [01:29<02:03,  2.13s/it]
 43%|####3     | 43/100 [01:31<02:01,  2.13s/it]
 44%|####4     | 44/100 [01:33<01:59,  2.13s/it]
 45%|####5     | 45/100 [01:36<01:57,  2.13s/it]
 46%|####6     | 46/100 [01:38<01:54,  2.13s/it]
 47%|####6     | 47/100 [01:40<01:53,  2.13s/it]
 48%|####8     | 48/100 [01:42<01:50,  2.13s/it]
 49%|####9     | 49/100 [01:44<01:48,  2.13s/it]
 50%|#####     | 50/100 [01:46<01:46,  2.13s/it]
 51%|#####1    | 51/100 [01:48<01:44,  2.13s/it]
 52%|#####2    | 52/100 [01:50<01:42,  2.13s/it]
 53%|#####3    | 53/100 [01:53<01:39,  2.13s/it]
 54%|#####4    | 54/100 [01:55<01:37,  2.12s/it]
 55%|#####5    | 55/100 [01:57<01:35,  2.12s/it]
 56%|#####6    | 56/100 [01:59<01:33,  2.13s/it]
 57%|#####6    | 57/100 [02:01<01:31,  2.13s/it]
 58%|#####8    | 58/100 [02:03<01:29,  2.13s/it]
 59%|#####8    | 59/100 [02:05<01:27,  2.13s/it]
 60%|######    | 60/100 [02:07<01:25,  2.13s/it]
 61%|######1   | 61/100 [02:10<01:23,  2.13s/it]
 62%|######2   | 62/100 [02:12<01:20,  2.13s/it]
 63%|######3   | 63/100 [02:14<01:18,  2.13s/it]
 64%|######4   | 64/100 [02:16<01:16,  2.13s/it]
 65%|######5   | 65/100 [02:18<01:14,  2.13s/it]
 66%|######6   | 66/100 [02:20<01:12,  2.13s/it]
 67%|######7   | 67/100 [02:22<01:10,  2.13s/it]
 68%|######8   | 68/100 [02:24<01:08,  2.13s/it]
 69%|######9   | 69/100 [02:27<01:06,  2.13s/it]
 70%|#######   | 70/100 [02:29<01:03,  2.13s/it]
 71%|#######1  | 71/100 [02:31<01:01,  2.13s/it]
 72%|#######2  | 72/100 [02:33<00:59,  2.13s/it]
 73%|#######3  | 73/100 [02:35<00:57,  2.13s/it]
 74%|#######4  | 74/100 [02:37<00:55,  2.13s/it]
 75%|#######5  | 75/100 [02:39<00:53,  2.13s/it]
 76%|#######6  | 76/100 [02:42<00:51,  2.13s/it]
 77%|#######7  | 77/100 [02:44<00:49,  2.13s/it]
 78%|#######8  | 78/100 [02:46<00:46,  2.13s/it]
 79%|#######9  | 79/100 [02:48<00:44,  2.13s/it]
 80%|########  | 80/100 [02:50<00:42,  2.13s/it]
 81%|########1 | 81/100 [02:52<00:40,  2.13s/it]
 82%|########2 | 82/100 [02:54<00:38,  2.13s/it]
 83%|########2 | 83/100 [02:56<00:36,  2.13s/it]
 84%|########4 | 84/100 [02:59<00:34,  2.13s/it]
 85%|########5 | 85/100 [03:01<00:31,  2.13s/it]
 86%|########6 | 86/100 [03:03<00:29,  2.13s/it]
 87%|########7 | 87/100 [03:05<00:27,  2.13s/it]
 88%|########8 | 88/100 [03:07<00:25,  2.13s/it]
 89%|########9 | 89/100 [03:09<00:23,  2.13s/it]
 90%|######### | 90/100 [03:11<00:21,  2.13s/it]
 91%|#########1| 91/100 [03:14<00:19,  2.13s/it]
 92%|#########2| 92/100 [03:16<00:17,  2.13s/it]
 93%|#########3| 93/100 [03:18<00:14,  2.13s/it]
 94%|#########3| 94/100 [03:20<00:12,  2.13s/it]
 95%|#########5| 95/100 [03:22<00:10,  2.13s/it]
 96%|#########6| 96/100 [03:24<00:08,  2.13s/it]
 97%|#########7| 97/100 [03:26<00:06,  2.13s/it]
 98%|#########8| 98/100 [03:28<00:04,  2.14s/it]
 99%|#########9| 99/100 [03:31<00:02,  2.14s/it]
100%|##########| 100/100 [03:33<00:00,  2.13s/it]
100%|##########| 100/100 [03:33<00:00,  2.13s/it]
 - final iteration number:  100
 - final log10 cost value:  6.0
 - converged:  False
Done.
Execution time:  213.2270079846494  seconds
----------------------------------------
/volatile/Chaithya/actions-runner/_work/pysap/pysap/examples/pysap-mri/CalibrationlessReconstruction_GL_Cartesian.py:100: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()