Intra-scan motion

import torch
from cornucopia.utils.py import meshgrid_ij
from cornucopia import (
    IntraScanMotionTransform,
    ArrayCoilTransform,
    SmallIntraScanMotionTransform,
    RandomSlicewiseAffineTransform,
)
import  matplotlib.pyplot as plt

Generate a digital phantom

shape = [128, 128]
radius = torch.stack(meshgrid_ij(*[torch.arange(s).float() for s in shape]), -1)
radius -= (torch.as_tensor(shape).float() - 1) / 2
radius = radius.square().sum(-1).sqrt()

mag = torch.zeros_like(radius, dtype=torch.float32)
mag[radius < 48] = 1
mag[radius < 44] = 2
mag[radius < 24] = 3
mag = mag[None]  # channels dimension

plt.figure(figsize=(10, 10))
plt.imshow(mag.squeeze(), cmap='gray', interpolation='nearest')
plt.axis('off')
plt.title('Signal magnitude')
plt.colorbar()
plt.show()

Output figure

Build k-space from 4 shots acquired with different object position

trf = IntraScanMotionTransform(coils=ArrayCoilTransform())
sos = trf(mag)

plt.figure(figsize=(10, 10))
plt.imshow(sos.squeeze(), cmap='gray', interpolation='nearest')
plt.axis('off')
plt.title('Signal sum-of-squares')
plt.colorbar()
plt.show()

Output figure

We can switch to a random sampling pattern

trf = IntraScanMotionTransform(coils=ArrayCoilTransform(), pattern='random')
sos = trf(mag)

plt.figure(figsize=(10, 10))
plt.imshow(sos.squeeze(), cmap='gray', interpolation='nearest')
plt.axis('off')
plt.title('Signal sum-of-squares')
plt.colorbar()
plt.show()

Output figure

Same experiment, but this time motion happens across a "slice" axis (no FFT involved)

plt.show()
trf = IntraScanMotionTransform(coils=ArrayCoilTransform(), freq=False)
sos = trf(mag)

plt.figure(figsize=(10, 10))
plt.imshow(sos.squeeze(), cmap='gray', interpolation='nearest')
plt.axis('off')
plt.title('Signal sum-of-squares')
plt.colorbar()
plt.show()

Output figure

Finally, sample some small motion

trf = SmallIntraScanMotionTransform()

shape = [4, 4]
plt.rcParams["figure.figsize"] = (15, 15)

plt.figure(figsize=(10, 10))
for i in range(shape[0] * shape[1]):
    plt.subplot(*shape, i+1)
    plt.imshow(trf(mag).squeeze(), cmap='gray', interpolation='nearest')
    plt.axis('off')
plt.suptitle('Small motion')
plt.show()

Output figure

We also have a transform that applies different amounts of (through plane) motion to different slices. It can be useful to model motion occuring during the acquisition of a stack of slices.

It is quite similar to IntraScanMotionTransform(freq=False), but implemented in a much more efficient way.

trf = RandomSlicewiseAffineTransform()

shape = [4, 4]
plt.rcParams["figure.figsize"] = (15, 15)

plt.figure(figsize=(10, 10))
for i in range(shape[0] * shape[1]):
    plt.subplot(*shape, i+1)
    plt.imshow(trf(mag).squeeze(), cmap='gray', interpolation='nearest')
    plt.axis('off')
plt.suptitle('Slicewise motion')
plt.show()

Output figure

trf = RandomSlicewiseAffineTransform(spacing=4, subsample=2)

shape = [4, 4]
plt.rcParams["figure.figsize"] = (15, 15)

plt.figure(figsize=(10, 10))
for i in range(shape[0] * shape[1]):
    plt.subplot(*shape, i+1)
    plt.imshow(trf(mag).squeeze(), cmap='gray', interpolation='nearest')
    plt.axis('off')
plt.suptitle('Slicewise motion')
plt.show()

Output figure

# It is also possible to provide pre-computed motion trajectories 
# that get interpolated with cubic splines

from cornucopia.random import Fixed
rotations = Fixed(torch.rand([5, 1])*15)
translations = Fixed(torch.rand([5, 2])*0.1)
trf = RandomSlicewiseAffineTransform(
    rotations=rotations, translations=translations,
    bulk_translations=0, bulk_rotations=0,
)

shape = [4, 4]
plt.rcParams["figure.figsize"] = (15, 15)

plt.figure(figsize=(10, 10))
for i in range(shape[0] * shape[1]):
    plt.subplot(*shape, i+1)
    plt.imshow(trf(mag).squeeze(), cmap='gray', interpolation='nearest')
    plt.axis('off')
plt.suptitle('Slicewise motion')
plt.show()

Output figure