Multiple Coils and advanced Encoding Features#

In this notebook, an analytic simulation using coil-sensitivity weighting and arbitrary 3D slice positioning / orientation is demonstrated. Additionally, it is shown how to add complex noise-samples to the kspace with varying standard deviations.

The content is structured as follows:

  1. Imports

  2. Define Slice Position / Imaging parameters

  3. Define and Instantiate Simulator

  4. Construct Phantom and Illustrate Setup

  5. Run Simulation

  6. Display Results

Imports#

[1]:

import sys
sys.path.append("../")
sys.path.insert(0, "../../")

from copy import deepcopy
import base64
from IPython.display import display, HTML, clear_output
import imageio

import tensorflow as tf
tf.get_logger().setLevel('ERROR')
physical_devices = tf.config.list_physical_devices('GPU')
if physical_devices:
    print("Available GPUS: \n\t", "\n\t".join([str(_) for _ in physical_devices]))
    tf.config.experimental.set_visible_devices(physical_devices[0], 'GPU')
    tf.config.experimental.set_memory_growth(physical_devices[0], True)

import pyvista
from pint import Quantity
import numpy as np
from tqdm.notebook import tqdm
import matplotlib.pyplot as plt
%matplotlib inline

import cmrsim
import cmrseq
import local_functions
clear_output()

Define Slice Position / Imaging Parameters#

[2]:

# Define Fourier-Encoding parameters
field_of_view = Quantity([25, 25], "cm")
matrix_size = (151, 151)
res = field_of_view / matrix_size

# Define imaging slice parameters
slice_thickness = Quantity(6, "mm")
readout_direction = np.array([1., 0., 0.])
slice_normal = np.array([0., 1., 1.])
slice_normal /= np.linalg.norm(slice_normal)
slice_position_offset = Quantity(6., "cm")

# Define RF-exication parameters (coupled with seeding slice atm)
flip_angle=Quantity(np.pi/4, "rad")
pulse_duration = Quantity(0.6, "ms")

print(f"Resolution: {res :1.3f} \nSlice Position {slice_normal * slice_position_offset :1.3}")

Resolution: [0.166 0.166] centimeter
Slice Position [0.0 4.24 4.24] centimeter

Define and Instantiate Simulator#

To incorporate coil-sensitivities weighting in the simulation the corresponding CoilSensitivity contrast module is included. The coil-maps can be specified as an arbitrary 3D map with complex values, here an ideal bird-cage configuration with 8 receive elements is used which is provided as classmethod of the sensitivity module. On calling the module instance, a trilinear interpolation of the coil-maps is performed for each location specified in the r_vectors argument. Therefore, the repetition-axis of the result is expanded 8-fold (one image per coil-channel) automatically. Additionally, even for moving object the simulator/coil configuration does not have to be recomputed.

Later in this notebook, the 3D rendering of the phantom placed inside the module-sensitivity-maps is shown in combination with the imaging slice position.

[3]:

from typing import Tuple

class ExampleSimulation(cmrsim.analytic.AnalyticSimulation):
    def build(self, matrix_size: Tuple[int, int] = (150, 150), field_of_view: Tuple[float, float] = (0.25, 0.3),
              slice_normal: np.ndarray = np.array([0., 0., 1.]), slice_position: Quantity = Quantity([0., 0., 0.], "m"),
              readout_direction:  np.ndarray = np.array([1., 0., 0.])):
        n_kx, n_ky = matrix_size

        # Encoding definition / To enable result comparisons, the noise is set to 0 here
        encoding_module = cmrsim.analytic.encoding.EPI(field_of_view=field_of_view, sampling_matrix_size=matrix_size,
                                                       read_out_duration=0.4, blip_duration=0.01, k_space_segments=n_ky,
                                                       absolute_noise_std=0.)

        # For further explanation of the orientation matrix see below:
        encoding_module.set_orientation_matrix(slice_position=slice_position, slice_normal=slice_normal,
                                               readout_direction=readout_direction)

        # Signal module construction  (1 repetition - Can be set arbitrarily)
        spinecho_module = cmrsim.analytic.contrast.SpinEcho(echo_time=tf.constant([50., ]),
                                                            repetition_time=tf.constant([3000., ]),
                                                            expand_repetitions=True)
        # Coil sensitivity construction
        coil_module = cmrsim.analytic.contrast.CoilSensitivity.from_3d_birdcage(map_shape=(100, 100, 150, 8),
                                                        map_dimensions=Quantity([[-15, 15], [-15, 15], [-20, 20]], "cm").m_as("m"),
                                                        relative_coil_radius=1.,
                                                        device="GPU:0")


        # Forward model composition
        forward_model = cmrsim.analytic.CompositeSignalModel(spinecho_module, coil_module)

        # Reconstruction
        recon_module = cmrsim.reconstruction.Cartesian2D(sample_matrix_size=(n_kx, n_ky), padding=None)
        return forward_model, encoding_module, recon_module

