3D Neuroimaging cartesian reconstruction
Contents
Note
Click here to download the full example code or to run this example in your browser via Binder
3D Neuroimaging cartesian reconstruction#
Author: LElgueddari
In this tutorial we will reconstruct an MRI image from the sparse 3D kspace measurements.
Import neuroimaging data#
We use the toy datasets available in pysap, more specifically the 3D orange and the cartesian acquisition scheme.
Package import
from modopt.math.metrics import ssim
from mri.operators import FFT, WaveletN
from mri.operators.utils import convert_mask_to_locations
from mri.reconstructors import SingleChannelReconstructor
from pysap.data import get_sample_data
import pysap
# Third party import
from modopt.opt.linear import Identity
from modopt.opt.proximity import SparseThreshold
import numpy as np
import matplotlib.pyplot as plt
Loading input data and convert it into a single channel using Sum-Of-Squares
image = get_sample_data('3d-pmri').data
image = np.linalg.norm(image, axis=0)
# Obtain K-Space Cartesian Mask (straight line readout along z)
mask = get_sample_data("2d-poisson-disk-mask").data
mask = np.repeat(np.expand_dims(mask, axis=-1), image.shape[-1],
axis=-1)
View Input
plt.subplot(1, 2, 1)
plt.imshow(np.abs(image[..., 80]), cmap='gray')
plt.title("MRI Data")
plt.subplot(1, 2, 2)
plt.imshow(mask[..., 80], cmap='gray')
plt.title("K-space Sampling Mask")
plt.show()
Out:
/volatile/Chaithya/actions-runner/_work/pysap/pysap/examples/pysap-mri/cartesian_reconstruction_3d.py:49: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
plt.show()
Generate the kspace#
From the 3D Orange volume 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
# Generate the subsampled kspace
fourier_op = FFT(mask=mask, shape=image.shape)
kspace_data = fourier_op.op(image)
Zero order solution
image_rec0 = fourier_op.adj_op(kspace_data)
base_ssim = ssim(image_rec0, image)
plt.imshow(np.abs(image_rec0[..., 80]), cmap='gray')
plt.title('Gridded solution : SSIM = ' + str(np.around(base_ssim, 2)))
plt.show()
Out:
/volatile/Chaithya/actions-runner/_work/pysap/pysap/examples/pysap-mri/cartesian_reconstruction_3d.py:70: 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_scales=4,
dim=3,
padding_mode="periodization",
)
regularizer_op = SparseThreshold(Identity(), 2 * 1e-11, thresh_type="soft")
# Setup Reconstructor
reconstructor = SingleChannelReconstructor(
fourier_op=fourier_op,
linear_op=linear_op,
regularizer_op=regularizer_op,
gradient_formulation='synthesis',
verbose=1,
)
Out:
/volatile/Chaithya/actions-runner/_work/_tool/Python/3.10.12/x64/lib/python3.10/site-packages/pywt/_multilevel.py:43: UserWarning: Level value of 4 is too high: all coefficients will experience boundary effects.
warnings.warn(
WARNING: Making input data immutable.
Lipschitz constant is 1.1
The lipschitz constraint is satisfied
Start Reconstruction
image_rec, costs, metrics = reconstructor.reconstruct(
kspace_data=kspace_data,
optimization_alg='fista',
num_iterations=200,
)
recon_ssim = ssim(image_rec, image)
plt.imshow(np.abs(image_rec[..., 80]), cmap='gray')
plt.title('Iterative Reconstruction : SSIM = ' + str(np.around(recon_ssim, 2)))
plt.show()
Out:
WARNING: Making input data immutable.
/volatile/Chaithya/actions-runner/_work/_tool/Python/3.10.12/x64/lib/python3.10/site-packages/pywt/_multilevel.py:43: UserWarning: Level value of 4 is too high: all coefficients will experience boundary effects.
warnings.warn(
- mu: 2e-11
- lipschitz constant: 1.1
- data: (128, 128, 160)
- wavelet: <mri.operators.linear.wavelet.WaveletN object at 0x7f3a4e96cb80> - 4
- max iterations: 200
- image variable shape: (128, 128, 160)
- alpha variable shape: (2621440,)
----------------------------------------
Starting optimization...
0%| | 0/200 [00:00<?, ?it/s]
0%| | 1/200 [00:01<03:43, 1.13s/it]
1%|1 | 2/200 [00:02<03:49, 1.16s/it]
2%|1 | 3/200 [00:03<03:51, 1.18s/it]
2%|2 | 4/200 [00:04<03:47, 1.16s/it]
2%|2 | 5/200 [00:05<03:42, 1.14s/it]
3%|3 | 6/200 [00:06<03:41, 1.14s/it]
4%|3 | 7/200 [00:08<03:40, 1.14s/it]
4%|4 | 8/200 [00:09<03:39, 1.15s/it]
4%|4 | 9/200 [00:10<03:39, 1.15s/it]
5%|5 | 10/200 [00:11<03:39, 1.16s/it]
6%|5 | 11/200 [00:12<03:35, 1.14s/it]
6%|6 | 12/200 [00:13<03:32, 1.13s/it]
6%|6 | 13/200 [00:14<03:34, 1.15s/it]
7%|7 | 14/200 [00:16<03:36, 1.16s/it]
8%|7 | 15/200 [00:17<03:35, 1.17s/it]
8%|8 | 16/200 [00:18<03:30, 1.15s/it]
8%|8 | 17/200 [00:19<03:28, 1.14s/it]
9%|9 | 18/200 [00:20<03:30, 1.16s/it]
10%|9 | 19/200 [00:21<03:30, 1.16s/it]
10%|# | 20/200 [00:23<03:28, 1.16s/it]
10%|# | 21/200 [00:24<03:26, 1.15s/it]
11%|#1 | 22/200 [00:25<03:23, 1.15s/it]
12%|#1 | 23/200 [00:26<03:20, 1.13s/it]
12%|#2 | 24/200 [00:27<03:20, 1.14s/it]
12%|#2 | 25/200 [00:28<03:21, 1.15s/it]
13%|#3 | 26/200 [00:29<03:22, 1.16s/it]
14%|#3 | 27/200 [00:31<03:20, 1.16s/it]
14%|#4 | 28/200 [00:32<03:22, 1.18s/it]
14%|#4 | 29/200 [00:33<03:21, 1.18s/it]
15%|#5 | 30/200 [00:34<03:19, 1.18s/it]
16%|#5 | 31/200 [00:35<03:18, 1.17s/it]
16%|#6 | 32/200 [00:36<03:14, 1.15s/it]
16%|#6 | 33/200 [00:38<03:10, 1.14s/it]
17%|#7 | 34/200 [00:39<03:09, 1.14s/it]
18%|#7 | 35/200 [00:40<03:07, 1.14s/it]
18%|#8 | 36/200 [00:41<03:08, 1.15s/it]
18%|#8 | 37/200 [00:42<03:06, 1.14s/it]
19%|#9 | 38/200 [00:43<03:02, 1.13s/it]
20%|#9 | 39/200 [00:44<03:01, 1.12s/it]
20%|## | 40/200 [00:45<02:59, 1.12s/it]
20%|## | 41/200 [00:47<02:57, 1.12s/it]
21%|##1 | 42/200 [00:48<03:00, 1.14s/it]
22%|##1 | 43/200 [00:49<03:00, 1.15s/it]
22%|##2 | 44/200 [00:50<03:00, 1.16s/it]
22%|##2 | 45/200 [00:51<02:56, 1.14s/it]
23%|##3 | 46/200 [00:52<02:55, 1.14s/it]
24%|##3 | 47/200 [00:54<02:58, 1.17s/it]
24%|##4 | 48/200 [00:55<02:58, 1.18s/it]
24%|##4 | 49/200 [00:56<02:58, 1.18s/it]
25%|##5 | 50/200 [00:57<02:57, 1.18s/it]
26%|##5 | 51/200 [00:58<02:52, 1.16s/it]
26%|##6 | 52/200 [00:59<02:48, 1.14s/it]
26%|##6 | 53/200 [01:00<02:48, 1.15s/it]
27%|##7 | 54/200 [01:02<02:48, 1.16s/it]
28%|##7 | 55/200 [01:03<02:50, 1.18s/it]
28%|##8 | 56/200 [01:04<02:49, 1.18s/it]
28%|##8 | 57/200 [01:05<02:45, 1.16s/it]
29%|##9 | 58/200 [01:06<02:42, 1.14s/it]
30%|##9 | 59/200 [01:07<02:41, 1.15s/it]
30%|### | 60/200 [01:09<02:44, 1.17s/it]
30%|### | 61/200 [01:10<02:44, 1.18s/it]
31%|###1 | 62/200 [01:11<02:44, 1.19s/it]
32%|###1 | 63/200 [01:12<02:42, 1.19s/it]
32%|###2 | 64/200 [01:13<02:39, 1.17s/it]
32%|###2 | 65/200 [01:15<02:35, 1.15s/it]
33%|###3 | 66/200 [01:16<02:32, 1.14s/it]
34%|###3 | 67/200 [01:17<02:34, 1.16s/it]
34%|###4 | 68/200 [01:18<02:32, 1.16s/it]
34%|###4 | 69/200 [01:19<02:30, 1.15s/it]
35%|###5 | 70/200 [01:20<02:30, 1.16s/it]
36%|###5 | 71/200 [01:21<02:29, 1.16s/it]
36%|###6 | 72/200 [01:23<02:28, 1.16s/it]
36%|###6 | 73/200 [01:24<02:26, 1.16s/it]
37%|###7 | 74/200 [01:25<02:24, 1.15s/it]
38%|###7 | 75/200 [01:26<02:23, 1.14s/it]
38%|###8 | 76/200 [01:27<02:23, 1.15s/it]
38%|###8 | 77/200 [01:28<02:20, 1.14s/it]
39%|###9 | 78/200 [01:29<02:18, 1.13s/it]
40%|###9 | 79/200 [01:31<02:17, 1.13s/it]
40%|#### | 80/200 [01:32<02:17, 1.15s/it]
40%|#### | 81/200 [01:33<02:16, 1.15s/it]
41%|####1 | 82/200 [01:34<02:15, 1.15s/it]
42%|####1 | 83/200 [01:35<02:14, 1.15s/it]
42%|####2 | 84/200 [01:36<02:12, 1.15s/it]
42%|####2 | 85/200 [01:38<02:13, 1.16s/it]
43%|####3 | 86/200 [01:39<02:12, 1.16s/it]
44%|####3 | 87/200 [01:40<02:09, 1.14s/it]
44%|####4 | 88/200 [01:41<02:08, 1.14s/it]
44%|####4 | 89/200 [01:42<02:10, 1.17s/it]
45%|####5 | 90/200 [01:43<02:09, 1.18s/it]
46%|####5 | 91/200 [01:45<02:09, 1.19s/it]
46%|####6 | 92/200 [01:46<02:08, 1.19s/it]
46%|####6 | 93/200 [01:47<02:04, 1.16s/it]
47%|####6 | 94/200 [01:48<02:02, 1.16s/it]
48%|####7 | 95/200 [01:49<02:02, 1.17s/it]
48%|####8 | 96/200 [01:50<02:02, 1.18s/it]
48%|####8 | 97/200 [01:52<02:02, 1.19s/it]
49%|####9 | 98/200 [01:53<02:01, 1.19s/it]
50%|####9 | 99/200 [01:54<02:00, 1.19s/it]
50%|##### | 100/200 [01:55<01:56, 1.17s/it]
50%|##### | 101/200 [01:56<01:53, 1.15s/it]
51%|#####1 | 102/200 [01:57<01:53, 1.15s/it]
52%|#####1 | 103/200 [01:59<01:51, 1.15s/it]
52%|#####2 | 104/200 [02:00<01:49, 1.15s/it]
52%|#####2 | 105/200 [02:01<01:50, 1.16s/it]
53%|#####3 | 106/200 [02:02<01:48, 1.15s/it]
54%|#####3 | 107/200 [02:03<01:46, 1.14s/it]
54%|#####4 | 108/200 [02:04<01:44, 1.13s/it]
55%|#####4 | 109/200 [02:05<01:45, 1.16s/it]
55%|#####5 | 110/200 [02:07<01:43, 1.15s/it]
56%|#####5 | 111/200 [02:08<01:42, 1.15s/it]
56%|#####6 | 112/200 [02:09<01:41, 1.16s/it]
56%|#####6 | 113/200 [02:10<01:41, 1.16s/it]
57%|#####6 | 114/200 [02:11<01:40, 1.16s/it]
57%|#####7 | 115/200 [02:12<01:38, 1.16s/it]
58%|#####8 | 116/200 [02:14<01:36, 1.15s/it]
58%|#####8 | 117/200 [02:15<01:35, 1.15s/it]
59%|#####8 | 118/200 [02:16<01:34, 1.15s/it]
60%|#####9 | 119/200 [02:17<01:32, 1.14s/it]
60%|###### | 120/200 [02:18<01:30, 1.13s/it]
60%|###### | 121/200 [02:19<01:29, 1.13s/it]
61%|######1 | 122/200 [02:20<01:29, 1.15s/it]
62%|######1 | 123/200 [02:22<01:28, 1.14s/it]
62%|######2 | 124/200 [02:23<01:26, 1.14s/it]
62%|######2 | 125/200 [02:24<01:25, 1.14s/it]
63%|######3 | 126/200 [02:25<01:24, 1.14s/it]
64%|######3 | 127/200 [02:26<01:24, 1.16s/it]
64%|######4 | 128/200 [02:27<01:23, 1.16s/it]
64%|######4 | 129/200 [02:28<01:21, 1.14s/it]
65%|######5 | 130/200 [02:30<01:20, 1.15s/it]
66%|######5 | 131/200 [02:31<01:20, 1.17s/it]
66%|######6 | 132/200 [02:32<01:20, 1.18s/it]
66%|######6 | 133/200 [02:33<01:19, 1.19s/it]
67%|######7 | 134/200 [02:34<01:18, 1.18s/it]
68%|######7 | 135/200 [02:35<01:15, 1.16s/it]
68%|######8 | 136/200 [02:37<01:13, 1.14s/it]
68%|######8 | 137/200 [02:38<01:12, 1.15s/it]
69%|######9 | 138/200 [02:39<01:11, 1.16s/it]
70%|######9 | 139/200 [02:40<01:11, 1.18s/it]
70%|####### | 140/200 [02:41<01:10, 1.18s/it]
70%|####### | 141/200 [02:42<01:08, 1.16s/it]
71%|#######1 | 142/200 [02:44<01:06, 1.15s/it]
72%|#######1 | 143/200 [02:45<01:05, 1.15s/it]
72%|#######2 | 144/200 [02:46<01:05, 1.17s/it]
72%|#######2 | 145/200 [02:47<01:05, 1.18s/it]
73%|#######3 | 146/200 [02:48<01:04, 1.19s/it]
74%|#######3 | 147/200 [02:50<01:02, 1.19s/it]
74%|#######4 | 148/200 [02:51<01:00, 1.17s/it]
74%|#######4 | 149/200 [02:52<00:58, 1.15s/it]
75%|#######5 | 150/200 [02:53<00:57, 1.14s/it]
76%|#######5 | 151/200 [02:54<00:57, 1.17s/it]
76%|#######6 | 152/200 [02:55<00:55, 1.16s/it]
76%|#######6 | 153/200 [02:56<00:54, 1.15s/it]
77%|#######7 | 154/200 [02:58<00:53, 1.16s/it]
78%|#######7 | 155/200 [02:59<00:52, 1.16s/it]
78%|#######8 | 156/200 [03:00<00:51, 1.16s/it]
78%|#######8 | 157/200 [03:01<00:49, 1.16s/it]
79%|#######9 | 158/200 [03:02<00:48, 1.15s/it]
80%|#######9 | 159/200 [03:03<00:47, 1.15s/it]
80%|######## | 160/200 [03:05<00:46, 1.16s/it]
80%|######## | 161/200 [03:06<00:44, 1.14s/it]
81%|########1 | 162/200 [03:07<00:43, 1.14s/it]
82%|########1 | 163/200 [03:08<00:42, 1.14s/it]
82%|########2 | 164/200 [03:09<00:41, 1.15s/it]
82%|########2 | 165/200 [03:10<00:40, 1.15s/it]
83%|########2 | 166/200 [03:11<00:39, 1.15s/it]
84%|########3 | 167/200 [03:13<00:37, 1.15s/it]
84%|########4 | 168/200 [03:14<00:36, 1.15s/it]
84%|########4 | 169/200 [03:15<00:35, 1.16s/it]
85%|########5 | 170/200 [03:16<00:34, 1.16s/it]
86%|########5 | 171/200 [03:17<00:33, 1.15s/it]
86%|########6 | 172/200 [03:18<00:32, 1.15s/it]
86%|########6 | 173/200 [03:19<00:31, 1.17s/it]
87%|########7 | 174/200 [03:21<00:30, 1.18s/it]
88%|########7 | 175/200 [03:22<00:29, 1.19s/it]
88%|########8 | 176/200 [03:23<00:28, 1.19s/it]
88%|########8 | 177/200 [03:24<00:26, 1.17s/it]
89%|########9 | 178/200 [03:25<00:25, 1.15s/it]
90%|########9 | 179/200 [03:26<00:24, 1.16s/it]
90%|######### | 180/200 [03:28<00:23, 1.16s/it]
90%|######### | 181/200 [03:29<00:22, 1.18s/it]
91%|#########1| 182/200 [03:30<00:21, 1.18s/it]
92%|#########1| 183/200 [03:31<00:19, 1.16s/it]
92%|#########2| 184/200 [03:32<00:18, 1.15s/it]
92%|#########2| 185/200 [03:33<00:17, 1.15s/it]
93%|#########3| 186/200 [03:35<00:16, 1.18s/it]
94%|#########3| 187/200 [03:36<00:15, 1.18s/it]
94%|#########3| 188/200 [03:37<00:14, 1.19s/it]
94%|#########4| 189/200 [03:38<00:13, 1.19s/it]
95%|#########5| 190/200 [03:39<00:11, 1.17s/it]
96%|#########5| 191/200 [03:41<00:10, 1.15s/it]
96%|#########6| 192/200 [03:42<00:09, 1.14s/it]
96%|#########6| 193/200 [03:43<00:08, 1.16s/it]
97%|#########7| 194/200 [03:44<00:06, 1.16s/it]
98%|#########7| 195/200 [03:45<00:05, 1.15s/it]
98%|#########8| 196/200 [03:46<00:04, 1.16s/it]
98%|#########8| 197/200 [03:47<00:03, 1.16s/it]
99%|#########9| 198/200 [03:49<00:02, 1.16s/it]
100%|#########9| 199/200 [03:50<00:01, 1.16s/it]
100%|##########| 200/200 [03:51<00:00, 1.15s/it]
100%|##########| 200/200 [03:51<00:00, 1.16s/it]
- final iteration number: 200
- final log10 cost value: 6.0
- converged: False
Done.
Execution time: 231.41732934722677 seconds
----------------------------------------
/volatile/Chaithya/actions-runner/_work/pysap/pysap/examples/pysap-mri/cartesian_reconstruction_3d.py:105: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
plt.show()
Total running time of the script: ( 4 minutes 15.820 seconds)
