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

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

import os
os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3"
os.environ['CUDA_VISIBLE_DEVICES'] = "3"

import tensorflow as tf
gpu = tf.config.get_visible_devices("GPU")
if gpu:
    tf.config.experimental.set_memory_growth(gpu[0], True)

from copy import deepcopy
import base64

# 3rd Party dependencies
from pint import Quantity
from tqdm.notebook import tqdm
import pyvista
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display, HTML, clear_output

%matplotlib inline

# Project library cmrsim
import cmrsim
import local_functions

Define Slice Position / Imaging Parameters

# 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.

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

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.

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

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

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="gray")
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="gray")
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

# 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="black")
    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"))

import imageio
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}" />'))

/usr/local/lib/python3.11/dist-packages/pyvista/core/filters/data_set.py:3483: PyVistaDeprecationWarning: probe filter is deprecated and will be removed in a future version.
            Use sample filter instead.
            If using `mesh1.probe(mesh2)`, use `mesh2.sample(mesh1)`.

  warnings.warn(

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.

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))

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

WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1712954396.740806 3231999 device_compiler.h:186] Compiled cluster using XLA!  This line is logged at most once for the lifetime of the process.

Run Simulation:         |XXXXXXXXXXXXXXX|1581304/1581304        |XXXXXXXXXXXXXXX|151/151 CPU times: user 53.9 s, sys: 9.15 s, total: 1min 3s
Wall time: 1min

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)