[4]:

build_kwargs=dict(matrix_size=matrix_size, field_of_view=field_of_view.m_as("m").tolist(),
                  slice_normal=slice_normal, slice_position=slice_normal*slice_position_offset,
                  readout_direction=readout_direction)
simulator = ExampleSimulation(name='simulation_example', build_kwargs=build_kwargs)

clear_output()

Construct Phantom and Illustrate Setup#

Set up spherical phantom#

The local function definition is used to create a pyvista mesh which contains a spherical pahntom with air-inclusions. The values for MR-relevant properties are manually set after selecting a slice from the phantom. The phantom slice is re-sampled to a resolution at least 5-fold (per direction) finer than the imaging resolution to avoid discretization artefacts (inverse-crime effects). Furthermore, the regular grid positions are slightly pertubed to further reduce discretization artefacts.

[5]:

phantom_ = local_functions.create_spherical_phantom(dimensions=(61, 61, 61), spacing=(0.004, 0.004, 0.004), big_radius=0.12)
phantom_["in_mesh"] = np.ones(phantom_.points.shape[0])
phantom = cmrsim.datasets.RegularGridDataset.from_unstructured_grid(phantom_, pixel_spacing=Quantity([1, 1, 1], "mm"),
                                                                    padding=Quantity([5, 5, 5], "cm"))
clear_output()

Select Slice according to definition at the top of the notebook

[6]:

voxels = phantom.select_slice(slice_normal, spacing=Quantity([0.25, 0.25, 2], "mm"),
                              slice_position=slice_normal*slice_position_offset,
                              field_of_view=Quantity([*field_of_view.m_as("m"), slice_thickness.m_as("m")], "m"),
                              readout_direction=readout_direction, in_mps=False)
voxels = voxels.threshold(0.1, "in_mesh")

Illustrate Slice-Position#

The particle positions fed into the simulation can be defined in the global coordinate system. Internally, all analytic encoding modules are assuming MPS/Slice-coordinates by default. If all contrast modules assume global coordinates, it is possible to set the slice-position of the encoding module, such that the corresponding transformation is applied before spatially integrating the object.

Below, the global slice-configuration as well as the MPS-transformation is illustrated

[7]:

mat = simulator.encoding_module.ori_matrix.read_value().numpy()
r_back = np.einsum('ij, nj -> ni', mat[:3, :3].T, voxels.points - mat[:, 3])

pyvista.close_all()
pyvista.start_xvfb()
plotter = pyvista.Plotter(off_screen=True, window_size=(1200, 400), shape=(1, 2))
plotter.add_mesh(voxels, color="white")
plotter.add_mesh(phantom_, cmap="plasma", opacity=0.1)
local_functions.add_custom_axes(plotter, scale=0.3)
box = local_functions.transformed_box(np.eye(3, 3), slice_normal, readout_direction, slice_normal * slice_position_offset,
                                      extend=Quantity([*field_of_view.m_as("m"), slice_thickness.m_as("m")], "m"))
plotter.add_mesh(box, color="orange", opacity=0.5)
plotter.add_title("Slice Position / Sliced Phantom \n in global coorinates", font_size=8)
plotter.subplot(0, 1)
plotter.add_mesh(pyvista.PolyData(r_back), color="white")
local_functions.add_custom_axes(plotter, scale=0.3)
box = local_functions.transformed_box(np.eye(3, 3), np.eye(3, 3)[2], np.eye(3, 3)[0], Quantity([0, 0, 0], "m"),
                                      extend=Quantity([*field_of_view.m_as("m"), slice_thickness.m_as("m")], "m"))
