Basics of performing Analytic Simulations#

This notebooks aims to show only the basic concepts of analytic simulations inside the cmrsim frame work. The incorporation of motion is shown in the next notebook.

It is subdivided into the following subsections:

- Define the simulation - Define a digital phantom - Run the simulation - Auxiliary functions to save the simulation config

Imports#

[3]:

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)

[8]:

# 3rd Party dependencies
from pint import Quantity
import pyvista
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline

# Project library cmrsim
import cmrsim
import local_functions

Define the simulation#

As described in the readme and the documentation of the SimulationBase-class, there is the convenient way of subclassing the abstract Base-class. All subclasses need to implement the the _build method, where the composition of building blocks is performed. As the simulation parameters are defined as tf.Variables their values are easily modifyable at any time, even with graph execution.

The _build member-function must return a triple containing an instance of a subclassed EncodingBase module, a CompositeSignalModel (composing all signal operator building blocks) and an instance of a ReconBase module. If the simulation is meant to return the k-space signal, the third entry of the returned triple can be set to None.

The following ExampleSimulation is setup as a single shot (EPI) spin-echo, considering the T2* effects (TE-t, TE+t) with a simple inverse FFT as reconstruction.

[5]:

from typing import Tuple


class ExampleSimulation(cmrsim.analytic.AnalyticSimulation):
    def build(self, matrix_size: Tuple[int, int] = (80, 100),
              field_of_view: Tuple[float, float] = (0.25, 0.3)):
        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=1,
                                                       absolute_noise_std=0.)
        # Signal module construction
        spinecho_module = cmrsim.analytic.contrast.SpinEcho(echo_time=50, repetition_time=10000, expand_repetitions=True)

        centered_times = tf.abs(encoding_module.get_sampling_times() - tf.reduce_max(encoding_module.get_sampling_times()) / 2)
        t2_star_module = cmrsim.analytic.contrast.t2_star.StaticT2starDecay(centered_times, expand_repetitions=False)

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

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

To run the simulation, a Simulator-object needs to be instantiated.

[6]:

field_of_view = Quantity((0.25, 0.3), "m")
matrix_size = (80, 100)
simulator = ExampleSimulation(name='simulation_example', build_kwargs=dict(matrix_size=matrix_size, field_of_view=field_of_view.m_as("m").tolist()))

CMRSeq DeprecationWarning: /tmp/ipykernel_3220715/888904239.py:10 : Consider using the cmrsim.analytic.encoding.GenericEncoding class incombination with a cmrseq.sequence definition

Now display what the requirements for the digital phantoms are, i.e. which physical properties must each material point have assigned:

[7]:

simulator.forward_model.required_quantities, simulator.get_k_space_shape()

[7]:

(('M0', 'T1', 'T2', 'T2star'),
 <tf.Tensor: shape=(2,), dtype=int32, numpy=array([   1, 8000], dtype=int32)>)

Define a digital phantom#

In the following cell, a 2D object is constructed by loading a semantic mask (labels are integers) from the XCAT phantom.

The coordinates are assumed to lie on a regular-grid, therefore the utility of cmrsim.utils.coordinates is an easy way to calculate them.

As described in the README, the digital phantom is defined as OrderedDictionary of numpy arrays, where the keys are named after the entries of the required quantities. Furthermore, the assign_to_numerical_label_map is used to assign values according to the semantic label of the input map.

[9]:

phantom_mesh = local_functions.create_cylinder_phantom()
clear_output()
phantom_mesh

[9]:

Header

Data Arrays

UnstructuredGrid	Information
N Cells	96456
N Points	130204
X Bounds	-1.000e-01, 9.900e-02
Y Bounds	-1.000e-01, 9.900e-02
Z Bounds	-5.000e-03, 2.500e-03
N Arrays	5

Name	Field	Type	N Comp	Min	Max
class	Points	int32	1	-1.000e+00	3.000e+00
M0	Points	float64	1	0.000e+00	9.000e-01
T1	Points	float64	1	0.000e+00	1.000e+03
T2	Points	float64	1	0.000e+00	2.500e+02
T2star	Points	float64	1	0.000e+00	5.000e+01

[10]:

pyvista.close_all()
pyvista.start_xvfb()
plotter = pyvista.Plotter(window_size=(1600, 300), shape=(1, 4))
for idx, sc in enumerate(["class", "M0", "T1", "T2"]):
    plotter.subplot(0, idx)
    plotter.add_mesh(phantom_mesh.copy(), scalars=sc)
    local_functions.add_custom_axes(plotter, scale=0.2)
    box = local_functions.transformed_box(reference_basis=np.eye(3, 3), slice_normal=np.array([0., 0., 1.]), readout_direction=np.array([1., 0., 0.]),
                                          slice_position=Quantity([0., 0., 0.], "m"), extend=Quantity([*field_of_view.m_as("m"), 0.01], "m"))
    plotter.add_mesh(box, color="white", opacity=0.5)
plotter.show(jupyter_backend="static")

CMRSeq Warning : Optimal rotation is not uniquely or poorly defined for the given sets of vectors.

../../_images/example_gallery_analytic_simulation_1_basic_functionality_12_1.png

To pass the phantom into the simulator, the a dictionary is used to instantiate a costum tf.data.Dataset implemented in the BaseDataset class. The dictionary contains numpy arrays corresponding to the required quantities listed above.

[11]:

