Note
Go to the end to download the full example code. or to run this example in your browser via Binder
Learning sampling pattern with decimation#
An example using PyTorch to showcase learning k-space sampling patterns with decimation.
This example showcases the auto-differentiation capabilities of the NUFFT operator with respect to the k-space trajectory in MRI-nufft.
Hereafter we learn the k-space sample locations \(\mathbf{K}\) using the following cost function:
\[\mathbf{\hat{K}} = arg \min_{\mathbf{K}} || \mathcal{F}_\mathbf{K}^* D_\mathbf{K} \mathcal{F}_\mathbf{K} \mathbf{x} - \mathbf{x} ||_2^2\]
where \(\mathcal{F}_\mathbf{K}\) is the forward NUFFT operator, \(D_\mathbf{K}\) is the density compensator for trajectory \(\mathbf{K}\), and \(\mathbf{x}\) is the MR image which is also the target image to be reconstructed.
Additionally, in order to converge faster, we also learn the trajectory in a multi-resolution fashion. This is done by first optimizing x8 times decimated trajectory locations, called control points. After a fixed number of iterations (5 in this example), these control points are upscaled by a factor of 2. Note that the NUFFT operator always holds linearly interpolated version of the control points as k-space sampling trajectory.
Note
This example can run on a binder instance as it is purely CPU based backend (finufft), and is restricted to a 2D single coil toy case.
Warning
This example only showcases the auto-differentiation capabilities, the learned sampling pattern is not scanner compliant as the gradients required to implement it violate the hardware constraints. In practice, a projection \(\Pi_\mathcal{Q}(\mathbf{K})\) onto the scanner constraints set \(\mathcal{Q}\) is recommended (see [Cha+16]). This is implemented in the proprietary SPARKLING package [Cha+22]. Users are encouraged to contact the authors if they want to use it.
import time
import brainweb_dl as bwdl
import joblib
import matplotlib.pyplot as plt
import numpy as np
import tempfile as tmp
import torch
from PIL import Image, ImageSequence
from tqdm import tqdm
from mrinufft import get_operator
from mrinufft.trajectories import initialize_2D_radial
/volatile/github-ci-mind-inria/action-runner/_work/_tool/Python/3.10.14/x64/lib/python3.10/site-packages/cupy/_environment.py:487: UserWarning:
--------------------------------------------------------------------------------
CuPy may not function correctly because multiple CuPy packages are installed
in your environment:
cupy-cuda11x, cupy-cuda12x
Follow these steps to resolve this issue:
1. For all packages listed above, run the following command to remove all
existing CuPy installations:
$ pip uninstall <package_name>
If you previously installed CuPy via conda, also run the following:
$ conda uninstall cupy
2. Install the appropriate CuPy package.
Refer to the Installation Guide for detailed instructions.
https://docs.cupy.dev/en/stable/install.html
--------------------------------------------------------------------------------
warnings.warn(f'''
Utils#
Model class#
Note
While we are only learning the NUFFT operator, we still need the gradient wrt_data=True to have all the gradients computed correctly. See [GRC23] for more details.
class Model(torch.nn.Module):
def __init__(
self,
inital_trajectory,
img_size=(256, 256),
start_decim=8,
interpolation_mode="linear",
):
super(Model, self).__init__()
self.control = torch.nn.Parameter(
data=torch.Tensor(inital_trajectory[:, ::start_decim]),
requires_grad=True,
)
self.current_decim = start_decim
self.interpolation_mode = interpolation_mode
sample_points = inital_trajectory.reshape(-1, inital_trajectory.shape[-1])
self.operator = get_operator("finufft", wrt_data=True, wrt_traj=True)(
sample_points,
shape=img_size,
squeeze_dims=False,
)
self.img_size = img_size
def _interpolate(self, traj, factor=2):
"""Torch interpolate function to upsample the trajectory"""
return torch.nn.functional.interpolate(
traj.moveaxis(1, -1),
scale_factor=factor,
mode=self.interpolation_mode,
align_corners=True,
).moveaxis(-1, 1)
def get_trajectory(self):
"""Function to get trajectory, which is interpolated version of control points."""
traj = self.control.clone()
for i in range(np.log2(self.current_decim).astype(int)):
traj = self._interpolate(traj)
return traj.reshape(-1, traj.shape[-1])
def upscale(self, factor=2):
"""Upscaling the model.
In this step, the number of control points are doubled and interpolated.
"""
self.control = torch.nn.Parameter(
data=self._interpolate(self.control),
requires_grad=True,
)
self.current_decim /= factor
def forward(self, x):
traj = self.get_trajectory()
self.operator.samples = traj
# Simulate the acquisition process
kspace = self.operator.op(x)
adjoint = self.operator.adj_op(kspace).abs()
return adjoint / torch.mean(adjoint)
State plotting#
def plot_state(axs, image, traj, recon, control_points=None, loss=None, save_name=None):
axs = axs.flatten()
# Upper left reference image
axs[0].imshow(np.abs(image[0]), cmap="gray")
axs[0].axis("off")
axs[0].set_title("MR Image")
# Upper right trajectory
axs[1].scatter(*traj.T, s=0.5)
if control_points is not None:
axs[1].scatter(*control_points.T, s=1, color="r")
axs[1].legend(
["Trajectory", "Control points"], loc="right", bbox_to_anchor=(2, 0.6)
)
axs[1].grid(True)
axs[1].set_title("Trajectory")
axs[1].set_xlim(-0.5, 0.5)
axs[1].set_ylim(-0.5, 0.5)
axs[1].set_aspect("equal")
# Down left reconstructed image
axs[2].imshow(np.abs(recon[0][0].detach().cpu().numpy()), cmap="gray")
axs[2].axis("off")
axs[2].set_title("Reconstruction")
# Down right loss evolution
if loss is not None:
axs[3].plot(loss)
axs[3].set_ylim(0, None)
axs[3].grid("on")
axs[3].set_title("Loss")
plt.subplots_adjust(hspace=0.3)
# Save & close
if save_name is not None:
plt.savefig(save_name, bbox_inches="tight")
plt.close()
else:
plt.show()
Optimizer upscaling#
def upsample_optimizer(optimizer, new_optimizer, factor=2):
"""Upsample the optimizer."""
for old_group, new_group in zip(optimizer.param_groups, new_optimizer.param_groups):
for old_param, new_param in zip(old_group["params"], new_group["params"]):
# Interpolate optimizer states
if old_param in optimizer.state:
for key in optimizer.state[old_param].keys():
if isinstance(optimizer.state[old_param][key], torch.Tensor):
old_state = optimizer.state[old_param][key]
if old_state.ndim == 0:
new_state = old_state
else:
new_state = torch.nn.functional.interpolate(
old_state.moveaxis(1, -1),
scale_factor=factor,
mode="linear",
).moveaxis(-1, 1)
new_optimizer.state[new_param][key] = new_state
else:
new_optimizer.state[new_param][key] = optimizer.state[
old_param
][key]
return new_optimizer
Data preparation#
A single image to train the model over. Note that in practice we would use a whole dataset instead (e.g. fastMRI).
volume = np.flip(bwdl.get_mri(4, "T1"), axis=(0, 1, 2))
image = torch.from_numpy(volume[-80, ...].astype(np.float32))[None]
image = image / torch.mean(image)
A basic radial trajectory with an acceleration factor of 8.
AF = 8
initial_traj = initialize_2D_radial(image.shape[1] // AF, image.shape[2]).astype(
np.float32
)
Trajectory learning#
Initialisation#
model = Model(initial_traj, img_size=image.shape[1:])
model = model.eval()
/volatile/github-ci-mind-inria/action-runner/_work/_tool/Python/3.10.14/x64/lib/python3.10/site-packages/mrinufft/_utils.py:94: UserWarning: Samples will be rescaled to [-pi, pi), assuming they were in [-0.5, 0.5)
warnings.warn(
The image obtained before learning the sampling pattern is highly degraded because of the acceleration factor and simplicity of the trajectory.
initial_recons = model(image)
fig, axs = plt.subplots(1, 3, figsize=(9, 3))
plot_state(axs, image, initial_traj, initial_recons)
/volatile/github-ci-mind-inria/action-runner/_work/_tool/Python/3.10.14/x64/lib/python3.10/site-packages/mrinufft/_utils.py:94: UserWarning: Samples will be rescaled to [-pi, pi), assuming they were in [-0.5, 0.5)
warnings.warn(
/volatile/github-ci-mind-inria/action-runner/_work/_tool/Python/3.10.14/x64/lib/python3.10/site-packages/finufft/_interfaces.py:329: UserWarning: Argument `data` does not satisfy the following requirement: C. Copying array (this may reduce performance)
warnings.warn(f"Argument `{name}` does not satisfy the following requirement: {prop}. Copying array (this may reduce performance)")
Training loop#
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
model.train()
losses = []
image_files = []
while model.current_decim >= 1:
with tqdm(range(30), unit="steps") as tqdms:
for i in tqdms:
out = model(image)
loss = torch.nn.functional.mse_loss(out, image[None, None])
numpy_loss = (loss.detach().cpu().numpy(),)
tqdms.set_postfix({"loss": numpy_loss})
losses.append(numpy_loss)
optimizer.zero_grad()
loss.backward()
optimizer.step()
with torch.no_grad():
# Clamp the value of trajectory between [-0.5, 0.5]
for param in model.parameters():
param.clamp_(-0.5, 0.5)
# Generate images for gif
filename = f"{tmp.NamedTemporaryFile().name}.png"
plt.clf()
fig, axs = plt.subplots(2, 2, figsize=(10, 10), num=1)
plot_state(
axs,
image,
model.get_trajectory().detach().cpu().numpy(),
out,
model.control.detach().cpu().numpy(),
losses,
save_name=filename,
)
image_files.append(filename)
if model.current_decim == 1:
break
else:
model.upscale()
optimizer = upsample_optimizer(
optimizer, torch.optim.Adam(model.parameters(), lr=1e-3)
)
0%| | 0/30 [00:00<?, ?steps/s]/volatile/github-ci-mind-inria/action-runner/_work/_tool/Python/3.10.14/x64/lib/python3.10/site-packages/mrinufft/_utils.py:94: UserWarning: Samples will be rescaled to [-pi, pi), assuming they were in [-0.5, 0.5)
warnings.warn(
/volatile/github-ci-mind-inria/action-runner/_work/_tool/Python/3.10.14/x64/lib/python3.10/site-packages/finufft/_interfaces.py:329: UserWarning: Argument `data` does not satisfy the following requirement: C. Copying array (this may reduce performance)
warnings.warn(f"Argument `{name}` does not satisfy the following requirement: {prop}. Copying array (this may reduce performance)")
/volatile/github-ci-mind-inria/action-runner/_work/mri-nufft/mri-nufft/examples/example_learn_samples_multires.py:259: UserWarning: Using a target size (torch.Size([1, 1, 1, 256, 256])) that is different to the input size (torch.Size([1, 1, 256, 256])). This will likely lead to incorrect results due to broadcasting. Please ensure they have the same size.
loss = torch.nn.functional.mse_loss(out, image[None, None])
0%| | 0/30 [00:00<?, ?steps/s, loss=(array(0.74174887, dtype=float32),)]/volatile/github-ci-mind-inria/action-runner/_work/_tool/Python/3.10.14/x64/lib/python3.10/site-packages/mrinufft/operators/autodiff.py:98: UserWarning: Casting complex values to real discards the imaginary part (Triggered internally at ../aten/src/ATen/native/Copy.cpp:305.)
grad_traj = torch.transpose(torch.sum(grad_traj, dim=1), 0, 1).to(
3%|▎ | 1/30 [00:00<00:14, 1.96steps/s, loss=(array(0.74174887, dtype=float32),)]
3%|▎ | 1/30 [00:00<00:14, 1.96steps/s, loss=(array(0.6624898, dtype=float32),)]
7%|▋ | 2/30 [00:00<00:13, 2.04steps/s, loss=(array(0.6624898, dtype=float32),)]
7%|▋ | 2/30 [00:00<00:13, 2.04steps/s, loss=(array(0.5745859, dtype=float32),)]
10%|█ | 3/30 [00:01<00:15, 1.79steps/s, loss=(array(0.5745859, dtype=float32),)]
10%|█ | 3/30 [00:01<00:15, 1.79steps/s, loss=(array(0.48269865, dtype=float32),)]
13%|█▎ | 4/30 [00:02<00:13, 1.90steps/s, loss=(array(0.48269865, dtype=float32),)]
13%|█▎ | 4/30 [00:02<00:13, 1.90steps/s, loss=(array(0.39541537, dtype=float32),)]
17%|█▋ | 5/30 [00:02<00:12, 1.96steps/s, loss=(array(0.39541537, dtype=float32),)]
17%|█▋ | 5/30 [00:02<00:12, 1.96steps/s, loss=(array(0.32777992, dtype=float32),)]
20%|██ | 6/30 [00:03<00:11, 2.02steps/s, loss=(array(0.32777992, dtype=float32),)]
20%|██ | 6/30 [00:03<00:11, 2.02steps/s, loss=(array(0.2930662, dtype=float32),)]
23%|██▎ | 7/30 [00:03<00:11, 2.07steps/s, loss=(array(0.2930662, dtype=float32),)]
23%|██▎ | 7/30 [00:03<00:11, 2.07steps/s, loss=(array(0.2890394, dtype=float32),)]
27%|██▋ | 8/30 [00:03<00:10, 2.11steps/s, loss=(array(0.2890394, dtype=float32),)]
27%|██▋ | 8/30 [00:03<00:10, 2.11steps/s, loss=(array(0.30443317, dtype=float32),)]
30%|███ | 9/30 [00:04<00:09, 2.13steps/s, loss=(array(0.30443317, dtype=float32),)]
30%|███ | 9/30 [00:04<00:09, 2.13steps/s, loss=(array(0.3226806, dtype=float32),)]
33%|███▎ | 10/30 [00:04<00:09, 2.15steps/s, loss=(array(0.3226806, dtype=float32),)]
33%|███▎ | 10/30 [00:04<00:09, 2.15steps/s, loss=(array(0.33095258, dtype=float32),)]
37%|███▋ | 11/30 [00:05<00:08, 2.15steps/s, loss=(array(0.33095258, dtype=float32),)]
37%|███▋ | 11/30 [00:05<00:08, 2.15steps/s, loss=(array(0.32629323, dtype=float32),)]
40%|████ | 12/30 [00:06<00:09, 1.86steps/s, loss=(array(0.32629323, dtype=float32),)]
40%|████ | 12/30 [00:06<00:09, 1.86steps/s, loss=(array(0.3123067, dtype=float32),)]
43%|████▎ | 13/30 [00:06<00:08, 1.96steps/s, loss=(array(0.3123067, dtype=float32),)]
43%|████▎ | 13/30 [00:06<00:08, 1.96steps/s, loss=(array(0.29476863, dtype=float32),)]
47%|████▋ | 14/30 [00:06<00:07, 2.03steps/s, loss=(array(0.29476863, dtype=float32),)]
47%|████▋ | 14/30 [00:06<00:07, 2.03steps/s, loss=(array(0.27900735, dtype=float32),)]
50%|█████ | 15/30 [00:07<00:07, 2.09steps/s, loss=(array(0.27900735, dtype=float32),)]
50%|█████ | 15/30 [00:07<00:07, 2.09steps/s, loss=(array(0.268334, dtype=float32),)]
53%|█████▎ | 16/30 [00:07<00:06, 2.13steps/s, loss=(array(0.268334, dtype=float32),)]
53%|█████▎ | 16/30 [00:07<00:06, 2.13steps/s, loss=(array(0.2633488, dtype=float32),)]
57%|█████▋ | 17/30 [00:08<00:06, 2.15steps/s, loss=(array(0.2633488, dtype=float32),)]
57%|█████▋ | 17/30 [00:08<00:06, 2.15steps/s, loss=(array(0.2626493, dtype=float32),)]
60%|██████ | 18/30 [00:08<00:05, 2.16steps/s, loss=(array(0.2626493, dtype=float32),)]
60%|██████ | 18/30 [00:08<00:05, 2.16steps/s, loss=(array(0.2641289, dtype=float32),)]
63%|██████▎ | 19/30 [00:09<00:05, 2.18steps/s, loss=(array(0.2641289, dtype=float32),)]
63%|██████▎ | 19/30 [00:09<00:05, 2.18steps/s, loss=(array(0.26640078, dtype=float32),)]
67%|██████▋ | 20/30 [00:09<00:04, 2.19steps/s, loss=(array(0.26640078, dtype=float32),)]
67%|██████▋ | 20/30 [00:09<00:04, 2.19steps/s, loss=(array(0.26791817, dtype=float32),)]
70%|███████ | 21/30 [00:10<00:04, 2.17steps/s, loss=(array(0.26791817, dtype=float32),)]
70%|███████ | 21/30 [00:10<00:04, 2.17steps/s, loss=(array(0.2674212, dtype=float32),)]
73%|███████▎ | 22/30 [00:10<00:04, 1.88steps/s, loss=(array(0.2674212, dtype=float32),)]
73%|███████▎ | 22/30 [00:10<00:04, 1.88steps/s, loss=(array(0.26476443, dtype=float32),)]
77%|███████▋ | 23/30 [00:11<00:03, 1.95steps/s, loss=(array(0.26476443, dtype=float32),)]
77%|███████▋ | 23/30 [00:11<00:03, 1.95steps/s, loss=(array(0.26061833, dtype=float32),)]
80%|████████ | 24/30 [00:11<00:02, 2.01steps/s, loss=(array(0.26061833, dtype=float32),)]
80%|████████ | 24/30 [00:11<00:02, 2.01steps/s, loss=(array(0.2559579, dtype=float32),)]
83%|████████▎ | 25/30 [00:12<00:02, 2.06steps/s, loss=(array(0.2559579, dtype=float32),)]
83%|████████▎ | 25/30 [00:12<00:02, 2.06steps/s, loss=(array(0.2515955, dtype=float32),)]
87%|████████▋ | 26/30 [00:12<00:01, 2.09steps/s, loss=(array(0.2515955, dtype=float32),)]
87%|████████▋ | 26/30 [00:12<00:01, 2.09steps/s, loss=(array(0.24800554, dtype=float32),)]
90%|█████████ | 27/30 [00:13<00:01, 2.10steps/s, loss=(array(0.24800554, dtype=float32),)]
90%|█████████ | 27/30 [00:13<00:01, 2.10steps/s, loss=(array(0.2453642, dtype=float32),)]
93%|█████████▎| 28/30 [00:13<00:00, 2.11steps/s, loss=(array(0.2453642, dtype=float32),)]
93%|█████████▎| 28/30 [00:13<00:00, 2.11steps/s, loss=(array(0.24350482, dtype=float32),)]
97%|█████████▋| 29/30 [00:14<00:00, 2.07steps/s, loss=(array(0.24350482, dtype=float32),)]
97%|█████████▋| 29/30 [00:14<00:00, 2.07steps/s, loss=(array(0.24218965, dtype=float32),)]
100%|██████████| 30/30 [00:14<00:00, 2.08steps/s, loss=(array(0.24218965, dtype=float32),)]
100%|██████████| 30/30 [00:14<00:00, 2.06steps/s, loss=(array(0.24218965, dtype=float32),)]
0%| | 0/30 [00:00<?, ?steps/s]
0%| | 0/30 [00:00<?, ?steps/s, loss=(array(0.24114653, dtype=float32),)]
3%|▎ | 1/30 [00:00<00:19, 1.46steps/s, loss=(array(0.24114653, dtype=float32),)]
3%|▎ | 1/30 [00:00<00:19, 1.46steps/s, loss=(array(0.24032444, dtype=float32),)]
7%|▋ | 2/30 [00:01<00:15, 1.80steps/s, loss=(array(0.24032444, dtype=float32),)]
7%|▋ | 2/30 [00:01<00:15, 1.80steps/s, loss=(array(0.23924257, dtype=float32),)]
10%|█ | 3/30 [00:01<00:14, 1.92steps/s, loss=(array(0.23924257, dtype=float32),)]
10%|█ | 3/30 [00:01<00:14, 1.92steps/s, loss=(array(0.23780563, dtype=float32),)]
13%|█▎ | 4/30 [00:02<00:12, 2.00steps/s, loss=(array(0.23780563, dtype=float32),)]
13%|█▎ | 4/30 [00:02<00:12, 2.00steps/s, loss=(array(0.23599169, dtype=float32),)]
17%|█▋ | 5/30 [00:02<00:12, 2.04steps/s, loss=(array(0.23599169, dtype=float32),)]
17%|█▋ | 5/30 [00:02<00:12, 2.04steps/s, loss=(array(0.23385708, dtype=float32),)]
20%|██ | 6/30 [00:03<00:11, 2.06steps/s, loss=(array(0.23385708, dtype=float32),)]
20%|██ | 6/30 [00:03<00:11, 2.06steps/s, loss=(array(0.23151523, dtype=float32),)]
23%|██▎ | 7/30 [00:03<00:11, 2.07steps/s, loss=(array(0.23151523, dtype=float32),)]
23%|██▎ | 7/30 [00:03<00:11, 2.07steps/s, loss=(array(0.22911337, dtype=float32),)]
27%|██▋ | 8/30 [00:03<00:10, 2.09steps/s, loss=(array(0.22911337, dtype=float32),)]
27%|██▋ | 8/30 [00:04<00:10, 2.09steps/s, loss=(array(0.22678916, dtype=float32),)]
30%|███ | 9/30 [00:04<00:09, 2.11steps/s, loss=(array(0.22678916, dtype=float32),)]
30%|███ | 9/30 [00:04<00:09, 2.11steps/s, loss=(array(0.22464427, dtype=float32),)]
33%|███▎ | 10/30 [00:04<00:09, 2.13steps/s, loss=(array(0.22464427, dtype=float32),)]
33%|███▎ | 10/30 [00:04<00:09, 2.13steps/s, loss=(array(0.22273009, dtype=float32),)]
37%|███▋ | 11/30 [00:05<00:10, 1.88steps/s, loss=(array(0.22273009, dtype=float32),)]
37%|███▋ | 11/30 [00:05<00:10, 1.88steps/s, loss=(array(0.22105053, dtype=float32),)]
40%|████ | 12/30 [00:06<00:09, 1.96steps/s, loss=(array(0.22105053, dtype=float32),)]
40%|████ | 12/30 [00:06<00:09, 1.96steps/s, loss=(array(0.21956299, dtype=float32),)]
43%|████▎ | 13/30 [00:06<00:08, 2.01steps/s, loss=(array(0.21956299, dtype=float32),)]
43%|████▎ | 13/30 [00:06<00:08, 2.01steps/s, loss=(array(0.21819022, dtype=float32),)]
47%|████▋ | 14/30 [00:06<00:07, 2.04steps/s, loss=(array(0.21819022, dtype=float32),)]
47%|████▋ | 14/30 [00:07<00:07, 2.04steps/s, loss=(array(0.21684906, dtype=float32),)]
50%|█████ | 15/30 [00:07<00:07, 2.04steps/s, loss=(array(0.21684906, dtype=float32),)]
50%|█████ | 15/30 [00:07<00:07, 2.04steps/s, loss=(array(0.21548082, dtype=float32),)]
53%|█████▎ | 16/30 [00:07<00:06, 2.06steps/s, loss=(array(0.21548082, dtype=float32),)]
53%|█████▎ | 16/30 [00:07<00:06, 2.06steps/s, loss=(array(0.21406442, dtype=float32),)]
57%|█████▋ | 17/30 [00:08<00:06, 2.08steps/s, loss=(array(0.21406442, dtype=float32),)]
57%|█████▋ | 17/30 [00:08<00:06, 2.08steps/s, loss=(array(0.21261415, dtype=float32),)]
60%|██████ | 18/30 [00:08<00:05, 2.03steps/s, loss=(array(0.21261415, dtype=float32),)]
60%|██████ | 18/30 [00:08<00:05, 2.03steps/s, loss=(array(0.2111668, dtype=float32),)]
63%|██████▎ | 19/30 [00:09<00:05, 2.05steps/s, loss=(array(0.2111668, dtype=float32),)]
63%|██████▎ | 19/30 [00:09<00:05, 2.05steps/s, loss=(array(0.20976183, dtype=float32),)]
67%|██████▋ | 20/30 [00:10<00:05, 1.82steps/s, loss=(array(0.20976183, dtype=float32),)]
67%|██████▋ | 20/30 [00:10<00:05, 1.82steps/s, loss=(array(0.20842579, dtype=float32),)]
70%|███████ | 21/30 [00:10<00:04, 1.91steps/s, loss=(array(0.20842579, dtype=float32),)]
70%|███████ | 21/30 [00:10<00:04, 1.91steps/s, loss=(array(0.20716691, dtype=float32),)]
73%|███████▎ | 22/30 [00:11<00:04, 1.98steps/s, loss=(array(0.20716691, dtype=float32),)]
73%|███████▎ | 22/30 [00:11<00:04, 1.98steps/s, loss=(array(0.20597577, dtype=float32),)]
77%|███████▋ | 23/30 [00:11<00:03, 2.03steps/s, loss=(array(0.20597577, dtype=float32),)]
77%|███████▋ | 23/30 [00:11<00:03, 2.03steps/s, loss=(array(0.20482925, dtype=float32),)]
80%|████████ | 24/30 [00:11<00:02, 2.08steps/s, loss=(array(0.20482925, dtype=float32),)]
80%|████████ | 24/30 [00:11<00:02, 2.08steps/s, loss=(array(0.20370199, dtype=float32),)]
83%|████████▎ | 25/30 [00:12<00:02, 2.12steps/s, loss=(array(0.20370199, dtype=float32),)]
83%|████████▎ | 25/30 [00:12<00:02, 2.12steps/s, loss=(array(0.20257255, dtype=float32),)]
87%|████████▋ | 26/30 [00:12<00:01, 2.15steps/s, loss=(array(0.20257255, dtype=float32),)]
87%|████████▋ | 26/30 [00:12<00:01, 2.15steps/s, loss=(array(0.20143044, dtype=float32),)]
90%|█████████ | 27/30 [00:13<00:01, 2.17steps/s, loss=(array(0.20143044, dtype=float32),)]
90%|█████████ | 27/30 [00:13<00:01, 2.17steps/s, loss=(array(0.20027772, dtype=float32),)]
93%|█████████▎| 28/30 [00:13<00:00, 2.18steps/s, loss=(array(0.20027772, dtype=float32),)]
93%|█████████▎| 28/30 [00:13<00:00, 2.18steps/s, loss=(array(0.19912508, dtype=float32),)]
97%|█████████▋| 29/30 [00:14<00:00, 2.17steps/s, loss=(array(0.19912508, dtype=float32),)]
97%|█████████▋| 29/30 [00:14<00:00, 2.17steps/s, loss=(array(0.19798627, dtype=float32),)]
100%|██████████| 30/30 [00:14<00:00, 1.91steps/s, loss=(array(0.19798627, dtype=float32),)]
100%|██████████| 30/30 [00:14<00:00, 2.01steps/s, loss=(array(0.19798627, dtype=float32),)]
0%| | 0/30 [00:00<?, ?steps/s]
0%| | 0/30 [00:00<?, ?steps/s, loss=(array(0.19687007, dtype=float32),)]
3%|▎ | 1/30 [00:00<00:13, 2.20steps/s, loss=(array(0.19687007, dtype=float32),)]
3%|▎ | 1/30 [00:00<00:13, 2.20steps/s, loss=(array(0.1958951, dtype=float32),)]
7%|▋ | 2/30 [00:00<00:12, 2.17steps/s, loss=(array(0.1958951, dtype=float32),)]
7%|▋ | 2/30 [00:00<00:12, 2.17steps/s, loss=(array(0.19494593, dtype=float32),)]
10%|█ | 3/30 [00:01<00:12, 2.17steps/s, loss=(array(0.19494593, dtype=float32),)]
10%|█ | 3/30 [00:01<00:12, 2.17steps/s, loss=(array(0.194014, dtype=float32),)]
13%|█▎ | 4/30 [00:01<00:11, 2.17steps/s, loss=(array(0.194014, dtype=float32),)]
13%|█▎ | 4/30 [00:01<00:11, 2.17steps/s, loss=(array(0.19309294, dtype=float32),)]
17%|█▋ | 5/30 [00:02<00:11, 2.17steps/s, loss=(array(0.19309294, dtype=float32),)]
17%|█▋ | 5/30 [00:02<00:11, 2.17steps/s, loss=(array(0.19218084, dtype=float32),)]
20%|██ | 6/30 [00:02<00:11, 2.17steps/s, loss=(array(0.19218084, dtype=float32),)]
20%|██ | 6/30 [00:02<00:11, 2.17steps/s, loss=(array(0.19128107, dtype=float32),)]
23%|██▎ | 7/30 [00:03<00:10, 2.17steps/s, loss=(array(0.19128107, dtype=float32),)]
23%|██▎ | 7/30 [00:03<00:10, 2.17steps/s, loss=(array(0.19040057, dtype=float32),)]
27%|██▋ | 8/30 [00:03<00:10, 2.17steps/s, loss=(array(0.19040057, dtype=float32),)]
27%|██▋ | 8/30 [00:03<00:10, 2.17steps/s, loss=(array(0.1895481, dtype=float32),)]
30%|███ | 9/30 [00:04<00:09, 2.17steps/s, loss=(array(0.1895481, dtype=float32),)]
30%|███ | 9/30 [00:04<00:09, 2.17steps/s, loss=(array(0.18873104, dtype=float32),)]
33%|███▎ | 10/30 [00:04<00:10, 1.89steps/s, loss=(array(0.18873104, dtype=float32),)]
33%|███▎ | 10/30 [00:04<00:10, 1.89steps/s, loss=(array(0.18795341, dtype=float32),)]
37%|███▋ | 11/30 [00:05<00:09, 1.96steps/s, loss=(array(0.18795341, dtype=float32),)]
37%|███▋ | 11/30 [00:05<00:09, 1.96steps/s, loss=(array(0.18721464, dtype=float32),)]
40%|████ | 12/30 [00:05<00:08, 2.02steps/s, loss=(array(0.18721464, dtype=float32),)]
40%|████ | 12/30 [00:05<00:08, 2.02steps/s, loss=(array(0.18651058, dtype=float32),)]
43%|████▎ | 13/30 [00:06<00:08, 2.06steps/s, loss=(array(0.18651058, dtype=float32),)]
43%|████▎ | 13/30 [00:06<00:08, 2.06steps/s, loss=(array(0.18583554, dtype=float32),)]
47%|████▋ | 14/30 [00:06<00:07, 2.10steps/s, loss=(array(0.18583554, dtype=float32),)]
47%|████▋ | 14/30 [00:06<00:07, 2.10steps/s, loss=(array(0.18518463, dtype=float32),)]
50%|█████ | 15/30 [00:07<00:07, 2.12steps/s, loss=(array(0.18518463, dtype=float32),)]
50%|█████ | 15/30 [00:07<00:07, 2.12steps/s, loss=(array(0.18455419, dtype=float32),)]
53%|█████▎ | 16/30 [00:07<00:06, 2.15steps/s, loss=(array(0.18455419, dtype=float32),)]
53%|█████▎ | 16/30 [00:07<00:06, 2.15steps/s, loss=(array(0.18394282, dtype=float32),)]
57%|█████▋ | 17/30 [00:08<00:06, 2.16steps/s, loss=(array(0.18394282, dtype=float32),)]
57%|█████▋ | 17/30 [00:08<00:06, 2.16steps/s, loss=(array(0.18335065, dtype=float32),)]
60%|██████ | 18/30 [00:08<00:05, 2.16steps/s, loss=(array(0.18335065, dtype=float32),)]
60%|██████ | 18/30 [00:08<00:05, 2.16steps/s, loss=(array(0.18277875, dtype=float32),)]
63%|██████▎ | 19/30 [00:09<00:05, 1.86steps/s, loss=(array(0.18277875, dtype=float32),)]
63%|██████▎ | 19/30 [00:09<00:05, 1.86steps/s, loss=(array(0.18222785, dtype=float32),)]
67%|██████▋ | 20/30 [00:09<00:05, 1.94steps/s, loss=(array(0.18222785, dtype=float32),)]
67%|██████▋ | 20/30 [00:09<00:05, 1.94steps/s, loss=(array(0.18169734, dtype=float32),)]
70%|███████ | 21/30 [00:10<00:04, 2.00steps/s, loss=(array(0.18169734, dtype=float32),)]
70%|███████ | 21/30 [00:10<00:04, 2.00steps/s, loss=(array(0.18118508, dtype=float32),)]
73%|███████▎ | 22/30 [00:10<00:03, 2.04steps/s, loss=(array(0.18118508, dtype=float32),)]
73%|███████▎ | 22/30 [00:10<00:03, 2.04steps/s, loss=(array(0.18068832, dtype=float32),)]
77%|███████▋ | 23/30 [00:11<00:03, 2.06steps/s, loss=(array(0.18068832, dtype=float32),)]
77%|███████▋ | 23/30 [00:11<00:03, 2.06steps/s, loss=(array(0.1802041, dtype=float32),)]
80%|████████ | 24/30 [00:11<00:02, 2.08steps/s, loss=(array(0.1802041, dtype=float32),)]
80%|████████ | 24/30 [00:11<00:02, 2.08steps/s, loss=(array(0.17973003, dtype=float32),)]
83%|████████▎ | 25/30 [00:12<00:02, 2.09steps/s, loss=(array(0.17973003, dtype=float32),)]
83%|████████▎ | 25/30 [00:12<00:02, 2.09steps/s, loss=(array(0.17926493, dtype=float32),)]
87%|████████▋ | 26/30 [00:12<00:01, 2.10steps/s, loss=(array(0.17926493, dtype=float32),)]
87%|████████▋ | 26/30 [00:12<00:01, 2.10steps/s, loss=(array(0.17880903, dtype=float32),)]
90%|█████████ | 27/30 [00:12<00:01, 2.10steps/s, loss=(array(0.17880903, dtype=float32),)]
90%|█████████ | 27/30 [00:12<00:01, 2.10steps/s, loss=(array(0.17836313, dtype=float32),)]
93%|█████████▎| 28/30 [00:13<00:00, 2.12steps/s, loss=(array(0.17836313, dtype=float32),)]
93%|█████████▎| 28/30 [00:13<00:00, 2.12steps/s, loss=(array(0.17792809, dtype=float32),)]
97%|█████████▋| 29/30 [00:14<00:00, 1.83steps/s, loss=(array(0.17792809, dtype=float32),)]
97%|█████████▋| 29/30 [00:14<00:00, 1.83steps/s, loss=(array(0.17750451, dtype=float32),)]
100%|██████████| 30/30 [00:14<00:00, 1.91steps/s, loss=(array(0.17750451, dtype=float32),)]
100%|██████████| 30/30 [00:14<00:00, 2.05steps/s, loss=(array(0.17750451, dtype=float32),)]
0%| | 0/30 [00:00<?, ?steps/s]
0%| | 0/30 [00:00<?, ?steps/s, loss=(array(0.17709179, dtype=float32),)]
3%|▎ | 1/30 [00:00<00:13, 2.14steps/s, loss=(array(0.17709179, dtype=float32),)]
3%|▎ | 1/30 [00:00<00:13, 2.14steps/s, loss=(array(0.17668465, dtype=float32),)]
7%|▋ | 2/30 [00:00<00:13, 2.14steps/s, loss=(array(0.17668465, dtype=float32),)]
7%|▋ | 2/30 [00:00<00:13, 2.14steps/s, loss=(array(0.17627303, dtype=float32),)]
10%|█ | 3/30 [00:01<00:12, 2.13steps/s, loss=(array(0.17627303, dtype=float32),)]
10%|█ | 3/30 [00:01<00:12, 2.13steps/s, loss=(array(0.17585906, dtype=float32),)]
13%|█▎ | 4/30 [00:01<00:12, 2.13steps/s, loss=(array(0.17585906, dtype=float32),)]
13%|█▎ | 4/30 [00:01<00:12, 2.13steps/s, loss=(array(0.17544556, dtype=float32),)]
17%|█▋ | 5/30 [00:02<00:11, 2.13steps/s, loss=(array(0.17544556, dtype=float32),)]
17%|█▋ | 5/30 [00:02<00:11, 2.13steps/s, loss=(array(0.17503598, dtype=float32),)]
20%|██ | 6/30 [00:02<00:11, 2.13steps/s, loss=(array(0.17503598, dtype=float32),)]
20%|██ | 6/30 [00:02<00:11, 2.13steps/s, loss=(array(0.17463395, dtype=float32),)]
23%|██▎ | 7/30 [00:03<00:10, 2.13steps/s, loss=(array(0.17463395, dtype=float32),)]
23%|██▎ | 7/30 [00:03<00:10, 2.13steps/s, loss=(array(0.17424281, dtype=float32),)]
27%|██▋ | 8/30 [00:04<00:12, 1.79steps/s, loss=(array(0.17424281, dtype=float32),)]
27%|██▋ | 8/30 [00:04<00:12, 1.79steps/s, loss=(array(0.1738656, dtype=float32),)]
30%|███ | 9/30 [00:04<00:11, 1.88steps/s, loss=(array(0.1738656, dtype=float32),)]
30%|███ | 9/30 [00:04<00:11, 1.88steps/s, loss=(array(0.17350477, dtype=float32),)]
33%|███▎ | 10/30 [00:04<00:10, 1.94steps/s, loss=(array(0.17350477, dtype=float32),)]
33%|███▎ | 10/30 [00:05<00:10, 1.94steps/s, loss=(array(0.17316222, dtype=float32),)]
37%|███▋ | 11/30 [00:05<00:09, 1.99steps/s, loss=(array(0.17316222, dtype=float32),)]
37%|███▋ | 11/30 [00:05<00:09, 1.99steps/s, loss=(array(0.1728391, dtype=float32),)]
40%|████ | 12/30 [00:05<00:08, 2.03steps/s, loss=(array(0.1728391, dtype=float32),)]
40%|████ | 12/30 [00:05<00:08, 2.03steps/s, loss=(array(0.17253551, dtype=float32),)]
43%|████▎ | 13/30 [00:06<00:08, 2.06steps/s, loss=(array(0.17253551, dtype=float32),)]
43%|████▎ | 13/30 [00:06<00:08, 2.06steps/s, loss=(array(0.17225075, dtype=float32),)]
47%|████▋ | 14/30 [00:06<00:07, 2.08steps/s, loss=(array(0.17225075, dtype=float32),)]
47%|████▋ | 14/30 [00:06<00:07, 2.08steps/s, loss=(array(0.17198272, dtype=float32),)]
50%|█████ | 15/30 [00:07<00:07, 2.09steps/s, loss=(array(0.17198272, dtype=float32),)]
50%|█████ | 15/30 [00:07<00:07, 2.09steps/s, loss=(array(0.17172849, dtype=float32),)]
53%|█████▎ | 16/30 [00:07<00:06, 2.10steps/s, loss=(array(0.17172849, dtype=float32),)]
53%|█████▎ | 16/30 [00:07<00:06, 2.10steps/s, loss=(array(0.17148453, dtype=float32),)]
57%|█████▋ | 17/30 [00:08<00:06, 2.10steps/s, loss=(array(0.17148453, dtype=float32),)]
57%|█████▋ | 17/30 [00:08<00:06, 2.10steps/s, loss=(array(0.17124742, dtype=float32),)]
60%|██████ | 18/30 [00:08<00:06, 1.84steps/s, loss=(array(0.17124742, dtype=float32),)]
60%|██████ | 18/30 [00:09<00:06, 1.84steps/s, loss=(array(0.17101422, dtype=float32),)]
63%|██████▎ | 19/30 [00:09<00:05, 1.93steps/s, loss=(array(0.17101422, dtype=float32),)]
63%|██████▎ | 19/30 [00:09<00:05, 1.93steps/s, loss=(array(0.1707828, dtype=float32),)]
67%|██████▋ | 20/30 [00:09<00:05, 1.99steps/s, loss=(array(0.1707828, dtype=float32),)]
67%|██████▋ | 20/30 [00:09<00:05, 1.99steps/s, loss=(array(0.17055179, dtype=float32),)]
70%|███████ | 21/30 [00:10<00:04, 2.04steps/s, loss=(array(0.17055179, dtype=float32),)]
70%|███████ | 21/30 [00:10<00:04, 2.04steps/s, loss=(array(0.17032059, dtype=float32),)]
73%|███████▎ | 22/30 [00:10<00:03, 2.07steps/s, loss=(array(0.17032059, dtype=float32),)]
73%|███████▎ | 22/30 [00:10<00:03, 2.07steps/s, loss=(array(0.17008899, dtype=float32),)]
77%|███████▋ | 23/30 [00:11<00:03, 2.11steps/s, loss=(array(0.17008899, dtype=float32),)]
77%|███████▋ | 23/30 [00:11<00:03, 2.11steps/s, loss=(array(0.16985726, dtype=float32),)]
80%|████████ | 24/30 [00:11<00:02, 2.12steps/s, loss=(array(0.16985726, dtype=float32),)]
80%|████████ | 24/30 [00:11<00:02, 2.12steps/s, loss=(array(0.16962597, dtype=float32),)]
83%|████████▎ | 25/30 [00:12<00:02, 2.15steps/s, loss=(array(0.16962597, dtype=float32),)]
83%|████████▎ | 25/30 [00:12<00:02, 2.15steps/s, loss=(array(0.16939586, dtype=float32),)]
87%|████████▋ | 26/30 [00:12<00:01, 2.17steps/s, loss=(array(0.16939586, dtype=float32),)]
87%|████████▋ | 26/30 [00:12<00:01, 2.17steps/s, loss=(array(0.16916767, dtype=float32),)]
90%|█████████ | 27/30 [00:13<00:01, 2.18steps/s, loss=(array(0.16916767, dtype=float32),)]
90%|█████████ | 27/30 [00:13<00:01, 2.18steps/s, loss=(array(0.16894199, dtype=float32),)]
93%|█████████▎| 28/30 [00:13<00:01, 1.90steps/s, loss=(array(0.16894199, dtype=float32),)]
93%|█████████▎| 28/30 [00:13<00:01, 1.90steps/s, loss=(array(0.16871896, dtype=float32),)]
97%|█████████▋| 29/30 [00:14<00:00, 1.98steps/s, loss=(array(0.16871896, dtype=float32),)]
97%|█████████▋| 29/30 [00:14<00:00, 1.98steps/s, loss=(array(0.16849855, dtype=float32),)]
100%|██████████| 30/30 [00:14<00:00, 2.05steps/s, loss=(array(0.16849855, dtype=float32),)]
100%|██████████| 30/30 [00:14<00:00, 2.04steps/s, loss=(array(0.16849855, dtype=float32),)]
# Make a GIF of all images.
imgs = [Image.open(img) for img in image_files]
imgs[0].save(
"mrinufft_learn_traj_multires.gif",
save_all=True,
append_images=imgs[1:],
optimize=False,
duration=2,
loop=0,
)
Results#
model.eval()
final_recons = model(image)
final_traj = model.get_trajectory().detach().cpu().numpy()
fig, axs = plt.subplots(1, 3, figsize=(9, 3))
plot_state(axs, image, final_traj, final_recons)
plt.show()
The learned trajectory above improves the reconstruction quality as compared to the initial trajectory shown above. Note of course that the reconstructed image is far from perfect because of the documentation rendering constraints. In order to improve the results one can start by training it for more than just 5 iterations per decimation level. Also density compensation should be used, even though it was avoided here for CPU compliance. Check out Learn Sampling pattern to know more.
References#
[Cha+16]
N. Chauffert, P. Weiss, J. Kahn and P. Ciuciu, “A Projection Algorithm for Gradient Waveforms Design in Magnetic Resonance Imaging,” in IEEE Transactions on Medical Imaging, vol. 35, no. 9, pp. 2026-2039, Sept. 2016, doi: 10.1109/TMI.2016.2544251.
[Cha+22]
G. R. Chaithya, P. Weiss, G. Daval-Frérot, A. Massire, A. Vignaud and P. Ciuciu, “Optimizing Full 3D SPARKLING Trajectories for High-Resolution Magnetic Resonance Imaging,” in IEEE Transactions on Medical Imaging, vol. 41, no. 8, pp. 2105-2117, Aug. 2022, doi: 10.1109/TMI.2022.3157269.
[GRC23]
Chaithya GR, and Philippe Ciuciu. 2023. “Jointly Learning Non-Cartesian k-Space Trajectories and Reconstruction Networks for 2D and 3D MR Imaging through Projection” Bioengineering 10, no. 2: 158. https://doi.org/10.3390/bioengineering10020158
Total running time of the script: (1 minutes 5.071 seconds)