plotter.add_mesh(box, color="orange", opacity=0.5)
plotter.add_title("Isochromates in MPS coorinates", font_size=8)
plotter.screenshot("temp.png")

b64 = base64.b64encode(open("temp.png",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))

/scratch/jweine/cmrsim/notebooks/analytic_simulation/../local_functions.py:43: UserWarning: Optimal rotation is not uniquely or poorly defined for the given sets of vectors.
  rot, _ = Rotation.align_vectors(reference_basis, new_basis)

Illustrate 3D Coil-sensitivity#

[11]:

# Load the coil map from the simulator instance:
coil_maps = simulator.forward_model.coil_sensitivities.coil_maps

# For plotting split into absolute/phase and create pyvista meshes
coilmaps_magn = {f"mag_coil_{i}": np.abs(coil_maps[..., i]) for i in range(coil_maps.shape[-1])}
coilmaps_phase = {f"phi_coil_{i}": np.angle(coil_maps[..., i]) * 0 for i in range(coil_maps.shape[-1])}
coil_maps = dict(**coilmaps_magn, **coilmaps_phase)
coil_dataset = cmrsim.datasets.RegularGridDataset.from_numpy(data=coil_maps, extent=Quantity([60, 60, 90], "cm"))
cylinder = pyvista.Cylinder(center=(0., 0., 0.), direction=(0., 0., 1.), radius=0.15, height=0.4)
cylinder = coil_dataset.mesh.clip_surface(cylinder.extract_surface())
full_sphere = coil_dataset.mesh.probe(phantom_)

# Create the images
images = []
for coil_id in tqdm(range(8), desc="Rendering images:"):
    pyvista.close_all()
    pyvista.start_xvfb()
    plotter = pyvista.Plotter(off_screen=True, window_size=(300, 300))
    plotter.add_mesh(voxels, color="white")
    plotter.add_mesh(full_sphere, scalars=f"mag_coil_{coil_id}", cmap="plasma", opacity=0.9)
    plotter.add_mesh(cylinder, scalars=f"mag_coil_{coil_id}", cmap="plasma", show_scalar_bar=False, opacity=0.8)
    local_functions.add_custom_axes(plotter, scale=0.3)
    plotter.add_title(f"Coil #{coil_id}", font_size=8)
    images.append(plotter.screenshot(f"coil_rendering.png"))

ims = np.concatenate(np.concatenate(np.stack(images).reshape(2, 4, 300, 300, 3), axis=1), axis=1)
imageio.v3.imwrite("coil_rendering.png", ims)
b64 = base64.b64encode(open("coil_rendering.png",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))

Run Simulation#

First assign MR-relevant physical properties to the particles. The spatial pertubation of particles is contained in the RegularGridDataset class, which is the reason to wrap the pyvista mesh with an instance of this class. Finally, we create an instance of an AnalyticDataset to pass the input to the simulator.

[10]:

n_particles = voxels.points.shape[0]
voxels["M0"]=np.ones(n_particles, dtype=np.complex64)
voxels["T1"]=np.ones(n_particles, dtype=np.float32) * 1000
voxels["T2"]=np.ones(n_particles, dtype=np.float32) * 200

slice_dataset = cmrsim.datasets.RegularGridDataset(voxels)
properties = slice_dataset.get_phantom_def(filter_by="in_mesh", keys=["M0", "T1", "T2"],  perturb_positions_std=Quantity(0.005, "mm").m_as("m"))

dataset = cmrsim.datasets.AnalyticDataset(properties, filter_inputs=True, expand_dimension=True)
print([[f"{k}: {v.shape}" for k, v in di.items()] for di in dataset(batchsize=1000).take(1)])

[['r_vectors: (1000, 1, 1, 3)', 'M0: (1000, 1, 1)', 'T1: (1000, 1, 1)', 'T2: (1000, 1, 1)']]

