GPU Non-cartesian 2D reconstruction: Undecimated Wavelet
Contents
Note
Click here to download the full example code or to run this example in your browser via Binder
GPU Non-cartesian 2D reconstruction: Undecimated Wavelet#
Author: Chaithya G R
In this tutorial we will reconstruct an MRI image from non-cartesian kspace measurements.
Import neuroimaging data#
We use the toy datasets available in pysap, more specifically a 2D brain slice and the acquisition cartesian scheme.
Package import
from mri.operators import NonCartesianFFT, WaveletUD2
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
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()
Out:
/volatile/Chaithya/actions-runner/_work/pysap/pysap/examples/pysap-mri/non_cartesian_reconstruction.py:48: 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 radial 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 = NonCartesianFFT(samples=kspace_loc, shape=image.shape, implementation='gpuNUFFT')
kspace_obs = fourier_op.op(image)[0]
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, 2)))
plt.show()
Out:
/volatile/Chaithya/actions-runner/_work/pysap/pysap/examples/pysap-mri/non_cartesian_reconstruction.py:71: 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 = WaveletUD2(
wavelet_id=24,
nb_scale=4,
)
regularizer_op = SparseThreshold(Identity(), 2e-8, 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:
WARNING: Making input data immutable.
Lipschitz constant is 17.815955850645466
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=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:
WARNING: Making input data immutable.
- mu: 2e-08
- lipschitz constant: 17.815955850645466
- data: (512, 512)
- wavelet: <mri.operators.linear.wavelet.WaveletUD2 object at 0x7f3a4ede7610> - 4
- max iterations: 200
- image variable shape: (512, 512)
- alpha variable shape: (2621440,)
----------------------------------------
Starting optimization...
0%| | 0/200 [00:00<?, ?it/s]
0%| | 1/200 [00:01<05:34, 1.68s/it]
1%|1 | 2/200 [00:03<05:32, 1.68s/it]
2%|1 | 3/200 [00:05<05:30, 1.68s/it]
2%|2 | 4/200 [00:06<05:27, 1.67s/it]
2%|2 | 5/200 [00:08<05:25, 1.67s/it]
3%|3 | 6/200 [00:10<05:23, 1.67s/it]
4%|3 | 7/200 [00:11<05:21, 1.66s/it]
4%|4 | 8/200 [00:13<05:19, 1.66s/it]
4%|4 | 9/200 [00:15<05:17, 1.66s/it]
5%|5 | 10/200 [00:16<05:16, 1.66s/it]
6%|5 | 11/200 [00:18<05:15, 1.67s/it]
6%|6 | 12/200 [00:20<05:13, 1.67s/it]
6%|6 | 13/200 [00:21<05:11, 1.67s/it]
7%|7 | 14/200 [00:23<05:10, 1.67s/it]
8%|7 | 15/200 [00:25<05:08, 1.67s/it]
8%|8 | 16/200 [00:26<05:07, 1.67s/it]
8%|8 | 17/200 [00:28<05:05, 1.67s/it]
9%|9 | 18/200 [00:30<05:03, 1.67s/it]
10%|9 | 19/200 [00:31<05:01, 1.67s/it]
10%|# | 20/200 [00:33<05:00, 1.67s/it]
10%|# | 21/200 [00:35<04:58, 1.67s/it]
11%|#1 | 22/200 [00:36<04:57, 1.67s/it]
12%|#1 | 23/200 [00:38<04:55, 1.67s/it]
12%|#2 | 24/200 [00:40<04:54, 1.67s/it]
12%|#2 | 25/200 [00:41<04:52, 1.67s/it]
13%|#3 | 26/200 [00:43<04:50, 1.67s/it]
14%|#3 | 27/200 [00:45<04:49, 1.67s/it]
14%|#4 | 28/200 [00:46<04:49, 1.68s/it]
14%|#4 | 29/200 [00:48<04:47, 1.68s/it]
15%|#5 | 30/200 [00:50<04:46, 1.69s/it]
16%|#5 | 31/200 [00:51<04:44, 1.69s/it]
16%|#6 | 32/200 [00:53<04:43, 1.69s/it]
16%|#6 | 33/200 [00:55<04:41, 1.69s/it]
17%|#7 | 34/200 [00:56<04:39, 1.69s/it]
18%|#7 | 35/200 [00:58<04:39, 1.69s/it]
18%|#8 | 36/200 [01:00<04:37, 1.69s/it]
18%|#8 | 37/200 [01:01<04:35, 1.69s/it]
19%|#9 | 38/200 [01:03<04:33, 1.69s/it]
20%|#9 | 39/200 [01:05<04:31, 1.69s/it]
20%|## | 40/200 [01:07<04:30, 1.69s/it]
20%|## | 41/200 [01:08<04:28, 1.69s/it]
21%|##1 | 42/200 [01:10<04:26, 1.69s/it]
22%|##1 | 43/200 [01:12<04:24, 1.69s/it]
22%|##2 | 44/200 [01:13<04:23, 1.69s/it]
22%|##2 | 45/200 [01:15<04:21, 1.69s/it]
23%|##3 | 46/200 [01:17<04:20, 1.69s/it]
24%|##3 | 47/200 [01:18<04:17, 1.68s/it]
24%|##4 | 48/200 [01:20<04:14, 1.67s/it]
24%|##4 | 49/200 [01:22<04:13, 1.68s/it]
25%|##5 | 50/200 [01:23<04:12, 1.68s/it]
26%|##5 | 51/200 [01:25<04:11, 1.69s/it]
26%|##6 | 52/200 [01:27<04:10, 1.69s/it]
26%|##6 | 53/200 [01:28<04:08, 1.69s/it]
27%|##7 | 54/200 [01:30<04:07, 1.69s/it]
28%|##7 | 55/200 [01:32<04:05, 1.69s/it]
28%|##8 | 56/200 [01:34<04:03, 1.69s/it]
28%|##8 | 57/200 [01:35<04:01, 1.69s/it]
29%|##9 | 58/200 [01:37<03:59, 1.69s/it]
30%|##9 | 59/200 [01:39<03:58, 1.69s/it]
30%|### | 60/200 [01:40<03:56, 1.69s/it]
30%|### | 61/200 [01:42<03:54, 1.69s/it]
31%|###1 | 62/200 [01:44<03:53, 1.69s/it]
32%|###1 | 63/200 [01:45<03:51, 1.69s/it]
32%|###2 | 64/200 [01:47<03:50, 1.69s/it]
32%|###2 | 65/200 [01:49<03:48, 1.69s/it]
33%|###3 | 66/200 [01:50<03:46, 1.69s/it]
34%|###3 | 67/200 [01:52<03:44, 1.69s/it]
34%|###4 | 68/200 [01:54<03:42, 1.69s/it]
34%|###4 | 69/200 [01:56<03:41, 1.69s/it]
35%|###5 | 70/200 [01:57<03:39, 1.69s/it]
36%|###5 | 71/200 [01:59<03:38, 1.69s/it]
36%|###6 | 72/200 [02:01<03:36, 1.69s/it]
36%|###6 | 73/200 [02:02<03:35, 1.70s/it]
37%|###7 | 74/200 [02:04<03:33, 1.70s/it]
38%|###7 | 75/200 [02:06<03:32, 1.70s/it]
38%|###8 | 76/200 [02:07<03:30, 1.70s/it]
38%|###8 | 77/200 [02:09<03:29, 1.70s/it]
39%|###9 | 78/200 [02:11<03:27, 1.70s/it]
40%|###9 | 79/200 [02:13<03:26, 1.71s/it]
40%|#### | 80/200 [02:14<03:25, 1.71s/it]
40%|#### | 81/200 [02:16<03:24, 1.72s/it]
41%|####1 | 82/200 [02:18<03:22, 1.72s/it]
42%|####1 | 83/200 [02:19<03:21, 1.72s/it]
42%|####2 | 84/200 [02:21<03:19, 1.72s/it]
42%|####2 | 85/200 [02:23<03:18, 1.73s/it]
43%|####3 | 86/200 [02:25<03:17, 1.73s/it]
44%|####3 | 87/200 [02:26<03:15, 1.73s/it]
44%|####4 | 88/200 [02:28<03:14, 1.73s/it]
44%|####4 | 89/200 [02:30<03:12, 1.74s/it]
45%|####5 | 90/200 [02:32<03:10, 1.74s/it]
46%|####5 | 91/200 [02:33<03:08, 1.73s/it]
46%|####6 | 92/200 [02:35<03:06, 1.73s/it]
46%|####6 | 93/200 [02:37<03:05, 1.73s/it]
47%|####6 | 94/200 [02:38<03:03, 1.73s/it]
48%|####7 | 95/200 [02:40<03:01, 1.73s/it]
48%|####8 | 96/200 [02:42<03:00, 1.73s/it]
48%|####8 | 97/200 [02:44<02:58, 1.73s/it]
49%|####9 | 98/200 [02:45<02:56, 1.73s/it]
50%|####9 | 99/200 [02:47<02:54, 1.73s/it]
50%|##### | 100/200 [02:49<02:52, 1.73s/it]
50%|##### | 101/200 [02:51<02:50, 1.73s/it]
51%|#####1 | 102/200 [02:52<02:49, 1.73s/it]
52%|#####1 | 103/200 [02:54<02:48, 1.73s/it]
52%|#####2 | 104/200 [02:56<02:46, 1.73s/it]
52%|#####2 | 105/200 [02:58<02:44, 1.73s/it]
53%|#####3 | 106/200 [02:59<02:43, 1.73s/it]
54%|#####3 | 107/200 [03:01<02:41, 1.74s/it]
54%|#####4 | 108/200 [03:03<02:40, 1.74s/it]
55%|#####4 | 109/200 [03:05<02:38, 1.74s/it]
55%|#####5 | 110/200 [03:06<02:36, 1.74s/it]
56%|#####5 | 111/200 [03:08<02:34, 1.74s/it]
56%|#####6 | 112/200 [03:10<02:33, 1.74s/it]
56%|#####6 | 113/200 [03:11<02:31, 1.74s/it]
57%|#####6 | 114/200 [03:13<02:30, 1.75s/it]
57%|#####7 | 115/200 [03:15<02:28, 1.75s/it]
58%|#####8 | 116/200 [03:17<02:26, 1.75s/it]
58%|#####8 | 117/200 [03:18<02:24, 1.74s/it]
59%|#####8 | 118/200 [03:20<02:22, 1.74s/it]
60%|#####9 | 119/200 [03:22<02:21, 1.74s/it]
60%|###### | 120/200 [03:24<02:19, 1.74s/it]
60%|###### | 121/200 [03:25<02:17, 1.74s/it]
61%|######1 | 122/200 [03:27<02:15, 1.74s/it]
62%|######1 | 123/200 [03:29<02:14, 1.74s/it]
62%|######2 | 124/200 [03:31<02:12, 1.75s/it]
62%|######2 | 125/200 [03:32<02:11, 1.75s/it]
63%|######3 | 126/200 [03:34<02:09, 1.75s/it]
64%|######3 | 127/200 [03:36<02:07, 1.75s/it]
64%|######4 | 128/200 [03:38<02:05, 1.75s/it]
64%|######4 | 129/200 [03:39<02:03, 1.74s/it]
65%|######5 | 130/200 [03:41<02:01, 1.74s/it]
66%|######5 | 131/200 [03:43<01:59, 1.74s/it]
66%|######6 | 132/200 [03:45<01:57, 1.73s/it]
66%|######6 | 133/200 [03:46<01:56, 1.73s/it]
67%|######7 | 134/200 [03:48<01:54, 1.73s/it]
68%|######7 | 135/200 [03:50<01:52, 1.73s/it]
68%|######8 | 136/200 [03:51<01:50, 1.73s/it]
68%|######8 | 137/200 [03:53<01:48, 1.73s/it]
69%|######9 | 138/200 [03:55<01:47, 1.73s/it]
70%|######9 | 139/200 [03:57<01:45, 1.73s/it]
70%|####### | 140/200 [03:58<01:43, 1.73s/it]
70%|####### | 141/200 [04:00<01:42, 1.73s/it]
71%|#######1 | 142/200 [04:02<01:40, 1.73s/it]
72%|#######1 | 143/200 [04:04<01:38, 1.73s/it]
72%|#######2 | 144/200 [04:05<01:36, 1.73s/it]
72%|#######2 | 145/200 [04:07<01:35, 1.73s/it]
73%|#######3 | 146/200 [04:09<01:33, 1.73s/it]
74%|#######3 | 147/200 [04:11<01:31, 1.73s/it]
74%|#######4 | 148/200 [04:12<01:29, 1.73s/it]
74%|#######4 | 149/200 [04:14<01:28, 1.73s/it]
75%|#######5 | 150/200 [04:16<01:26, 1.73s/it]
76%|#######5 | 151/200 [04:17<01:24, 1.73s/it]
76%|#######6 | 152/200 [04:19<01:23, 1.73s/it]
76%|#######6 | 153/200 [04:21<01:21, 1.73s/it]
77%|#######7 | 154/200 [04:23<01:19, 1.73s/it]
78%|#######7 | 155/200 [04:24<01:18, 1.73s/it]
78%|#######8 | 156/200 [04:26<01:16, 1.75s/it]
78%|#######8 | 157/200 [04:28<01:15, 1.75s/it]
79%|#######9 | 158/200 [04:30<01:13, 1.75s/it]
80%|#######9 | 159/200 [04:31<01:11, 1.75s/it]
80%|######## | 160/200 [04:33<01:09, 1.74s/it]
80%|######## | 161/200 [04:35<01:07, 1.74s/it]
81%|########1 | 162/200 [04:37<01:06, 1.74s/it]
82%|########1 | 163/200 [04:38<01:04, 1.74s/it]
82%|########2 | 164/200 [04:40<01:02, 1.74s/it]
82%|########2 | 165/200 [04:42<01:00, 1.74s/it]
83%|########2 | 166/200 [04:44<00:59, 1.74s/it]
84%|########3 | 167/200 [04:45<00:57, 1.74s/it]
84%|########4 | 168/200 [04:47<00:55, 1.74s/it]
84%|########4 | 169/200 [04:49<00:53, 1.74s/it]
85%|########5 | 170/200 [04:51<00:52, 1.74s/it]
86%|########5 | 171/200 [04:52<00:50, 1.74s/it]
86%|########6 | 172/200 [04:54<00:48, 1.73s/it]
86%|########6 | 173/200 [04:56<00:46, 1.73s/it]
87%|########7 | 174/200 [04:57<00:44, 1.70s/it]
88%|########7 | 175/200 [04:59<00:42, 1.69s/it]
88%|########8 | 176/200 [05:01<00:40, 1.67s/it]
88%|########8 | 177/200 [05:02<00:38, 1.69s/it]
89%|########9 | 178/200 [05:04<00:36, 1.67s/it]
90%|########9 | 179/200 [05:06<00:34, 1.67s/it]
90%|######### | 180/200 [05:07<00:33, 1.66s/it]
90%|######### | 181/200 [05:09<00:31, 1.68s/it]
91%|#########1| 182/200 [05:11<00:30, 1.67s/it]
92%|#########1| 183/200 [05:12<00:28, 1.67s/it]
92%|#########2| 184/200 [05:14<00:27, 1.69s/it]
92%|#########2| 185/200 [05:16<00:25, 1.71s/it]
93%|#########3| 186/200 [05:18<00:23, 1.70s/it]
94%|#########3| 187/200 [05:19<00:22, 1.72s/it]
94%|#########3| 188/200 [05:21<00:20, 1.73s/it]
94%|#########4| 189/200 [05:23<00:19, 1.74s/it]
95%|#########5| 190/200 [05:25<00:17, 1.75s/it]
96%|#########5| 191/200 [05:26<00:15, 1.76s/it]
96%|#########6| 192/200 [05:28<00:14, 1.77s/it]
96%|#########6| 193/200 [05:30<00:12, 1.77s/it]
97%|#########7| 194/200 [05:32<00:10, 1.78s/it]
98%|#########7| 195/200 [05:33<00:08, 1.78s/it]
98%|#########8| 196/200 [05:35<00:07, 1.78s/it]
98%|#########8| 197/200 [05:37<00:05, 1.78s/it]
99%|#########9| 198/200 [05:39<00:03, 1.78s/it]
100%|#########9| 199/200 [05:41<00:01, 1.78s/it]
100%|##########| 200/200 [05:42<00:00, 1.78s/it]
100%|##########| 200/200 [05:42<00:00, 1.71s/it]
- final iteration number: 200
- final log10 cost value: 6.0
- converged: False
Done.
Execution time: 342.864671215415 seconds
----------------------------------------
/volatile/Chaithya/actions-runner/_work/pysap/pysap/examples/pysap-mri/non_cartesian_reconstruction.py:105: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
plt.show()
Total running time of the script: ( 6 minutes 33.842 seconds)