phantom_dict = dict(
    M0=phantom_mesh["M0"].astype(np.complex64),
    T1=phantom_mesh["T1"].astype(np.float32),
    T2=phantom_mesh["T2"].astype(np.float32),
    T2star=phantom_mesh["T2star"].astype(np.float32),
    r_vectors=phantom_mesh.points.astype(np.float32))

print("Phantom Definition: \t", [f"{k}: {v.shape}" for k, v in phantom_dict.items()])
dataset = cmrsim.datasets.AnalyticDataset(phantom_dict, filter_inputs=True, expand_dimension=True)

Phantom Definition:      ['M0: (130204,)', 'T1: (130204,)', 'T2: (130204,)', 'T2star: (130204,)', 'r_vectors: (130204, 3)']

Run the simulation#

To start the simulation, use the __call__ function of the simulator instance {() - operator}.

[12]:

%time simulation_result = simulator(dataset, voxel_batch_size=10000)
simulation_result.shape

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

Run Simulation:         |XXXXXXXXXXXXXXX|127548/127548  |XXXXXXXXXXXXXXX|1/1 CPU times: user 3.46 s, sys: 699 ms, total: 4.16 s
Wall time: 4.31 s

[12]:

TensorShape([1, 100, 80])

[13]:

f, a = plt.subplots(1, 2, figsize=(8, 6))
magnitude_image = np.squeeze(np.abs(simulation_result))
phase_image = np.squeeze(np.angle(simulation_result))
img1 = a[0].imshow(magnitude_image, cmap="gray")
img2 = a[1].imshow(phase_image, cmap="twilight")
[plt.colorbar(im, ax=ax, location="bottom", pad=0.01, shrink=0.9) for im, ax in zip([img1, img2], a)]
[_.axis("off") for _ in a.flatten()]
[_.set_title(t, fontsize=20) for _, t in zip(a, ("Magnitude", "Phase"))]
f.tight_layout()

../../_images/example_gallery_analytic_simulation_1_basic_functionality_17_0.png

Auxilliary functions#

1. Saving and Loading#

The next cell shows how the simulation with all its configurations can be saved as tf.train.Checkpoint, and afterwards be restored from the checkpoint by calling the classmethod SimulationBase.from_checkpoint(...). To ensure that the parameters are equal, the simulation is redone and the results are compared.

[14]:

import os

# Specify saving path
chkpt = './model_checkpoints/example_checkpoint'
os.makedirs(os.path.dirname(chkpt), exist_ok=True)

# Save simulation to checkpoint and delete old instance
simulator.save(checkpoint_path=chkpt)
del simulator

# construct new intance from checkpoint
new_simulator = ExampleSimulation.from_checkpoint(chkpt)

CMRSeq DeprecationWarning: /tmp/ipykernel_3220715/888904239.py:10 : Consider using the cmrsim.analytic.encoding.GenericEncoding class incombination with a cmrseq.sequence definition

[15]:

%time simulation_load_result = new_simulator(dataset, voxel_batch_size=5000)

print("If this is equal to zero, the simulated results are identical:\n\t Difference of simulated k-space ->",
      tf.reduce_sum(simulation_result - simulation_load_result))

Run Simulation:         |XXXXXXXXXXXXXXX|127548/127548  |XXXXXXXXXXXXXXX|1/1 CPU times: user 3.05 s, sys: 303 ms, total: 3.35 s
Wall time: 2.99 s
If this is equal to zero, the simulated results are identical:
         Difference of simulated k-space -> tf.Tensor((0.073531434-3.5832636e-07j), shape=(), dtype=complex64)

2. Manipulate Parameters#

After instantitation of simulator instances, the variables contained in the submodules can be changed by directly accesing the tf.Variable.assign() method. To show if the results differ after setting different Echo-Times the simulated images are compared again:

[16]:

new_simulator.forward_model.spin_echo.TE.assign([10,])  # in milliseconds
%time n = new_simulator(dataset, voxel_batch_size=5000)
print('Mean absolute difference for different TE:', tf.reduce_mean(tf.abs(simulation_load_result - n)))

new_simulator.forward_model.spin_echo.TE.assign([50,])  # in milliseconds
%time n = new_simulator(dataset, voxel_batch_size=5000)
print('Mean absolute difference for same TE:', tf.reduce_mean(tf.abs(simulation_load_result - n)))

Run Simulation:         |XXXXXXXXXXXXXXX|127548/127548  |XXXXXXXXXXXXXXX|1/1 CPU times: user 1.42 s, sys: 139 ms, total: 1.56 s
Wall time: 1.16 s
Mean absolute difference for different TE: tf.Tensor(1.522993, shape=(), dtype=float32)
Run Simulation:         |XXXXXXXXXXXXXXX|127548/127548  |XXXXXXXXXXXXXXX|1/1 CPU times: user 1.4 s, sys: 152 ms, total: 1.55 s
Wall time: 1.14 s
Mean absolute difference for same TE: tf.Tensor(0.0, shape=(), dtype=float32)

3. Display summary as JSON#

The SimulationBase class includes a functionality to create a report of the simulation configuration (aka submodule Variables) as dictionary. This dictionary can be used to generate a JSON string

[17]:

import json
from IPython.display import JSON

summary_dict = new_simulator.configuration_summary
JSON(json.dumps(summary_dict))

CMRSeq Warning : JSON expects JSONable dict or list, not JSON strings

[17]:

<IPython.core.display.JSON object>