Bloch Simulation #

Note

Notebooks demonstrating the content described below can be found by navigating to Examples in the navigation bar at the top of the page

1. Magnetization Rotation #

As described in the input defintion section, cmrsim assumes every particle/iso-chromate to bear the property magnetization with the entries \((M_{xy}^{+}, M_{xy}^{-}, M_{z})\), representing the magnetization in the rotating frame. The actual simulation consists of the three common operations performed on these magnetization vectors, namely:

Rotation around the x-axis due to RF-application
Rotation around the z-axis due to Gradient-application
Relaxation according to the particle properties (\(T1, T2\) )

The application of these three operations are implemented in the cmrsim.bloch.BaseBlochOperator. Therefore, subclassing the BaseBlochOperator allows to calculate the values for rotation and relaxation according to your needs, while reusing the base-implementation. For the most general case of a Bloch-simulation, the GeneralBlochOperator performs a trapezoidal numerical integration of the phase

\[\phi(t) = \gamma \int_{0}^{t} \vec{G}(t') \vec{r}(t') dt',\]

and the flip angle per time-step:

\[\alpha(t) = \gamma \int_{0}^{t} B_{1}(t') dt'.\]

Effects, that alter the phase per time-increment can be specified as submodule (see section 3.4). Motion of particles is incorporated by providing an instance of Trajectory module (see section 2), to the simulation, that increments the particle position according to its update rule. See section 3 to learn how to instantiate and call the GeneralBlochOperator. This section also explains how the acquired signal is saved and accessed after simulation.

2. The Trajectory Module #

The trajectory module contains the logic determining the next particle position given the current one. CMRsim assumes these modules to be subclasses of the tf.Module class. The only requirements that a trajectory module needs to meet for being compatible with the general bloch simulation module are:

Implement the increment_particles method with the signature
(self, r_old: tf.Tensor, dt: tf.Tensor, **kwargs)-> (r_new: tf.Tensor, additional_lookups: dict)
The increment_particles method must be compatible with the tf.function() decorator, even better if the jit_compile=True argument is feasible.

In short, the tf.function() decorator creates a computational graph allowing faster GPU execution within Tensorflow2. In many cases, just using the tensorflow functionality similarly to how one would implement it with NumPy is sufficient to meet criterion (ii). For more information about the tf.function() decorator consult the tensorflow documentation.

The simplest trajectory module would just return the identy operation for the old position. This is how Bloch simulations with static isochromats are defined in CMRsim. Particles moving in a field of spatially varying velocity can be represented by a module that looks up the velocity at the current position and adds an positional increment of \(\delta \vec{r} = \vec{v}\delta t\). This look up (with nearest neighbour as well as trilinear interpolation) is already implemented in tensorflow meeting the criterion (ii) and can be used out of the box. To see how this works checkout the implementation of e.g. :code:’cmrsim.trajectory.FlowTrajectory’. To see what Trajectory Modules are already shipped with cmrsim checkout the API documentation and the examples section.

3. The GeneralBlochOperator #

3.1 Instantiation #

The fundamental idea of the GeneralBlochOperator is to provide a convenient way of running a bloch-simulation for a sequence that is defined on a (not necessarily uniform) temporal raster containing at least one of the channels (RF, Gradients, ADC). For detailed information checkout the API-documentation. In short, all channels are specified as NumPy arrays and can be handed to the Module constructor as keyword argument. If not provided the simulation will just not perform rotations/signal accumulation. The recommended way of defining the sequences is using the CMRseq package, which is also featured in the example notebooks. In principle every Numpy array will be usable though.

Note

As the computational graph is constructed on calling the module, not specifying the channels does not add the overhead of performing conditional checks on each time-step.

Instantiation example

Assume the gridded waveforms to be defined previously. The gradient definition contains the channels (x, y, z), the adc is defined with channels (is_on, phase) and the RF waveform is of dtype complex. The meaning of n_reps is expained in 3.2

time_raster.shape, grad_grid.shape, rf_grid.shape, adc_grid.grid
>>> (n_reps, t), (n_reps, t, 2), (n_reps, t), (n_reps, t, 2)

bloch_mod = cmrsim.bloch.GeneralBlochOperator(name="example",
                                              device="GPU:0,
                                              gamma=system_specs.gamma_rad.m_as("rad/mT/ms"),
                                              time_grid=time_raster,
                                              gradient_waveforms=grad_grid,
                                              rf_waveforms=rf_grid,
                                              adc_events=adc_grid)

3.2 Parallel vs Sequential execution #

The arguments containing the waveforms have a leading axis n_reps, which is used to efficiently execute multiple repetitions. This means the values in the arrays can change but the temporal raster hast to be identical. Examples for this would be varying phase-encoding gradients over consecutive TRs or Diffusion/Velocity encoding waveforms with varying scaling. The reason for this being efficient is, the internal working of the tf function and its memory layout on the GPU (which should not be a concern for most people). If no repeating elements are present, this axis has the axis 1 and the following explanations are irrelevant.

Executing the simulation is done by calling the instantiated module. For a detailed description of the calling signature look into the API documentation. Code examples for calling the module are given below as well as in the example notebooks. In short, the call function takes the inital magnetizations, inital positions and all required particle properties alongside the trajectory module instance and returns the resulting magnetization vectors as well as the new positions.

Parallel: In some simulation experiments it might be desirable to not compute particle updates over and over again for varying sequence parameters (e.g. VENC) so the GeneralBlochOperator allows to specify the boolean run_parallel argument. If true the integration is evaluated for a single set of trajectories but for all waveforms. The corresponding signal is then accumulated as described in 3.4.

Call with parallel execution

For memory reasons batch the digital phantom and iterate over it

properties = ...  # From input definition
n_samples = tf.shape(module.gradient_waveforms)[1]
n_repetitions = tf.shape(module.gradient_waveforms)[0]
for phantom_batch in tqdm(tf.data.Dataset.from_tensor_slices(properties).batch(int(3e5))):
    m, r = bloch_mod(trajectory_module=trajectory_module, run_parallel=True, **phantom_batch)

Sequential: If the simulation aims at investigating the magnetization history or if re-seeding / dropping of particles before every rf-pulse pulse is required, the repetitions should be run sequentially. To accumulate the signal according to the correct repetition, the calling signature includes the argument repetition_index.

Call with sequential execution

After each repetition the new magnetization is passed to the module again

properties = ...  # From input definition
n_samples = tf.shape(module.gradient_waveforms)[1]
n_repetitions = tf.shape(module.gradient_waveforms)[0]

for phantom_batch in tqdm(tf.data.Dataset.from_tensor_slices(properties).batch(int(3e5))):
    r = phantom_batch.pop("initial_position")
    m = phantom_batch.pop("magnetization")
    for rep_idx in tqdm(range(n_repetitions)):
        m, r = module(trajectory_module=trajectory_module, initial_position=r, magnetization=m,
                      repetition_index=rep_idx, run_parallel=False, **phantom_batch)

3.3 Signal accumulation #

For each step of the time raster, the gridded adc input must specify if the step should be recorded or not. If the the value of the (is_on) channel is greater than zero, the sum over all isochromates transverse magnetization weighted by its effective magnetization is evaluated and added to the accumulator bloch_module.time_signal_acc[rep_idx][step_idx]. As the fourier integral can be decomposed into sums of batches, it is trivial to split the phantom and progressively summing up the signal.

After execution the signal can be directy accessed, and its value must be manually reset before re-using the simulation module.

3.4 Submodules #

Currently, submodule adding a phase-contribution to the gradient application are supported. The available submodules are listed in API documentation as cmrsim.bloch.submodules. A list of submodules can be passed as argument on instantiation of the GeneralBlochOperator object. For an example using off-resonance in a BSSFP resulting in banding artifacts checkout the example notebooks.