BlochSimulators

BlochSimulators is a Julia package for performing Bloch simulations within the context of Magnetic Resonance Imaging. It provides efficient implementations of both the Isochromat Summation model and the Extended Phase Graph (EPG) model to simulate MR signals resulting from custom pulse sequences and k-space trajectories. Simulations can be deployed on various computational resources, with strong support for CUDA-compatible GPU acceleration to achieve high runtime performance. The package is well-suited for simulating dictionaries for MR Fingerprinting or performing forward model evaluations for MR-STAT. For detailed information, please refer to the Documentation.

Installation

BlochSimulators.jl is registered in the General Julia registry. To install the package, press ] in the Julia REPL to enter package mode, followed by either

pkg> add BlochSimulators (if you want to use the package as-is)

or

pkg> dev BlochSimulators (if you want to make modificatios to the source code).

Basic Usage

Here's a minimal example simulating the magnetization for a single voxel using the EPG model with a FISP sequence:

using BlochSimulators

# Define tissue properties (T1=1.0s, T2=0.1s)
parameters = @parameters T1=1.0 T2=0.1

# Use a predefined FISP sequence (replace with your actual sequence definition)
# Note: You might need to define or load a sequence struct first.
# This is a placeholder example.
struct MyFISPSequence <: EPGSimulator end # Placeholder
sequence = MyFISPSequence() # Placeholder

# Simulate magnetization
magnetization = simulate_magnetization(CPU1(), sequence, parameters)

println("Final magnetization (Mx, My, Mz): ", magnetization)

Please see the examples directory and the documentation for more complete examples, including sequence definitions and signal simulation.

Examples

The examples folder contains several example sequence structs as well as k-space trajectory structs. Users of BlochSimulators.jl are encouraged to modify these examples or assemble their own structs with custom sequence/k-space trajectory implementations. Example Julia code on how to use a custom sequence to simulate an MR Fingerprinting dictionary is available here. Another example in which both a sequence and k-space trajectory struct are used to simulate MR signal is available here.

Citation

See CITATION.bib for details on how to cite this work.

---

Developer Guide

This section provides an overview of the package architecture and design decisions to help contributors understand how the code is organized and how to extend it.

Package Architecture Overview

BlochSimulators is built around a separation of concerns design philosophy that cleanly divides the simulation workflow into distinct, composable components:

1. Tissue Properties → What are we simulating? (T₁, T₂, B₁, B₀, etc.)
2. Sequences → How do we manipulate magnetization? (RF pulses, gradients, timing)
3. Operators → Low-level physics engines (Isochromat vs EPG models)
4. Trajectories → How do we encode spatial information? (k-space sampling)
5. Computational Resources → Where do we run? (CPU single/multi-threaded, GPU, distributed)

This modular design allows researchers to mix and match components (e.g., different sequences with different trajectories) without reimplementing core physics.

Core Design Patterns

1. Type-Driven Dispatch and Compile-Time Optimization

The package heavily uses Julia's type system to achieve high performance:

• Parametric types encode critical information at compile time. For example, EPGSimulator{T,Ns} includes both the precision T (Float32/Float64) and the number of configuration states Ns as type parameters.
• Ns as a type parameter is essential: Ns indicates the size of the EPG configuration state matrix Ω (3×Ns). This state matrix is accessed frequently and is the core of EPG simulations, so we cannot have it heap-allocated. Having its size known at compile time enables stack allocation (via StaticArrays) and is critical for GPU kernel compilation—the GPU compiler needs to know array sizes to generate efficient code. This is why we dispatch on Ns as a type parameter rather than a field value.
• Val{Ns} is used in struct fields when Ns needs to be stored but also accessible at compile time for dispatch.

2. Informal Interfaces via Abstract Types

Rather than formal traits (like in Rust), the package uses abstract types to define informal interfaces:

• BlochSimulator{T} → Subtypes must implement:

  • output_eltype(sequence) → What type of output (Complex, Real, Isochromat)? If you just want transverse magnetization, combine the x and y components as a complex number (x + iy). For some sequences where transverse magnetization always aligns with one axis, you can store as reals (reducing memory and computation). If you need full x,y,z components at each timepoint, use Isochromat as the output type.
  • output_size(sequence) → How many time points?
  • simulate_magnetization!(magnetization, sequence, state, parameters) → The actual physics implementation

