Cartesian Self Calibrating Reconstruction
Contents
Note
Click here to download the full example code or to run this example in your browser via Binder
Cartesian Self Calibrating Reconstruction#
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 SelfCalibrationReconstructor
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
cartesian_ref_image = get_sample_data('2d-pmri')
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/SelfCalibrationReconstruction_Cartesian.py:44: 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/SelfCalibrationReconstruction_Cartesian.py:65: 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,
)
regularizer_op = SparseThreshold(Identity(), 1.5e-8, thresh_type="soft")
# Setup Reconstructor
reconstructor = SelfCalibrationReconstructor(
fourier_op=fourier_op,
linear_op=linear_op,
regularizer_op=regularizer_op,
gradient_formulation='synthesis',
kspace_portion=0.01,
verbose=1,
)
Run the FISTA reconstruction and view results
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('Iterative Reconstruction : SSIM = ' + str(np.around(recon_ssim, 2)))
plt.show()
Out:
/volatile/Chaithya/actions-runner/_work/_tool/Python/3.10.12/x64/lib/python3.10/site-packages/mri/reconstructors/utils/extract_sensitivity_maps.py:86: UserWarning: Data Values seem to have rank 3 (>2). Using is_fft for now.
warnings.warn('Data Values seem to have rank '
WARNING: Making input data immutable.
Lipschitz constant is 1.0983815952872344
The lipschitz constraint is satisfied
WARNING: Making input data immutable.
- mu: 1.5e-08
- lipschitz constant: 1.0983815952872344
- data: (512, 512)
- wavelet: <mri.operators.linear.wavelet.WaveletN object at 0x7f3a4e93ae00> - 4
- max iterations: 200
- image variable shape: (512, 512)
- alpha variable shape: (291721,)
----------------------------------------
Starting optimization...
0%| | 0/200 [00:00<?, ?it/s]
0%| | 1/200 [00:01<03:24, 1.03s/it]
1%|1 | 2/200 [00:02<03:20, 1.01s/it]
2%|1 | 3/200 [00:03<03:17, 1.00s/it]
2%|2 | 4/200 [00:04<03:15, 1.00it/s]
2%|2 | 5/200 [00:04<03:14, 1.00it/s]
3%|3 | 6/200 [00:05<03:12, 1.01it/s]
4%|3 | 7/200 [00:06<03:11, 1.01it/s]
4%|4 | 8/200 [00:07<03:11, 1.00it/s]
4%|4 | 9/200 [00:08<03:10, 1.00it/s]
5%|5 | 10/200 [00:09<03:08, 1.01it/s]
6%|5 | 11/200 [00:10<03:07, 1.01it/s]
6%|6 | 12/200 [00:11<03:06, 1.01it/s]
6%|6 | 13/200 [00:12<03:05, 1.01it/s]
7%|7 | 14/200 [00:13<03:05, 1.01it/s]
8%|7 | 15/200 [00:14<03:04, 1.00it/s]
8%|8 | 16/200 [00:15<03:03, 1.00it/s]
8%|8 | 17/200 [00:16<03:02, 1.00it/s]
9%|9 | 18/200 [00:17<03:01, 1.00it/s]
10%|9 | 19/200 [00:18<03:00, 1.00it/s]
10%|# | 20/200 [00:19<02:59, 1.00it/s]
10%|# | 21/200 [00:20<02:58, 1.00it/s]
11%|#1 | 22/200 [00:21<02:57, 1.00it/s]
12%|#1 | 23/200 [00:22<02:56, 1.00it/s]
12%|#2 | 24/200 [00:23<02:55, 1.00it/s]
12%|#2 | 25/200 [00:24<02:54, 1.00it/s]
13%|#3 | 26/200 [00:25<02:53, 1.00it/s]
14%|#3 | 27/200 [00:26<02:52, 1.00it/s]
14%|#4 | 28/200 [00:27<02:51, 1.00it/s]
14%|#4 | 29/200 [00:28<02:50, 1.00it/s]
15%|#5 | 30/200 [00:29<02:49, 1.00it/s]
16%|#5 | 31/200 [00:30<02:48, 1.00it/s]
16%|#6 | 32/200 [00:31<02:47, 1.00it/s]
16%|#6 | 33/200 [00:32<02:46, 1.00it/s]
17%|#7 | 34/200 [00:33<02:45, 1.00it/s]
18%|#7 | 35/200 [00:34<02:44, 1.00it/s]
18%|#8 | 36/200 [00:35<02:43, 1.00it/s]
18%|#8 | 37/200 [00:36<02:42, 1.00it/s]
19%|#9 | 38/200 [00:37<02:41, 1.00it/s]
20%|#9 | 39/200 [00:38<02:40, 1.00it/s]
20%|## | 40/200 [00:39<02:39, 1.00it/s]
20%|## | 41/200 [00:40<02:38, 1.00it/s]
21%|##1 | 42/200 [00:41<02:37, 1.00it/s]
22%|##1 | 43/200 [00:42<02:36, 1.00it/s]
22%|##2 | 44/200 [00:43<02:35, 1.00it/s]
22%|##2 | 45/200 [00:44<02:34, 1.00it/s]
23%|##3 | 46/200 [00:45<02:33, 1.00it/s]
24%|##3 | 47/200 [00:46<02:33, 1.00s/it]
24%|##4 | 48/200 [00:47<02:32, 1.00s/it]
24%|##4 | 49/200 [00:48<02:30, 1.00it/s]
25%|##5 | 50/200 [00:49<02:29, 1.00it/s]
26%|##5 | 51/200 [00:50<02:28, 1.00it/s]
26%|##6 | 52/200 [00:51<02:27, 1.00it/s]
26%|##6 | 53/200 [00:52<02:26, 1.00it/s]
27%|##7 | 54/200 [00:53<02:25, 1.00it/s]
28%|##7 | 55/200 [00:54<02:24, 1.00it/s]
28%|##8 | 56/200 [00:55<02:23, 1.01it/s]
28%|##8 | 57/200 [00:56<02:22, 1.00it/s]
29%|##9 | 58/200 [00:57<02:21, 1.01it/s]
30%|##9 | 59/200 [00:58<02:20, 1.00it/s]
30%|### | 60/200 [00:59<02:19, 1.01it/s]
30%|### | 61/200 [01:00<02:18, 1.00it/s]
31%|###1 | 62/200 [01:01<02:17, 1.01it/s]
32%|###1 | 63/200 [01:02<02:16, 1.00it/s]
32%|###2 | 64/200 [01:03<02:15, 1.01it/s]
32%|###2 | 65/200 [01:04<02:14, 1.01it/s]
33%|###3 | 66/200 [01:05<02:13, 1.01it/s]
34%|###3 | 67/200 [01:06<02:12, 1.00it/s]
34%|###4 | 68/200 [01:07<02:11, 1.00it/s]
34%|###4 | 69/200 [01:08<02:10, 1.00it/s]
35%|###5 | 70/200 [01:09<02:09, 1.00it/s]
36%|###5 | 71/200 [01:10<02:08, 1.00it/s]
36%|###6 | 72/200 [01:11<02:07, 1.00it/s]
36%|###6 | 73/200 [01:12<02:06, 1.00it/s]
37%|###7 | 74/200 [01:13<02:05, 1.00it/s]
38%|###7 | 75/200 [01:14<02:04, 1.00it/s]
38%|###8 | 76/200 [01:15<02:03, 1.00it/s]
38%|###8 | 77/200 [01:16<02:02, 1.00it/s]
39%|###9 | 78/200 [01:17<02:01, 1.00it/s]
40%|###9 | 79/200 [01:18<02:00, 1.01it/s]
40%|#### | 80/200 [01:19<01:59, 1.00it/s]
40%|#### | 81/200 [01:20<01:58, 1.00it/s]
41%|####1 | 82/200 [01:21<01:57, 1.00it/s]
42%|####1 | 83/200 [01:22<01:56, 1.00it/s]
42%|####2 | 84/200 [01:23<01:55, 1.00it/s]
42%|####2 | 85/200 [01:24<01:54, 1.00it/s]
43%|####3 | 86/200 [01:25<01:53, 1.00it/s]
44%|####3 | 87/200 [01:26<01:52, 1.00it/s]
44%|####4 | 88/200 [01:27<01:51, 1.00it/s]
44%|####4 | 89/200 [01:28<01:50, 1.00it/s]
45%|####5 | 90/200 [01:29<01:49, 1.00it/s]
46%|####5 | 91/200 [01:30<01:48, 1.00it/s]
46%|####6 | 92/200 [01:31<01:47, 1.00it/s]
46%|####6 | 93/200 [01:32<01:46, 1.00it/s]
47%|####6 | 94/200 [01:33<01:45, 1.00it/s]
48%|####7 | 95/200 [01:34<01:44, 1.00it/s]
48%|####8 | 96/200 [01:35<01:43, 1.00it/s]
48%|####8 | 97/200 [01:36<01:42, 1.00it/s]
49%|####9 | 98/200 [01:37<01:41, 1.00it/s]
50%|####9 | 99/200 [01:38<01:40, 1.00it/s]
50%|##### | 100/200 [01:39<01:39, 1.00it/s]
50%|##### | 101/200 [01:40<01:38, 1.00it/s]
51%|#####1 | 102/200 [01:41<01:37, 1.00it/s]
52%|#####1 | 103/200 [01:42<01:36, 1.00it/s]
52%|#####2 | 104/200 [01:43<01:35, 1.00it/s]
52%|#####2 | 105/200 [01:44<01:35, 1.00s/it]
53%|#####3 | 106/200 [01:45<01:33, 1.00it/s]
54%|#####3 | 107/200 [01:46<01:32, 1.00it/s]
54%|#####4 | 108/200 [01:47<01:31, 1.00it/s]
55%|#####4 | 109/200 [01:48<01:30, 1.00it/s]
55%|#####5 | 110/200 [01:49<01:29, 1.00it/s]
56%|#####5 | 111/200 [01:50<01:28, 1.00it/s]
56%|#####6 | 112/200 [01:51<01:27, 1.00it/s]
56%|#####6 | 113/200 [01:52<01:26, 1.00it/s]
57%|#####6 | 114/200 [01:53<01:25, 1.00it/s]
57%|#####7 | 115/200 [01:54<01:24, 1.00it/s]
58%|#####8 | 116/200 [01:55<01:23, 1.00it/s]
58%|#####8 | 117/200 [01:56<01:22, 1.00it/s]
59%|#####8 | 118/200 [01:57<01:21, 1.00it/s]
60%|#####9 | 119/200 [01:58<01:20, 1.00it/s]
60%|###### | 120/200 [01:59<01:19, 1.00it/s]
60%|###### | 121/200 [02:00<01:18, 1.00it/s]
61%|######1 | 122/200 [02:01<01:17, 1.00it/s]
62%|######1 | 123/200 [02:02<01:16, 1.00it/s]
62%|######2 | 124/200 [02:03<01:15, 1.00it/s]
62%|######2 | 125/200 [02:04<01:14, 1.00it/s]
63%|######3 | 126/200 [02:05<01:13, 1.00it/s]
64%|######3 | 127/200 [02:06<01:12, 1.00it/s]
64%|######4 | 128/200 [02:07<01:11, 1.00it/s]
64%|######4 | 129/200 [02:08<01:10, 1.01it/s]
65%|######5 | 130/200 [02:09<01:09, 1.00it/s]
66%|######5 | 131/200 [02:10<01:08, 1.00it/s]
66%|######6 | 132/200 [02:11<01:07, 1.00it/s]
66%|######6 | 133/200 [02:12<01:06, 1.00it/s]
67%|######7 | 134/200 [02:13<01:05, 1.00it/s]
68%|######7 | 135/200 [02:14<01:04, 1.00it/s]
68%|######8 | 136/200 [02:15<01:03, 1.00it/s]
68%|######8 | 137/200 [02:16<01:02, 1.00it/s]
69%|######9 | 138/200 [02:17<01:01, 1.01it/s]
70%|######9 | 139/200 [02:18<01:00, 1.01it/s]
70%|####### | 140/200 [02:19<00:59, 1.00it/s]
70%|####### | 141/200 [02:20<00:58, 1.00it/s]
71%|#######1 | 142/200 [02:21<00:57, 1.00it/s]
72%|#######1 | 143/200 [02:22<00:56, 1.00it/s]
72%|#######2 | 144/200 [02:23<00:55, 1.00it/s]
72%|#######2 | 145/200 [02:24<00:54, 1.00it/s]
73%|#######3 | 146/200 [02:25<00:53, 1.00it/s]
74%|#######3 | 147/200 [02:26<00:52, 1.00it/s]
74%|#######4 | 148/200 [02:27<00:51, 1.00it/s]
74%|#######4 | 149/200 [02:28<00:50, 1.00it/s]
75%|#######5 | 150/200 [02:29<00:50, 1.00s/it]
76%|#######5 | 151/200 [02:30<00:48, 1.00it/s]
76%|#######6 | 152/200 [02:31<00:47, 1.00it/s]
76%|#######6 | 153/200 [02:32<00:46, 1.00it/s]
77%|#######7 | 154/200 [02:33<00:45, 1.00it/s]
78%|#######7 | 155/200 [02:34<00:45, 1.00s/it]
78%|#######8 | 156/200 [02:35<00:43, 1.00it/s]
78%|#######8 | 157/200 [02:36<00:42, 1.00it/s]
79%|#######9 | 158/200 [02:37<00:41, 1.00it/s]
80%|#######9 | 159/200 [02:38<00:41, 1.00s/it]
80%|######## | 160/200 [02:39<00:39, 1.00it/s]
80%|######## | 161/200 [02:40<00:38, 1.00it/s]
81%|########1 | 162/200 [02:41<00:37, 1.00it/s]
82%|########1 | 163/200 [02:42<00:36, 1.00it/s]
82%|########2 | 164/200 [02:43<00:35, 1.00it/s]
82%|########2 | 165/200 [02:44<00:34, 1.00it/s]
83%|########2 | 166/200 [02:45<00:33, 1.00it/s]
84%|########3 | 167/200 [02:46<00:32, 1.00it/s]
84%|########4 | 168/200 [02:47<00:31, 1.00it/s]
84%|########4 | 169/200 [02:48<00:30, 1.00it/s]
85%|########5 | 170/200 [02:49<00:29, 1.00it/s]
86%|########5 | 171/200 [02:50<00:28, 1.00it/s]
86%|########6 | 172/200 [02:51<00:27, 1.00it/s]
86%|########6 | 173/200 [02:52<00:26, 1.00it/s]
87%|########7 | 174/200 [02:53<00:25, 1.00it/s]
88%|########7 | 175/200 [02:54<00:24, 1.00it/s]
88%|########8 | 176/200 [02:55<00:23, 1.00it/s]
88%|########8 | 177/200 [02:56<00:22, 1.00it/s]
89%|########9 | 178/200 [02:57<00:21, 1.00it/s]
90%|########9 | 179/200 [02:58<00:20, 1.00it/s]
90%|######### | 180/200 [02:59<00:19, 1.00it/s]
90%|######### | 181/200 [03:00<00:19, 1.00s/it]
91%|#########1| 182/200 [03:01<00:18, 1.00s/it]
92%|#########1| 183/200 [03:02<00:17, 1.00s/it]
92%|#########2| 184/200 [03:03<00:15, 1.00it/s]
92%|#########2| 185/200 [03:04<00:14, 1.00it/s]
93%|#########3| 186/200 [03:05<00:13, 1.00it/s]
94%|#########3| 187/200 [03:06<00:12, 1.00it/s]
94%|#########3| 188/200 [03:07<00:11, 1.00it/s]
94%|#########4| 189/200 [03:08<00:10, 1.00it/s]
95%|#########5| 190/200 [03:09<00:10, 1.00s/it]
96%|#########5| 191/200 [03:10<00:09, 1.00s/it]
96%|#########6| 192/200 [03:11<00:08, 1.00s/it]
96%|#########6| 193/200 [03:12<00:07, 1.00s/it]
97%|#########7| 194/200 [03:13<00:05, 1.00it/s]
98%|#########7| 195/200 [03:14<00:04, 1.00it/s]
98%|#########8| 196/200 [03:15<00:03, 1.00it/s]
98%|#########8| 197/200 [03:16<00:02, 1.00it/s]
99%|#########9| 198/200 [03:17<00:01, 1.00it/s]
100%|#########9| 199/200 [03:18<00:00, 1.00it/s]
100%|##########| 200/200 [03:19<00:00, 1.00it/s]
100%|##########| 200/200 [03:19<00:00, 1.00it/s]
- final iteration number: 200
- final log10 cost value: 6.0
- converged: False
Done.
Execution time: 199.45882562408224 seconds
----------------------------------------
/volatile/Chaithya/actions-runner/_work/pysap/pysap/examples/pysap-mri/SelfCalibrationReconstruction_Cartesian.py:100: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
plt.show()
Total running time of the script: ( 3 minutes 51.225 seconds)