By assigning values to the absolute_noise_std member-variable of the encoding module, the standard deviation of complex valued gaussian noise-samples that is added to the simulated k-space samples is determined. If the value is < 0 no noise is added. If you want to add two indepented noise-samples with the same standard deviation, you can specifiy the same value multiple times. The different noise-instances are accesible by indexing the simulation result on the last non-image axis. (E.g. 1 axis (zero-based) for eight coils and a matrix size of (151, 151) -> (8, #noise, 151, 151))

[12]:

simulator.encoding_module.absolute_noise_std.assign([-1, 5, 10, 20])
%time images_all = simulator(dataset, voxel_batch_size=30000)

2022-12-22 17:25:26.963894: I tensorflow/core/grappler/optimizers/data/replicate_on_split.cc:32] Running replicate on split optimization
2022-12-22 17:25:29.471816: I tensorflow/compiler/xla/service/service.cc:173] XLA service 0x7fc1c800c4a0 initialized for platform CUDA (this does not guarantee that XLA will be used). Devices:
2022-12-22 17:25:29.471899: I tensorflow/compiler/xla/service/service.cc:181]   StreamExecutor device (0): NVIDIA TITAN RTX, Compute Capability 7.5
2022-12-22 17:25:29.521716: W tensorflow/compiler/tf2xla/kernels/assert_op.cc:38] Ignoring Assert operator linear_interpolation_lookup/Assert/Assert
2022-12-22 17:25:29.526699: I tensorflow/compiler/mlir/tensorflow/utils/dump_mlir_util.cc:268] disabling MLIR crash reproducer, set env var `MLIR_CRASH_REPRODUCER_DIRECTORY` to enable.

Run Simulation:         |---------------|0/1650589      |---------------|1/151

2022-12-22 17:25:29.677078: I tensorflow/tsl/platform/default/subprocess.cc:304] Start cannot spawn child process: No such file or directory
2022-12-22 17:25:29.777608: I tensorflow/compiler/jit/xla_compilation_cache.cc:477] Compiled cluster using XLA!  This line is logged at most once for the lifetime of the process.

Run Simulation:         |XXXXXXXXXXXXXX-|1650000/1650589        |XXXXXXXXXXXXXXX|151/151

2022-12-22 17:26:46.799080: W tensorflow/compiler/tf2xla/kernels/assert_op.cc:38] Ignoring Assert operator linear_interpolation_lookup/Assert/Assert

Run Simulation:         |XXXXXXXXXXXXXXX|1650589/1650589        |XXXXXXXXXXXXXXX|151/151 CPU times: user 1min 36s, sys: 16.9 s, total: 1min 53s
Wall time: 1min 20s

[13]:

plt.style.use(['dark_background'])
print(images_all.shape)
for noise_idx in range(images_all.shape[1]):
    images = np.squeeze(images_all[:, noise_idx])
    magnitude, phase = np.abs(images), np.angle(images)

    fig, axes = plt.subplots(2, images.shape[0], figsize=(15, 5))
    axes[0, 0].set_ylabel("Magnitude")
    axes[0, 0].set_ylabel("Phase")

    for ax, im in zip(axes[0], magnitude):
        art = ax.imshow(im, cmap="gray", origin="lower")
        plt.colorbar(art, ax=ax, location="bottom", orientation="horizontal", pad=0.01)
        ax.axis("off")

    for ax, im in zip(axes[1], phase):
        art = ax.imshow(im, cmap="twilight", origin="lower")
        plt.colorbar(art, ax=ax, location="bottom", orientation="horizontal", pad=0.01)
        ax.axis("off")
    fig.suptitle(f"Noise level {simulator.encoding_module.absolute_noise_std.read_value()[noise_idx]}")
    fig.tight_layout()
    fig.savefig(f"temp{noise_idx}.png")
    plt.close(fig)

imageio.v3.imwrite("temp.png", np.concatenate([imageio.v3.imread(f"temp{i}.png")
                                               for i in range(images_all.shape[1])], axis=0))
b64 = base64.b64encode(open("temp.png",'rb').read()).decode('ascii')
display(HTML(f'<img src="data:image/gif;base64,{b64}" />'))

(8, 4, 151, 151)