• AbstractTissueProperties{N,T} → Encapsulates voxel parameters (T₁, T₂, B₁, etc.)

  • Subtypes of FieldVector from StaticArrays for performance
  • Named fields improve code readability while maintaining vector operations
  • The @parameters macro provides a convenient constructor

• AbstractTrajectory{T} → Subtypes must implement:

  • nreadouts(trajectory), nsamplesperreadout(trajectory, idx)
  • to_sample_point(m, trajectory, readout_idx, sample_idx, parameters) → Propagate magnetization with spatial encoding

Why informal interfaces? Julia does not have formal interfaces (yet), but interfaces provide significant benefits for code organization and extensibility. By defining informal interfaces through abstract types and their required methods, developers implementing new sequences or trajectories know exactly what methods to define by looking at the abstract type's documentation.

3. Resource-Based Dispatch for Heterogeneous Computing

The package uses ComputationalResources.jl to abstract computational backends:

• CPU1() → Single-threaded execution
• CPUThreads() → Multi-threaded (respects JULIA_NUM_THREADS)
• CPUProcesses() → Distributed computing (works with DArrays)
• CUDALibs() → GPU execution

Key design decisions:

• The gpu(x) function (using Functors.jl) recursively moves structs to GPU, converting Arrays to CuArrays
• Sequences/trajectories must be marked with @functor to enable GPU transfer
• On GPU, simulations use a fixed THREADS_PER_BLOCK with WARPSIZE (hardcoded to 32 for CUDA devices) threads per voxel
• Type stability and zero-allocation in kernel code is critical for GPU performance

4. The Two-Level Simulation API

The package exposes two main entry points that differ in their level of abstraction:

Level 1: simulate_magnetization(resource, sequence, parameters)

• Simulates magnetization at echo times only, without spatial encoding
• Returns array of size (output_size(sequence), length(parameters))
• Used for MR Fingerprinting dictionaries where k-space trajectories aren't needed
• Each voxel's simulation is independent → embarrassingly parallel

Level 2: simulate_signal(resource, sequence, parameters, trajectory, coordinates, coil_sensitivities)

• Full signal simulation including spatial encoding and multi-coil reception
• First calls simulate_magnetization for all voxels
• Then propagates magnetization to all k-space sample points via trajectory
• Returns array of size (# samples per readout, # readouts, # coils)
• Used for realistic MRI signal simulations

This two-level design avoids code duplication and allows users to choose the appropriate abstraction level.

Code Organization

src/
├── BlochSimulators.jl           # Main module file, includes all components
├── interfaces/                   # Abstract types and informal interfaces
│   ├── tissueproperties.jl      # T₁T₂, T₁T₂B₁, etc. + @parameters macro
│   ├── sequences.jl             # BlochSimulator, IsochromatSimulator, EPGSimulator
│   ├── trajectories.jl          # AbstractTrajectory interface
│   └── coordinates.jl           # Spatial coordinate handling
├── operators/                    # Low-level physics implementations
│   ├── isochromat.jl            # Rodrigues rotation, decay/regrowth for isochromats
│   ├── epg.jl                   # EPG state manipulation (shift, RF, decay, etc.)
│   └── utils.jl                 # Shared utilities
├── simulate/                     # High-level simulation orchestration
│   ├── magnetization.jl         # simulate_magnetization + resource dispatch
│   └── signal.jl                # simulate_signal + trajectory integration
└── utils/                        # Cross-cutting concerns
    ├── gpu.jl                   # gpu() function, Functors integration
    └── precision.jl             # f32/f64 precision conversion

sequences/                        # Example sequence implementations
├── fisp2d.jl, fisp3d.jl         # Gradient-spoiled FISP (EPG-based)
├── pssfp2d.jl, pssfp3d.jl       # pSSFP (EPG-based)
├── generic2d.jl, generic3d.jl   # Generic isochromat simulators
└── adiabatic.jl                 # Adiabatic inversion sequences

trajectories/                     # Example trajectory implementations
├── cartesian.jl, radial.jl, spiral.jl
└── _abstract.jl                 # Shared trajectory utilities

Key principle: Examp