# Requires Python 3.9 or higher
pip install MRzeroCore

MR-zero is a framework for easy MRI sequence optimization and development of self-learning sequence development strategies. The vision is documented in this paper. These goals are backed by a modern Bloch simulation (see this paper). More material can be found in the literature.

Quick start

For a quick introduction, look at the MR-zero Playground!

An ever growing collection of Jupyter Notebook demonstrate core capabilities of MR-zero. These notebooks run on Google Colab - directly in your browser, no installation needed!

Alternatively, install MR-zero Core locally:

# Requires Python 3.9 or higher
pip install MRzeroCore

Simulate any Pulseq sequence in one line

MR-zero makes MRI sequence simulation easy - just one line for simulating .seq files:

import MRzeroCore as mr0

# That's it - automatic phantom download and simulation!
signal, ktraj_adc = mr0.util.simulate('your_sequence.seq')

Even simpler with PyPulseq - no need to worry about writing .seq files:

import pypulseq as pp
# Create sequence with PyPulseq
seq = build_your_sequence()
# ... build sequence ...
signal, ktraj_adc = mr0.util.simulate(seq)
# ... reconstruct, e.g.: with a NUFFT of signal and ktraj_adc

Further documentation

MR-zero uses PyTorch for computation, enabling GPU compute automatic differentiation with backpropagation. This means you can easily develop your own loss functions, sequence building code and more - and then optimize any input parameters efficiently with gradient descent. In the example above, this could mean to simply extend build_your_sequence() by sequence parameters like flip angles and writing an image based loss function for the reconstruction. Then you can lie back and let the computer find the best pulse train.

How this can be done will be explained in the following pages, where this documentation lists the most important ideas and applications of MR-zero. If you think that something is missing or misleading, pleas open an issue on GitHub or directly via the button on the top right on every documentation page.

Playground

Welcome to Playground MR0, a playground to share, vary and simulate MR sequences. MR sequences are written in the Pulseq standard using the pypulseq library. Pulseq files are simulated with the efficient Phase Distribution Graph Bloch simulation. Here we share links to example colabs that contain various MR sequences or let you upload your own seq file for simulation.

Many of the examples are build using PyPulseq and simulate the resulting .seq files with MR0. These .seq files could also be measured on any MRI scanner using a Pulseq interpreter.

Code and simulate PyPulseq

Notebook
Free Induction Decay	mr0_FID_seq.ipynb
Spin Echo CPMG	mr0_SE_CPMG_seq.ipynb
Stimulated Echo 3 pulses - 5 echoes	mr0_STE_3pulses_5echoes_seq.ipynb
FLASH 2D sequence	mr0_FLASH_2D_seq.ipynb
GRE EPI 2D sequence	mr0_EPI_2D_seq.ipynb
DWI SE EPI 2D sequence	mr0_DWI_SE_EPI.ipynb
Diffusion prepared STEAM	mr0_diffusion_prep_STEAM_2D_seq.ipynb
RARE 2D sequence	mr0_RARE_2D_seq.ipynb
TSE 2D sequence	mr0_TSE_2D_multi_shot_seq.ipynb
Interactive GRE to FLASH	mr0_GRE_to_FLASH.ipynb
balanced SSFP sequence	mr0_bSSFP_2D_seq.ipynb
DREAM STE for B0, B1, TxRx mapping	mr0_DREAM_STE_seq.ipynb
DREAM STID for B0, B1, TxRx mapping	mr0_DREAM_STID_seq.ipynb
Pulseq with RF shimming	pulseq_rf_shim.ipynb

Plot and simulate predifined .seq files

Notebook
Simulate pypulseq example sequences	mr0_pypulseq_exmpls_seq.ipynb
Simulate own uploaded seq files	mr0_upload_seq.ipynb

MR-zero optimization

Gradient descent optimizations using automatic differentiation by backpropagation. Some notebooks use pulseq-zero for optimizable sequence definitions with PyPulseq.

Notebook
IR FLASH 2D sequence for T1 mapping using a fit	mr0_opt_FLASH_2D_IR_Fit_T1.ipynb
IR FLASH 2D sequence for T1 mapping using a NN	mr0_opt_FLASH_2D_IR_voxelNN_T1.ipynb
FLASH flip angle opt. for PSF (with pulseq-zero)	Pulseq_zero_FLASH_FAopt_PSFtask.ipynb
TSE flip angle opt. for SAR (with pulseq-zero)	Pulseq_zero_TSE_FAopt_SARtask.ipynb
DESC with pulseq-zero	pulseq_zero_DESC_demo.ipynb

MR-double-zero optimization

Gradient-free optimization with nevergrad

Notebook
Ernst angle optimization	mr00_FLASH_2D_ernstAngle_opt.ipynb

MR plot wall of fame

famous historic plots recreated

MR0 example notebooks

The following sequences are examples of how to realize various tasks in MR-zero rather than demonstrations of specific MRI sequences.

Notebook
Pure MR0 FLASH	flash.ipynb
pulseq FLASH	pulseq_flash.ipynb
pulseq pTx FLASH	pulseq_sim_pTx.ipynb

Integration

MRI simulations are usually not used in isolation. Instead, simulated sequences are written with various tools, optimized flip angle trains are measured, quantified phantoms are viewed in different programs. The following chapters list currently existing integrations with MR-zero and show how to utilize them.

Pulseq Integration

MRzero Core provides seamless integration with the Pulseq standard, making it easy to simulate any Pulseq sequence file in just one line of code.

Get Started Immediately

Simplest possible - one line simulation:

import MRzeroCore as mr0

# That's it - simulate any .seq file with automatic phantom download!
signal, ktraj_adc = mr0.util.simulate(seq)  # Downloads standard phantom if not found

seq.plot(plot_now=False)
mr0.util.insert_signal_plot(seq=seq, signal =signal.numpy())  # show sequence plots with simulated signals inserted in ADC plot.
plt.show()

Full control with prepass-mainpass:

import MRzeroCore as mr0
obj_p = mr0.VoxelGridPhantom.brainweb('subject05.npz') #load a phantom object from file
obj_p = obj_p.interpolate(64, 64, 32).slices([15])
obj_p = obj_p.build()

seq0 = mr0.Sequence.import_file('your_sequence.se')  # Read in the Pulseq seq file

graph = mr0.compute_graph(seq0, obj_p, 2000, 1e-4)  # pre-pass
signal = mr0.execute_graph(graph, seq0, obj_p, print_progress=False) # main-pass

seq.plot(plot_now=False)
mr0.util.insert_signal_plot(seq=seq, signal =signal.numpy())
plt.show()

Try now: mr0_upload_seq.ipynb

Simplified Workflows

Utily functions

MRzero now includes streamlined utility functions for common simulation tasks:

1. One line - automatic everything

   # Quickest 2D brain phantom simulation 
   signal, ktraj_adc = mr0.util.simulate(seq)  # Downloads standard phantom automatically
   

2. Two lines - custom phantom parameters

   # Quick 2D brain phantom sim with modifications
   obj = mr0.util.load_phantom([128,128], B0_polynomial=(0,0,-150,1,1,0))  # Custom B0 map
   signal, ktraj_adc = mr0.util.simulate(seq, obj)
   

3. Three lines - full control

   # Custom phantom from URL/file with arbitrary modifications and high accuracy
   obj = mr0.util.load_phantom([128,128], B0_polynomial=(0,0,-150,1,1,0), url='numerical_brain_cropped.mat')  
   obj.B1 = obj.B1**3     # Very inhomogeneous B1+ field
   signal, ktraj_adc = mr0.util.simulate(seq, obj, accuracy=1e-5)
   

Direct PyPulseq Integration

import pypulseq as pp

# Create sequence with PyPulseq
seq = pp.Sequence()
# ... build sequence ...

# Simulate directly - no need to write/read .seq file!
signal, ktraj_adc = mr0.util.simulate(seq)  # seq can be PyPulseq object

For MATLAB Pulseq Users

See MATLAB integration

Examples - Ready-to-run in Google Colab

Notebook	Complexity
Simulate own uploaded seq files	Beginner	mr0_upload_seq.ipynb
FLASH 2D sequence	Beginner	mr0_FLASH_2D_seq.ipynb
GRE EPI 2D sequence	Intermediate	mr0_EPI_2D_seq.ipynb
DWI SE EPI 2D sequence	Advanced	mr0_DWI_SE_EPI.ipynb
TSE 2D sequence	Advanced	mr0_TSE_2D_multi_shot_seq.ipynb

Browse All 20+ Examples →

Further Features

Visualization

# Plot PyPulseq sequence plot with MR-zero signals
seq.plot(plot_now=False)
mr0.util.insert_signal_plot(seq=seq, signal =signal.numpy())
plt.show()

# Plot the detected k-space used by MR-zero
seq.plot_kspace_trajectory()

#Plot the phase distribution graph
graph = mr0.compute_graph(seq0, obj_p, 2000, 1e-4)  # pre-pass
graph.plot()

Simulation

# Simple simulation
signal = mr0.util.simulate(seq)

# With custom phantom
phantom = mr0.util.load_phantom(size=(64, 64))
signal, kspace = mr0.util.simulate(seq, phantom)

Getting Help

Quick help:

• Common issues: FAQ & Troubleshooting
• GitHub Issues: Report bugs or ask for features
• Github Discussions: Ask Questions

MATLAB Pulseq ↔ MR-zero Integration Guide

This guide shows how to seamlessly integrate MATLAB Pulseq with Python MR-zero for powerful sequence development and simulation workflows.

Quick Start for MATLAB Users

Already have MATLAB Pulseq sequences? Get started in 3 steps:

1. Export your sequence: seq.write('my_sequence.seq') in MATLAB
2. Upload to Colab: Use our (mr0_upload_seq.ipynb )
3. Simulate instantly: No Python installation required!

Integration Workflows

Workflow 1: MATLAB seq creation → Python Simulation

Use Case: Leverage MR-zero’s fast simulation while keeping your MATLAB development workflow

% In MATLAB: Create your sequence as usual
seq = mr.Sequence();
% ... build your sequence ...
seq.write('my_flash.seq');

Then switch to python and load the seq file.

# In Python/Colab: Simulate with one line
import MRzeroCore as mr0
seq = mr0.Sequence.import_file('my_flash.seq')
signal = mr0.util.simulate(seq)

You can just upload your seq file in mr0_upload_seq.ipynb

Workflow 2: MATLAB seq creation → Call Python from Matlab using sys command

see Pulseq demo simMR0

Need Help? Try our interactive examples first, or check the GitHub issues for community support.

Pulseq-zero

Frequently asked questions

Common questions and solutions for MRzero Core and Pulseq integration.

What’s the easiest way to simulate a sequence?

A: Code and export your sequence as Pulseq file. One line:

import MRzeroCore as mr0

# Automatic phantom download and simulation
signal, ktraj_adc = mr0.util.simulate('your_sequence.seq')

For complete Pulseq integration features, see: Pulseq Integration Guide →

Can I simulate PyPulseq sequences directly without writing .seq files?

A: Yes! This is new functionality:

import pypulseq as pp
import MRzeroCore as mr0

# Create sequence with PyPulseq
seq = pp.Sequence()
# ... build sequence ...

# Simulate directly - no .seq file needed!
signal, ktraj_adc = mr0.util.simulate(seq)  # seq is PyPulseq object

Note, this is will still generate a temporary seq file. For avoiding seq files completely, check out Pulseq-zero

How do I use custom phantoms with B0 field modifications?

A:

# Load phantom with polynomial B0 augmentation
obj = mr0.util.load_phantom(
    [128, 128], 
    B0_polynomial=(0, 0, -150, 1, 1, 0),  # (const, lin_x, lin_y, quad_x, quad_y, quad_xy)
    url='numerical_brain_cropped.mat'     # Custom phantom URL or local file
)

# Further customize
obj.B1 = obj.B1**3     # Inhomogeneous B1+ field
obj.T2 = 80e-3         # Custom relaxation times

# High-accuracy simulation
signal, ktraj_adc = mr0.util.simulate(seq, obj, accuracy=1e-5)

Does MRzero support ISMRMRD raw data format?

A: Current status: ISMRMRD raw data support is not yet available, but development is underway as a first pull request exists. So it can be expected in upcoming releases. Stay updated:

• Watch the GitHub repository for ISMRMRD support announcements
• Check release notes for new features

My .seq file won’t load / import fails

A: Pulseq 1.4.2 is used and tested by us the most. Pulseq 1.3 and 1.5 should work too, but might need some different load and plot functions.

See also in the comments in mr0_upload_seq.ipynb

Getting version compatibility errors

A: Most issues are resolved by:

• Using matching Pulseq versions between creation and simulation
• Updating to latest MRzeroCore: pip install --upgrade MRzeroCore
• Checking GitHub Issues for known compatibility problems

Simulation is too slow / runs out of memory

A: Reduce sequence resolution, go to single slice, reduce phantom resolution, reduce accuracy.

# Reduce phantom size for testing
phantom = mr0.util.load_phantom(size=(32, 32))  # Instead of (128, 128)

# Reduce simulation accuracy for speed
signal = mr0.util.simulate(seq, accuracy=1e-2)  # Default: 1e-3

# Use smaller sequence for testing
# Remove unnecessary repetitions during development

GPU not being used / CUDA errors

A:

import torch
print(f"CUDA available: {torch.cuda.is_available()}")
print(f"CUDA device: {torch.cuda.get_device_name(0) if torch.cuda.is_available() else 'None'}")

# MR-zero uses GPU when seq0 and obj is using cuda
```python
graph = mr0.compute_graph(seq0, obj, 200, 1e-4)
# Minimum emitted and latent signal of 0.01 (compared per repetition to strongest state)
signal = mr0.execute_graph(graph, seq0.cuda(), obj.cuda(), 0.01, 0.01)

How does MRzero handle slice-selective pulses?

A: Current limitation: MRzero uses instantaneous pulses and currently ignores slice selection gradients. This means:

# In PyPulseq - slice selective pulse
rf, gz, gzr = pp.make_sinc_pulse(
    flip_angle=90*np.pi/180, 
    slice_thickness=5e-3,  # This is ignored in MRzero
    duration=1e-3
)

To have correct gradient behavior (spoiling, diffusion), our interpreter puts the instantanoues pulse in the center of the pulse, and correctly plays out the gradients. Thus, if you have assymetric pulses, you have to be careful, and better replace them by block pulses.

Workarounds:

1. For 2D simulations: Use thin phantom slices to approximate slice selection
2. For comparison with Scans: Consider this limitation when comparing to real measurements, apparent flip angles might be smaller due to slice profile effects.

How do I choose the right accuracy parameter?

A: The accuracy parameter controls simulation precision vs speed. When in doubt rerun with higher accuracy.

# High accuracy (slow, precise)
signal = mr0.util.simulate(seq, accuracy=1e-8)  # Research-quality

# Fast accuracy (quick testing)
signal = mr0.util.simulate(seq, accuracy=1e-2)  # Development/testing

# Very fast (rough estimates)
signal = mr0.util.simulate(seq, accuracy=1e-1)  # Initial debugging only

What does the accuracy parameter actually control?

A: It sets the min_emitted_signal and min_latent_signal thresholds in the PDG simulation:

• Lower values = more computation but higher precision
• Higher values = faster simulation but potential artifacts
• For sequences with small signals (like multi-echo), use lower values

Does MRzero support magnetization transfer effects?

A: Current status: MRzero currently does not model MT effects. The simulation assumes:

• Single pool (free water) only
• No bound pool interactions
• No MT saturation effects

How accurate is diffusion simulation in MRzero?

A: MRzero handles currently only isotropic diffusion, using the EPG equations from Weigel et al.:

How are relaxation times handled?

A:

# Set realistic values in phantom
phantom.T1 = 1.0    # seconds
phantom.T2 = 80e-3  # seconds  
phantom.T2dash = 50e-3  # T2' dephasing effects with 1/T2* = 1/T2 + 1/T2'

# All relaxation effects are included in simulation
# including T2* decay during readout

Tips:

• Use literature values for tissue types
• T2’ does NOT include intravoxel B0 inhomogeneity effects, this can be achieved by high enough B0 map resolution.

Sequence created in MATLAB doesn’t work in Python

A: Check Pulseq version compatibility

# Try different import approaches
seq = mr0.Sequence.import_file('file.seq')  # Standard since Pulseq 1.4, shoudl work for 1.5
seq = mr0.Sequence.from_seq_file('file.seq')  # Alternative for Pulseq 1.3 

Resources:

• GitHub Issues: Report bugs
• Examples: Interactive examples
• API Docs: Complete reference
• Pulseq Integration: Integration guide

Still stuck? Try the interactive examples first - they solve most common issues!

NIfTI phantoms

In previous versions, MR-zero used ad-hoc formats for its simulation phantoms. It stored the right combination of maps in numpys compressed file format. Starting 2026, we fully switch to a new file format with a precise specification. It is a combination of the widely used NIfTI file format an a JSON configuration file. Exact definitions of the contents of this configuration file makes writing loaders for different applications / programming languages easy. At the same time, care was taken to allow great flexibility and human readability of the files.

You can see the specification here. A closer look at the coordinate system is shown here.

To see how to use phantoms, you can either look at examples or view the API documentation.

BrainWeb NIfTI phantoms

Phantoms are built from BrainWeb data and stored in the NIfTI format. This data is not included directly. Instead, a BrainWeb downloader is provided by mr0, which automatically downloads the necessary data and combines it with literature values for the physical properties at different field strenghts. It is provided by `mr0.generate_brainweb_phantoms:

import MRzeroCore as mr0

# Stores the downloaded data in a `cache` folder in the given directory.
# All generated phantoms are emitted as a series of folders.
mr0.generate_brainweb_phantoms("output/brainweb")

These phantoms can be loaded as tissue dictionaries, where individual tissues can be modified before converting it into simulation data. For more information look at the NIfTI API documentation.

tissues = mr0.TissueDict.load("output/brainweb/subj05/subj05-3T.json")
# you can modify individual tissues - its a dictionary of `VoxelGridPhantom`s:
tissues.pop("fat")
tissues["gm"].T1 *= 1.2

# either merge all tissues into one VoxelGridPhantom with `tissues.combine()`,
# which can be useful for plotting the maps or slightly faster simulation
tissues.combine().plot()
# ...or build simulation data directly to get partial volume effects:
sim_data = tissues.build()

Custom phantoms for experimentation

Sometimes it is useful to generate phantoms consisting of individual voxels. This can help with educational tasks, analysing effective resolution and PSFs and more.

For a single voxel:

obj_p = mr0.CustomVoxelPhantom(
    pos=[[-0.25, -0.25, 0]],
    PD=[1.0],
    T1=[3.0],
    T2=[0.5],
    T2dash=[30e-3],
    D=[0.0],
    B0=0,
    voxel_size=0.1,
    voxel_shape="box"  # Can also be "sinc" or "gauss"
)

For a more complex structure try an L-shape:

obj_p = mr0.CustomVoxelPhantom(
    pos=[[-0.4, -0.4, 0], [-0.4, -0.2, 0],
         [-0.3, -0.2, 0], [-0.2, -0.2, 0], [-0.1, -0.2, 0]],
    PD=[1.0, 1.0, 0.5, 0.5, 0.5],
    T1=[1.0, 0.5, 0.5, 0.5, 2],
    T2=0.1,
    T2dash=0.1,
    D=0.0,
    B0=0,
    voxel_size=0.1,
    voxel_shape="box"
)

Predefined phantoms

Older versions of MR-zero Core stored phantoms as MATLAB .mat files or as NumPy .npz. These are deprecated in favour of the new NIfTI phantom standard, but still available online.

A simple 2D brain phantom as MATLAB .mat file

This phantom comes from a quantification measurement. It is 2D, 1Tx (single channel B1+) and contains the following data:

• PD (proton density)
• T1 relaxation
• T2 relaxation
• T2’ dephasing
• ADC (apparent diffusion coefficient)
• dB0 inhomogeneities
• B1+ transmit field
• B1- recieve field

The phantom is stored on GitHub: numerical_brain_cropped.mat. Load it with MR-zero:

# Load file with a shell command - works in interactive notebooks (e.g.: on Google Colab):
!wget https://github.com/MRsources/MRzero-Core/raw/main/documentation/playground_mr0/numerical_brain_cropped.mat &> /dev/null

# After downloading, load with MR-zero
import MRzeroCore as mr0
obj_p = mr0.VoxelGridPhantom.load_mat('numerical_brain_cropped.mat')
# You can modify the phantom first:
# change resolution, remove dB0 inhomogeneities, increase diffusion...
brain_phantom_res = 64
obj_p = obj_p.interpolate(brain_phantom_res, brain_phantom_res, 1)
obj_p.B0[:] = 0
obj_p.D *= 1.5
# Convert the phantom into sparse SimData
obj_p = obj_p.build()

Predefined 3T BrainWeb phantom

The old, pre-generated phantoms from BrainWeb data are replaced by the new NIfTI phantoms. There is still one available for backwards compatibility, stored on GitHub: subject05.npz It contains all maps except for recieve fields:

• PD (proton density)
• T1 relaxation
• T2 relaxation
• T2’ dephasing
• ADC (apparent diffusion coefficient)
• dB0 inhomogeneities (generated empirically)
• B1+ transmit field (generated empirically)
• B1- recieve field

# Load file with a shell command - works in interactive notebooks (e.g.: on Google Colab):
!wget https://github.com/MRsources/MRzero-Core/raw/main/documentation/playground_mr0/subject05.npz &> /dev/null

sz = [64, 64]
# (i) load a phantom object from file
obj_p = mr0.VoxelGridPhantom.brainweb('subject05.npz')
obj_p = obj_p.interpolate(sz[0], sz[1], 32).slices([15])
obj_p.size[2]=0.08

# ========================================================
# EXAMPLE: insert rectangular "Tumor" in the diffusion map
# ========================================================

# typical brain tumor ADC values are around ~1.5 * 10^-3 mm^2/s,
# which lies between GM/WM and CSF (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3000221)
# mr0 uses D in units of 10^-3 * mm ^2/s  this is the same as µm^2/ms

# construct tumor border region
for ii in range(15, 25):
    for jj in range(15, 25):
        obj_p.D[ii, jj] = torch.tensor(0.75)

# construct tumor filling
for ii in range(16, 24):
    for jj in range(16, 24):
        obj_p.D[ii, jj] = torch.tensor(1.5)

# Convert the phantom into sparse SimData
obj_p = obj_p.build()

NIfTI Phantom Specification

The “NIfTI phantom specification” describes the storage of data suitable for MR imaging simulations. A phantom consists of one or more NIfTI files for per-voxel data and a single JSON file that defines the phantom and references the NIfTI files.

1. Storage

The naming of files is a convention: implementations can but are not required to reject non-conforming names. The supported convention is:

📂 subj42
├ 📄 subj42.nii.gz
├ 📄 subj42_dB0.nii.gz
├ 📄 subj42_B1+.nii.gz
├ ...
├ 📄 subj42-3T.json
└ 📄 subj42-7T.json

The phantom name is subj42, used as directory name and file prefix. It is available in two variants: subj42-3T.json and subj42-7T.json. Contents of the JSON file are described in Section 4.

Per-voxel data is stored in <name>_<property>.nii(.gz) files; the density map omits the property postfix. See Section 2 for properties and Section 3 for the NIfTI format requirements. The density property is required for every tissue loaded from NIfTI (as opposed to a constant value).

2. Phantom Definition

Units

The JSON file contains a "units" field that explicitly specifies the unit of each property. Currently, only the default units listed below are supported; automatic conversion may be added in the future.

System

The JSON file contains a "system" field that describes the physical system of the MRI experiment.

The gyro value implicitly defines the nucleus: 42.5764 MHz/T corresponds to ^1H. Tissue properties (particularly relaxation times and dB0) are specific to the system they were measured at. A phantom defined at 3T should not be used in a 7T simulation without adjusting tissue values accordingly.

Tissue Properties

The following properties are defined by the NIfTI phantom specification. If stored as per-voxel NIfTI data, the file postfix should match the key (e.g.: subj42_dB0.nii.gz).

3. NIfTI Data

Per-voxel tissue properties are stored in .nii files following the NIfTI v1.1 specification, optionally gzip-compressed (.nii.gz).

• Each file contains a single property for all tissues
• Data must be 4-dimensional (use singleton dimensions for non-3D data)
  • Dimensions 1-3: spatial (size 1 if unused)
  • Dimension 4: tissue index
• All NIfTI files must share the same resolution and orientation
• Spatial data should follow the RAS+ convention (index 0: R, 1: A, 2: S, growing towards positive) to ensure correct orientation for tools ignoring the affine matrix
• The affine matrix must transform data into RAS+ using mm as units (as per NIfTI spec)

4. JSON Phantom File

The phantom is defined in a .json file with the following structure:

{
  "file_type": "nifti_phantom_v1",
  "units": {
    "gyro": "MHz/T",
    "B0": "T",
    "T1": "s",
    "T2": "s",
    "T2'": "s",
    "ADC": "10^-3 mm^2/s",
    "dB0": "Hz",
    "B1+": "rel",
    "B1-": "rel"
  },
  "system": {
    "gyro": 42.5764,
    "B0": 3.0
  },
  "tissues": {
    "gm": {
      "density": "subj42.nii.gz[0]",
      "T1": 1.56,
      "T2": 0.083,
      "T2'": 0.32,
      "ADC": 0.83,
      "dB0": "subj42_dB0.nii.gz[0]",
      "B1+": [
        "subj42_B1+.nii.gz[0]",
        "subj42_B1+.nii.gz[1]",
        "subj42_B1+.nii.gz[2]",
        "subj42_B1+.nii.gz[3]",
        "subj42_B1+.nii.gz[4]",
        "subj42_B1+.nii.gz[5]",
        "subj42_B1+.nii.gz[6]",
        "subj42_B1+.nii.gz[7]"
      ]
    },
    "wm": {
      "density": "subj42.nii.gz[1]",
      "T1": 0.83,
      "dB0": "subj42_dB0.nii.gz[0]"
    },
    "fat": {
      "density": "subj42.nii.gz[4]",
      "dB0": {
        "file": "subj42_dB0.nii.gz[0]",
        "func": "x - 420"
      }
    }
  }
}

• file_type: must be "nifti_phantom_v1".
• units: units for all properties (see Units).
• system: physical system parameters (see System).
• tissues: keys are arbitrary tissue names; values define tissue properties.

Each tissue must have a density of type file_ref (it defines the tissue’s shape). All other properties are optional and fall back to their default values. B1+ and B1- are lists of maps to support multi-channel data (a single-element list is valid).

Tissue Property Values

Each property in a tissue can be set via one of three methods:

file_ref

A file reference has the form <file_name>[<index>]. The file must be a .nii(.gz) in the same directory as the JSON file. The index selects a volume along the 4th dimension of the NIfTI file, allowing multiple tissues’ data to share a single file.

mapping_func

A mapping function is an equation evaluated per voxel. It supports basic arithmetic (+ - * / ( )) and the following variables:

Examples:

Coordinate System

• We always use RAS+ in a subject-aligned coordinate system
• NIfTI’s can store two orientations at once and do not specify which one to use
• MITK uses a LPS+ coordinate system and negates the xy affine entries on loading
• The scanner says data is in the SCANNER coordinate system, but this changes with sequence settings.
• Phantom z direction should always point in $B_0$ direction

Note

We assume that in measurement and FOV, MRI sequences are aligned to the subject.

When storing phantoms, we should always orient them to the subject-aligned RAS+ system (origin best at center of FOV but can be arbitrary). Correctly stored with sform_code == 2 and qform unused (qform_code == 0), which is the default for nibabel but check for simpleITK!

MR-zero

The MR-zero coordinate system is equivalent to NIfTI’s RAS+ with voxel coordinates at their centers. MR-zero (currently) ignores the affine matrix, except to extract the physical size of the NIfTI data. It generates the voxel coordinates with the following code:

pos_x, pos_y, pos_z = torch.meshgrid(
    size[0] * torch.fft.fftshift(torch.fft.fftfreq(shape[0])),
    size[1] * torch.fft.fftshift(torch.fft.fftfreq(shape[1])),
    size[2] * torch.fft.fftshift(torch.fft.fftfreq(shape[2])),
    indexing="ij"
)

There is no corner vs center here - voxels are points. The grid follows the typical FFT definition torch.fft.fftshift, which means they are placed as follows:

# Even number of samples, e.g.: N = 32
pos = [-16, -15, ..., -1, 0, 1, ..., 15] / 32 * size
# Odd number of samples, e.g.: N = 21
pos = [-10, ..., -1, 0, 1, ..., 10] / 21 * size

This is matching typical FFT impls and k-space encodings: a plain FFT reco will return a pixel-perfect image without half-voxel shift. (This also means that for even sample counts, the phantom extends are not exactly centered. On the other hand, there is always a voxel exactly at (0, 0, 0)).

NIfTI coordinate system

The “true” specification of NIfTI seems to be the C header: nifti1_h.pdf - search for section “3D IMAGE (VOLUME) ORIENTATION AND LOCATION IN SPACE”.

The specification says that NIfTI’s are in RAS+:

the continuous coordinates are referred to as (x,y,z). The voxel index coordinates […] are referred to as (i,j,k)
[…]
The (x,y,z) coordinates refer to the CENTER of a voxel. In methods 2 and 3, the (x,y,z) axes refer to a subject-based coordinate system, with +x = Right +y = Anterior +z = Superior. This is a right-handed coordinate system.
[…]
The i index varies most rapidly, j index next, k index slowest.

The 3 Methods this refers to are:

1. qform_code == 0 && sform_code == 0: compat for ANALYZE files, only supports scaling
2. qform_code > 0: use a quaternion for additional rotation
3. sform_code > 0: full affine matrix specifies rotation, scale, offset (and even skewing)

Problem
NIfTI’s allow qform and sform to co-exist and doesn’t specify which one to use. It states:

In this scheme, a dataset would originally be set up so that the Method 2 coordinates represent what the scanner reported. Later, a registration to some standard space can be computed and inserted in the header. Image display software can use either transform, depending on its purposes and needs.

nibabel defaults to use sform and in normal use never sets different transforms for sform and qform, expect when set manually.

Mapping

The affine matrices map $[i,j,k] \mapsto (x, y, z)$. Even though the standard says that $(x, y, z)$ always refer to the subject-based RAS+ system, the origin of the coordinate system is specified twice, by the qcode_form and the scode_form:

Conclusion

NIfTI’s written with nibabel always use (unless forced) sform_code == 2: affine matrix with arbitrary coordinates, since the scanner system is typically not known when creating phantoms from code. Similarly, when reading files, we usually have no way of transforming between scanner and subject coordinate systems, since this information is not stored in the NIfTI…

Warning

…except if both sform and qform are used and one of them selects scanner coordinates and the other one subject coordinates. We should not use this in our generated phantoms, but have to take care when loading NIfTI files that use this property:

Different viewers could decide differently which of the both systems to use - double check when debugging a wrong orientation!

This information is also found in the NIfTI FAQ - Q17 and Q19

State selection

Documentation on PDG State Parameters and Configuration

Overview

The Phase Distribution Graph (PDG) states are frequently initialized with default parameters. This document outlines the essential considerations and configurations required for optimal performance, especially when dealing with exotic sequences.

Standard PDG Configuration

Typically, PDG states are invoked using the following standard parameters:

graph = mr0.compute_graph(seq0, obj_ph)
signal = mr0.execute_graph(graph, seq0, obj_ph)

The standard values used then are here written explicitly:

graph = mr0.compute_graph(seq0, obj_ph, max_state_count=200, min_state_mag=0.0001)
signal = mr0.execute_graph(graph, seq0, obj_ph, min_emitted_signal=0.01, min_latent_signal=0.01)

Importance of State Selection Thresholds

For arbitrary sequences, careful examination of state selection is crucial. This process remains important for all users employing PDG for scientific purposes.

As outlined in the main paper, there are four main parameters to adjust:

1. Threshold Number of States (max_state_count)
2. Threshold Magnetization (min_state_mag)
3. Threshold Emitted Signal (min_emitted_signal)
4. Threshold Latent Signal (min_latent_signal)

To resolve issues of inaccurate simulations one can adjust the parameters as follows to get a more accurate simulation:

graph = mr0.compute_graph(seq0, obj_ph, max_state_count=5000, min_state_mag=1e-12)
signal = mr0.execute_graph(graph, seq0, obj_ph, min_emitted_signal=0.001, min_latent_signal=0.001)

This adjustment will result in a longer simulation duration, but more accurate results

We recommend periodically validating results with more precise simulations.

Each sequence may require unique adjustments to balance simulation speed and accuracy. It’s often beneficial to start with a higher state count for a more accurate baseline.

For sequences with poor timing or spoiling, state explosion can occur (Hennig 1991). The state selection thresholds will filter out these states, potentially leading to inaccurate simulation results.

API overview

All available functions (snake_case naming) / classes (CamelCase) are listed here. In the following pages, more information is listed. For detailed descriptions of each items, look at their Python docstrings (as written in the source code or shown by your Python IDE) as well.

import MRzeroCore as mr0

• Sequence building blocks
  • mr0.Sequence
  • mr0.chain
  • mr0.Repetition
  • mr0.Pulse
  • mr0.PulseUsage
• Simulation data
  • mr0.SimData
  • mr0.VoxelGridPhantom
  • mr0.TissueDict
  • mr0.CustomVoxelPhantom
• NIfTI phantoms
  • mr0.PhantomUnits
  • mr0.PhantomSystem
  • mr0.NiftiRef
  • mr0.NiftiMapping
  • mr0.NiftiTissue
  • mr0.NiftiPhantom
• Simulation
  • mr0.isochromat_sim
  • mr0.compute_graph
  • mr0.compute_graph_ext
  • mr0.Graph
  • mr0.execute_graph
• Reconstruction
  • mr0.reco_adjoint
• util module
  • mr0.util.get_signal_from_real_system
  • mr0.util.insert_signal_plot
  • mr0.util.pulseq_plot
  • mr0.util.pulseq_plot_142
  • mr0.util.pulseq_plot_pre14
  • mr0.util.imshow
  • mr0.util.load_phantom
  • mr0.util.simulate
  • mr0.util.simulate_2d
• mr0.sig_to_mrd
• mr0.pulseq_write_cartesian
• mr0.generate_brainweb_phantoms

MRI Sequence definitions

The MRzeroCore Sequence is a Python list of Repetitions, each of which is starting with an instantaneous Pulse. To construct sequences, you can

• load from .seq files: Pulseq Integration
• construct with pulseq-zero: Pulseq-zero
• defined in the MR-zero directly:

import MRzeroCore as mr0
import numpy as np

seq = mr0.Sequence()

# Iterate over the repetitions
for i in range(64):
    rep = seq.new_rep(2 + 64 + 2)
    
    # Set the pulse
    rep.pulse.usage = mr0.PulseUsage.EXCIT
    rep.pulse.angle = 5 * np.pi/180
    rep.pulse.phase = i**2 * 117*np.pi/180

    # Set encoding
    rep.gradm[:, :] = ...

    # Set timing
    rep.event_time[:] = ...

    # Set ADC
    rep.adc_usage[2:-2] = 1
    rep.adc_phase[2:-2] = np.pi - rep.pulse.phase

Sequence

Defines an MRI sequence. Derived from a Python list.

Sequence(repetitions, normlized_grads=True)

• Sequence.cpu(): move sequence data to CPU
• Sequence.cuda(device=None): move sequence data to selected or default CUDA device
• Sequence.device (property): the device of the repetition data
• Sequence.clone(): make a copy of the repetition with cloned
• Sequence.new_rep(event_count): appends a Repetition.zero(event_count) to the sequence
• Sequence.get_kspace(): Returns the kspace trajectory of the measured signal as a single tensor (ADC events only, concatenated repetitions)
• Sequence.get_full_kspace(): Returns the full kspace trajectory as a list of per-repetition tensors (containing all events)
• Sequence.plot_kspace_trajectory(): Plot the result of get_full_kspace()
• Sequence.get_contrast_mask(contrast): Returns a bool mask for the selected contrast which can be applied to the simulated signal or kspace
• Sequence.get_contrasts(): Returns a list of all used contrast IDs (apply for all repetitions)
• Sequence.shift_contrasts(offset): Shift all contrast IDs by offset (apply for all repetitions)
• get_duration(): Return the total duration of the sequence in seconds

Sequences can be imported from Pulseq .seq files or Siemens .dsv files:

Sequence.import_file(file_name, exact_trajectories, print_stats, default_shim, ref_voltage, resolution)

chain

Helper function to combine multiple sequences into a multi-contrast sequence.

chain(*sequences, oneshot=False)

Returns a new sequence object. If oneshot is set to:

• True, sequences are simply concatenated.
• False, the adc_usages are shifted so that every sub-sequence has their own set of IDs.

Repetition

Part of a Sequence, containing an RF pulse and a list of event_count events.

Repetition(pulse: Pulse, event_time, gradm, adc_phase, adc_usage)

• Repetition.event_count (property): number of events in this repetition
• Repetition.cpu(): move repetition data to CPU
• Repetition.cuda(device=None): move repetition data to selected or default CUDA device
• Repetition.device (property): the device of the repetition data
• Repetition.zero(event_count): create a zero’d repetition with the given number of events
• Repetition.clone(): make a copy of the repetition with cloned
• Repetition.get_contrasts(): return a list of all used adc_usage values (except zero)
• Repetition.shift_contrasts(offset): add offset to all adc_usage values (except zero)

Pulse

Contains the definition of an instantaneous RF Pulse.

Pulse(usage: PulseUsage, angle, phase, shim_array, selective)

• Pulse.cpu(): move pulse data to CPU
• Pulse.cuda(device=None): move pulse data to selected or default CUDA device
• Pulse.device (property): the device of the pulse data
• Pulse.zero(): alternative constructor which zero-initializes everything
• Pulse.clone(): make a copy of the pulse with cloned

PulseUsage

Used for computing the k-space trajectory or to select the pulse type in legacy exporters.

Simulation Phantoms

Phantoms based on BrainWeb data can be downloaded and generated easily, as explained in Generating Phantoms. The data required by the simulation defined by SimData. It should not be created directly, but rather by defining or loading a VoxelGridPhantom or CustomVoxelPhantom.

Example for loading a BrainWeb phantom, scaling it and using it as SimData on the GPU:

phantom = mr0.VoxelGridPhantom.load("subject05.npz")
phantom = phantom.interpolate(128, 128, 32).slices([16])
data = phantom.build().cuda()

Tip

A new type of phantoms is now officially supported! These are based on NIfTIs and allow partial volume effects, segmented phantoms and easy phantom reconfiguration without code changes. Take a look at their description and API.

Voxel Shape

Simulation data (and therefore all Phantoms) define a shape for their voxels. The following options are available:

Most cases should use "sinc":
Input maps, signals, and reconstructed images are all sampled as points. Extending this bandwidth-limited data to continuous functions means that sinc-shaped voxels are the correct form in basically all cases. Using box-shaped voxels means that a round trip of putting maps into the simulation and then quantifying them from the reconstructed images will inevitably blur them. Box shaped voxels do not have flat signal response like sinc shaped voxels have (up to the nyquist frequency). Point shaped voxels would be a good fit too, except that they would prohibit gradient spoiling (and are a bad description of actual tissue).

SimData

Contains the data used for simulation (together with the Sequence). The values are flattened tensors, which only contain voxels with a non-zero proton density.

# Usually not constructed directly but via the `.build()` function on a phantom
SimData(PD, T1, T2, T2dash, D, B0, B1, coil_sens, size, voxel_pos, nyquist, dephasing_func, recover_func=None, phantom_motion=None, voxel_motion=None, tissue_masks=None)

• SimData.cpu(): move data to the CPU
• SimData.cuda(): move data to the GPU
• SimData.device (property): the device of the repetition data
• SimData.recover(): reconstruct the VoxelGridPhantom or CustomVoxelPhantom that was used to build this SimData

VoxelGridPhantom

Phantom format where tissue parameters are cartesian grids of voxels. The data is stored in 3D or 4D tensors with indices [(channel), x, y, z].

VoxelGridPhantom(PD, T1, T2, T2dash, D, B0, B1, coil_sens, size, phantom_motion=None, voxel_motion=None, tissue_masks=None)

Parameters of this function are identical to those of SimData, except that they are stored in 3D or 4D tensors. They are indexed [(channel), x, y, z].

• VoxelGridPhantom.load(file_name): Loading phantom from a .npz file
• VoxelGridPhantom.load_mat(file_name, ...): Loading phantom from a .mat file
• VoxelGridPhantom.slices(slices): Removing all but the specified list of slices (in z direction)
• VoxelGridPhantom.scale_fft(): Resizing the phantom by doing FFT -> truncate -> FFt
• VoxelGridPhantom.interpolate(): Resizing the phantom by using torch.nn.functional.interpolate(..., mode='trilinear')
• VoxelGridPhantom.plot(): Plotting all maps of the phantom (center slice)
• VoxelGridPhantom.build(PD_threshold): Convert this phantom into SimData, sparsifying the arrays (flattening and masking where pd < PD_threshold)

TissueDict

A tissue dictionary can store multiple, overlapping tissues (which lead to partial volume effects). It is a dictionary of VoxelGridPhantoms and provides additional helper functions.

TissueDict is modelled to hold the data loaded from a NiftiPhantom. To load such a file, you can call TissueDict.load() with the config .json of the NIfTI phantom. Alternatively, you can create a NiftiPhantom and then load the phantom (and the volume data) with TissueDict.load() - this time passing a folder and the NiftiPhantom.

TissueDict(dict[str, VoxelGridPhantom])

• TissueDict.load(path, config=None): either load from .json or from a dir + NiftiPhantom
• TissueDict.save(path, gyro, B0): save as NIfTI phantom
• TissueDict.interpolate(): apply [`VoxelGridPhantom.interpolate()] to all tissues
• TissueDict.slices(): apply [`VoxelGridPhantom.slices()] to all tissues
• TissueDict.combine(): merge all tissues into a single VoxelGridPhantom
• TissueDict.build(PD_threshold, voxel_shape): convert into SimData

CustomVoxelPhantom

Phantom format where tissue parameters are 1D lists of voxels with arbitrary positions.

CustomVoxelPhantom(pos, PD=1.0, T1=1.5, T2=0.1, T2dash=0.05, D=1.0, B0=0.0, B1=1.0, voxel_size=0.1, voxel_shape="sinc")

The size of this phantom is computed from the extends of the voxel positions.
Useful for testing reconstruction, sub-voxel positioning, teaching and more.

Parameters of this function are identical to those of SimData, except that they can also be python lists (and are converted to tensors internally). In addition, it contains a convenience method for plotting the phantom by rendering the voxels (and their selected shape) at any resolution.

• CustomVoxelPhantom.generate_PD_map(): Voxels can be sinc- box- or gauss-shaped and at arbitrary positions. This function renders them into a fixed-resolution proton density map.
• CustomVoxelPhantom.plot(): Render and plot all maps of this phantom
• CustomVoxelPhantom.build(): Convert into SimData instance for simulation

NIfTI Phantoms

The classes shown here reflect the content of the .json config file structure. These files define how various NIfTIs are combined into a phantom. You can load a complete phantom directly with mr0.TissueDict.load. But you can also load only the configuration with mr0.NiftiPhantom.load. This allows you to modify it as needed before loading the phantom or saving the changes.

All classes are dataclasses, augmented with additional functions.

PhantomUnits

Stores the physical units used by the phantom. These are currently fixed to the values used by MR-zero (default values) and are stored for documentation (no conversion is implemented). However, different units and automatic conversion might be implemented in the future.

PhantomUnits(gyro, B0, T1, T2, T2dash, ADC, dB0, B1_tx, B1_rx)

• PhantomUnits.default(): return PhantomUnits with the default units given above
• PhantomUnits.from_dict(config): load from a dictionary and check if units are valid
• PhantomUnits.to_dict(): convert to a dictionary

PhantomSystem

Describes the physical properties of the MRI experiment.

PhantomSystem(gyro, B0)

• PhantomSystem.from_dict(config): load from a dictionary
• PhantomSystem.to_dict(): convert to a dictionary

NiftiRef

Reference to a NIfTI file, including the index to the tissue. NIfTIs can contain multiple 3D volumes, stacked in the 4th dimension. For a description look at the specification.

NiftiRef(file_name, tissue_index)

As given by the spec, NiftiRefs are stored as "path/to/nifti.nii.gz[<index>]" strings:

• NiftiRef.parse(config): parse from a string
• NiftiRef.to_str(): convert to a string

NiftiMapping

Combination of NiftiRef and a mapping function for simple modifications. For a description look at the specification.

NiftiMapping(file, func)

• NiftiMapping.parse(config): parse from a dictionary (calling NiftiRef.parse)
• NiftiMapping.to_dict(): convert to a dictionary

NiftiTissue

Definition of a single tissue. Each scalar property is given by a float, NiftiRef or NiftiMapping. B1_tx and B1_rx are lists to support multi-channel data. Units are given by PhantomUnits.

NiftiTissue(density, T1, T2, T2dash, ADC, dB0, B1_tx, B1_rx)

where Prop = float | NiftiRef | NiftiMapping

• NiftiTissue.default(): return a NiftiTissue with the default units given above - density must be passed
• NiftiTissue.from_dict(config): load from a dictionary
• NiftiTissue.to_dict(): convert to a dictionary

NiftiPhantom

Represents the configuration given by a phantom .json file.

NiftiPhantom(units, system, tissues)

The class attribute file_type is always "nifti_phantom_v1" and is not a constructor parameter.

• NiftiPhantom.default(gyro=42.5764, B0=3): return a NiftiPhantom without tissues
• NiftiPhantom.load(path): load from .json file, given by the path
• NiftiPhantom.save(path): save to a .json file, given by the path
• NiftiPhantom.from_dict(config): load from a dictionary
• NiftiPhantom.to_dict(): convert to a dictionary

Simulation

MR-zero provides two simulations:

• a purposefully simple isochromat based simulation to check for correctness
• a state-of-the-art PDG simulation, built for speed and accuracy.

Both simulations are differentiable and capable to run on any device supported by PyTorch (all CPUs and NVIDIA GPUs). While the isochromat simulation is rather slow and only there for confirmation purposes, PDG is one of the fastest currently available simulations. You can find more information about it in the Paper.

It is strongly adviced to use the PDG simulation - isochromat_sim only exists for testing purposes.

Phase Distribution Graph (PDG) simulation

Phase Distribution Graphs split the magnetization into multiple states, forming a Graph over the duration of the sequence. More details can be found in the PDG Paper, once published.

The simulation runs in two passes:

• The pre-pass (compute_graph) simulates a large number of states, but only roughly estimates the signal. This information is used to assess which states are important for the signal, and which can be skipped because they don’t contribute anything.
• The main-pass (execute_graph) uses the exact signal equation and takes all voxels into account, but simulates as few states as possible. It does so by using the information generated by the pre-pass.

The precision / speed tradeoff for both passes can be chosen as follows:

# Default parameters: max 200 states per repetition
# with at leat 0.0001 of maximum possible magnetization
graph = mr0.compute_graph(seq, data, 200, 1e-4)
# Minimum emitted and latent signal of 0.01 (normalized)
signal = mr0.execute_graph(graph, seq, data, 0.01, 0.01)

The pre-pass limits the state count by both the maximum number (higher = more states / slower + more precise simulaiton) of states and the minimum magnetization (lower = more states). In practice, it turns out that the second parameter is less surprising in how it influences the quality of the simulation, which means that the first parameter should be set to a high number and the second parameter should be increased if the pre-pass is too slow (usually the main-pass is more likely to be the bottleneck).

The main-pass has two parameters which are used as threshold to determine which states are used in the signal equation and which states are simulated. Usually they are both set to the same number. Lower numbers mean that more states are above the thresold and will be simulated. It is not useful to set min_emitted_signal to a value lower than min_latent_signal, as this is trying to include to states into the signal equation that are not simulated.

Examples:

# Very precise simulation
graph = mr0.compute_graph(seq, data, 10000, 1e-7)
signal = mr0.execute_graph(graph, seq, data, 1e-7, 1e-7)

# Average precision
graph = mr0.compute_graph(seq, data, 10000, 1e-4)
signal = mr0.execute_graph(graph, seq, data, 1e-3, 1e-4)

# FID only
graph = mr0.compute_graph(seq, data, 10000, 1e-3)
signal = mr0.execute_graph(graph, seq, data, 1, 0.1)

Optimal settings can always depend on the sequence itself. Some sequences are more demanding and require a very precise simulation, some are already decently described by only a few states.

compute_graph

This utility function computes the parameters for compute_graph_ext from the passed data. Returns the computed Graph.

compute_graph_ext

Compute the Phase Distribution Graph for the given sequence and settings. Runs a simulation of seq for a single voxel with the given properties. Discards states with magnetization below min_state_mag and the states with the lowest magnetization if surpassing max_state_count. Returns the computed Graph.

execute_graph

Takes the Sequence, the SimData and the Graph and calculates the measured ADC signal. Because of the work done by the pre-pass, only the minimal work needed in order to achieve the desired precision is executed. This precision can be tuned by the min_emitted_signal and min_latent_signal thresholds. Higher values lead to less states being simulated, which improves speed and reduces accuracy. A value of 1 will mean that 2 states will be simulated (z0 and one +), resulting in the FID signal. A value of 0 means that everything will be simulated that somehow contributes to the signal, even if vanishingly small.

Graph

A wrapper for the retuned graph, which provides two additional methods:

Graph.plot

Creates a matplotlib scatter plot of all PDG states. Does not create or show the figure, which allows to modify it, for example:

graph = compute_graph(seq, data)

plt.figure()
graph.plot(dephasing="k_x")
plt.grid()
plt.show()

Graph.plot_data

Returns the data used in Graph.plot.

isochromat_sim

Isochromat based Bloch simulation. The code is implemented for readability, not speed. This allows to verify its correctness, bugs are not expected.

Warning

This simulation

• is slow
• does not support diffusion
• always uses "box" - voxels, which is usually a bad choice and will introduce slight blurring. See Voxel Shape for more information. This choice was made for simplicity.

Reconstruction methods

NUFFT

Note

Currently, only one single reconstruction method exists here - despite more (like NUFFT or GRAPPA) being used in many scripts since a long time. We could integrate them here with MR-zero.

Nufft reconstruction is not integrated into mr0 directly, but can be realized with torchkbnufft and code similar to the following:

import torchkbnufft as tkbn

# Construct a sequence "seq" and simulate it, resulting in "signal"

kdata = signal.view(1, 1, -1)
ktraj = seq.get_kspace()[:, :2].T / (2 * np.pi)
im_size = (64, 64)

adjnufft_ob = tkbn.KbNufftAdjoint(im_size)
dcomp = tkbn.calc_density_compensation_function(ktraj, im_size=im_size)
reco = adjnufft_ob(kdata * dcomp, ktraj).view(im_size)

reco_adjoint

Adjoint reconstruction builds a very simple backwards version of the encoding / measurement operation. It essentially resembles a DFT and is thus slower and consumes more memory than an FFT, but can handle any readout trajectory.

util module

Additional helper functions.

import MRzeroCore as mr0
mr0.util.<function>()

get_signal_from_real_system

Wait for scanner TWIX files on the given path, and return if ready. Useful if a shared network folder between scanner and research PC is available.

Note

This is developed for internal use and might not be useful for others.

insert_signal_plot

Insert a measured signal into a currently open pypulseq plot. This is only supported by newer pulseq versions, which do not show (and close) the figure immediately if plot_now=False.

seq.plot(plot_now=False)  # option doesn't exist for older pulseq versions
mr0.util.insert_signal_plot(seq, signal)
plt.show()

pulseq_plot

Wrapper around insert_signal_plot, pulseq_plot_142 and pulseq_plot_pre14 which selects the appropriate plotting function based on the installed pulseq version.

pulseq_plot_142

Use only with pypulseq versions 1.4.x / tested with 1.4.2

Identical to pulseq plot(), but returns the generated axes! This is needed to insert the signal into the ADC plot. Modified from pypulseq 1.4.2 sequence.py - plot() - see pulseq documentation for details.

pulseq_plot_pre14

Use only with pypulseq versions before 1.4

Pulseq sequence plot with some aesthetical changes. Most important, it can insert the simulated signal into the ADC plot. Modified from pypulseq 1.2.0post1 sequence.py - plot() - see pulseq documentation for details.

imshow

Modified version of matplotlib imshow(), tailored for showing MR reconstructions. reco is expected to be an array with dimensions [x, y], [x, y, z] or [coil, x, y, z]

• plots images in RAS+
  • x is left-to-right
  • y is bottom-to-top
• plots 3d data (z) as a side-by-side grid of its slices
• applies quadratic coil-combination for multi-coil data

# Replace
import matplotlib.pyplot as plt
plt.imshow(reco, *args, **kwargs)

# with
import MRzeroCore as mr0
mr0.util.imshow(reco, *args, **kwargs)

load_phantom

Download a phantom from the MR-zero repository (or a custom URL) for getting started quickly.

DEFAULT_PHANTOM_URL = "https://github.com/MRsources/MRzero-Core/raw/main/documentation/playground_mr0/numerical_brain_cropped.mat"

simulate

Helper function to simulate using PDG.

simulate_2d

Caution

Deprecated - use load_phantom + simulate instead. For parameter descriptions, look at these functions.

Other methods

Miscellaneous methods provided by MR-zero. These are currently exported directly, but might be sorted into sub-modules if the list grows.

generate_brainweb_phantoms

Builds MR-zero phantoms with the data provided by BrainWeb. These phantoms are not stored in the repository; call this function once to generate them locally.

The configuration files with the used tissue parameters can be found here: brainweb_data.json

sig_to_mrd

Write the simulation output to an ISMRMD file; supported by many reconstruction tools.

pulseq_write_cartesian

Warning

Outdated - we suggest writing sequences with Pulseq(-zero) instead. Exporters might be removed in the future.

Export an MR-zero sequence to Pulseq. Since MR-zero sequences have no limitations in timing, slew rates, amplitudes and more, this exporter has to try to convert this into a sequence that can actually be measured. The exporter only supports cartesian trajectories due to its use of trapezoidal gradients. There is no guarantee that it produces the expected results.

Differentiable Rounding Helpers for Pulseq Zero Sequences

This module provides differentiable stand-ins for common rounding operations (ceil, floor, round) used in Pulseq sequences. They are designed to keep the optimization pipeline differentiable by passing gradients through unchanged, effectively treating these operations as the identity during backpropagation.

Design idea

• Forward pass: Behaves exactly like the corresponding PyTorch / NumPy rounding operation.
• Backward pass: Returns the incoming gradient unchanged.
• Use case: Allows differentiable optimization of otherwise discrete parameters.

Functional Wrappers

To use the wrapper functions:

from pulseqzero import round, ceil, floor

ceil(x)

Differentiable version of torch.ceil.

y = ceil(x)

floor(x)

Differentiable version of torch.floor.

y = floor(x)

round(x)

Differentiable version of torch.round.

y = round(x)

Literature

This is a list of all papers and conference abstracts (known to us) which reference MR-zero. If you think something is missing, please open an issue or pull request to extend the list.

Papers

Phase distribution graphs for fast, differentiable, and spatially encoded Bloch simulations of arbitrary MRI sequences.
Endres J, Weinmüller S, Dang HN, Zaiss M.
Magn Reson Med. 2024; https://doi.org/10.1002/mrm.30055

MR-zero meets FLASH – Controlling the transient signal decay in gradient- and rf-spoiled gradient echo sequences.
Weinmüller S, Endres J, Dang HN, Stollberger R, Zaiss M.
Magn Reson Med. 2024; https://doi.org/10.1002/mrm.30318

MRzero - Automated discovery of MRI sequences using supervised learning.
Loktyushin A, Herz K, Dang HN.
Magn Reson Med. 2021; https://doi.org/10.1002/mrm.28727

MR-zero meets RARE MRI: Joint optimization of refocusing flip angles and neural networks to minimize T2-induced blurring in spin echo sequences.
Dang, H, Endres, J, Weinmüller S et al.
Magn Reson Med. 2023; https://doi.org/10.1002/mrm.29710

Conference abstracts

Weinmüller S, Dang HN, Endres J, Glang F, Loktyushin A, Zaiss M.
A blurring-free 3D snapshot readout for fast CEST- or relaxation-prepared MRI.
In: 24. Jahrestagung Der Deutschen Sektion Der ISMRM. ; 2022:6-7.

ISMRM

Advances in MRzero: supervised learning of parallel imaging sequences including joint non-Cartesian trajectory and flip angle optimization.
Glang F, Loktyushin A, Herz K, Dang HN, Deshmane A, Weinmüller S, Zaiss M et al.
ISMRM 2021, abstract 4200; https://archive.ismrm.org/2021/4200.html

MRzero sequence generation using analytic signal equations as forward model and neural network reconstruction for efficient auto-encoding.
Weinmüller S, Dang HN, Loktyushin A, Glang F, Doerfler A, Maier A, Schölkopf B, Scheffler K, Zaiss M.
ISMRM 2021, abstract 1761; https://archive.ismrm.org/2021/1761.html

Overcoming System Imperfections Using End-to-End MR Sequence Design.
West D, Glang F, Endres J, Zaiss M, Hajnal J, Malik S.
ISMRM 2023, abstract 1242; https://archive.ismrm.org/2023/1242.html

Non-idealized system (NIS) optimization of EPI sequences at ultra-high field.
West D, Glang F, Endres J, Leitão D, McElroy S, Zaiss M, Hajnal J, Malik S.
ISMRM 2024, abstract 0528; https://archive.ismrm.org/2024/0528.html

ESMRMB 2025

Book of Abstracts ESMRMB 2025 Online 41st Annual Scientific Meeting 8–11 October 2025. https://doi.org/10.1007/s10334-025-01278-8

Breast tissue segmentation at 7T using Chemical Exchange Saturation Transfer (CEST)
Magda Duarte, Katharina Breininger, Katharina Tkotz4, Sebastian Bickelhaupt, Moritz Zaiss
abstract 041

Agentic MR sequence development - Leveraging LLMs with MR tools and tests for physics-informed sequence development
Moritz Zaiss, Jonathan Endres, Simon Weinmüller
abstract 046

Breast Digital Twin: simulating with MR-zero
Magda Duarte, Felix Dietz, Tobias Dornstetter, Jonathan Endres, Simon Weinmüller, Sebastian Bickelhaupt, Moritz Zaiss
abstract 055

End-to-End Optimization of Variable Flip Angle Schemes for Enhanced Apparent Resolution in 7T 3D Fast Spin Echo MRI
Peter Dawood, Martin Blaimer, Angelika Mennecke, Fraticelli Laura, Simon Weinmüller, Jonathan Endres, Peter Jakob, Moritz Zaiss
abstract 066

Highly accelerated snapshot readout for high-resolution CEST imaging
Martin Freudensprung, Patrick Liebig, Simon Weinmüller1, Magda Durate, Moritz Zaiss
abstract 076

Simulation of dynamic B0 with Phase Distribution Graphs
Jonathan Endres, Simon Weinmüller, Moritz Zaiss
abstract 119

Phase graph-based MRI simulation including off-resonant pulse response
Felix Dietz, Simon Weinmüller, Jonathan Endres, Moritz Zaiss
abstract 130

Towards B0 and B1 Robust Spiral Turbo Spin Echo Sequence at 7T
Simon Weinmüller, Jürgen Hennig, Felix Dietz, Peter Dawood, Jonathan Endres, Moritz Zaiss
abstract 233

Patient specific quantitative MR twins for synthetic previews and protocol planning
Deepak Charles, Simon Weinmüller, Fabian Wagner, Rainer Schneider, Jonathan Endres, Moritz Zaiss
abstract 359

Exact motion simulation vs. Retrospective application on static data
Jonathan Endres, Simon Weinmüller, Moritz Zaiss
abstract 445

Optimal control pulses and MIMOSA for CEST preparation at 7 T
Martin Freudensprung, Clemens Stilianu, Simon Weinmüller, Rudolf Stollberger, Moritz Zaiss
abstract 504

ESMRMB 2024

Book of Abstracts ESMRMB 2024 Online 40th Annual Scientific Meeting 2–5 October 2024. https://doi.org/10.1007/s10334-024-01191-6

Quantalizer: a quantitative digital twin to create and run tailored sequences.
Simon Weinmüller1, Jonathan Endres, Hoai Nam Dang, Junaid Rajput, Felix Glang, Andreas Maier, Moritz Zaiss
abstract 011

MR-zero-tailored dummy refocusing pulses in TSE sequences reduce artifacts.
Jacob Malich, Simon Weinmüller, Jonathan Endres, Peter Dawood, Moritz Zaiss
abstract 112

kT-points optimization using Pulseq and MR-zero at 7T.
Martin Freudensprung, Simon Weinmüller1, Jonathan Endres,,Armin Nagel, Moritz Zaiss
abstract 422

Joint reconstruction of EPI and FLASH using Phase Distribution Graphs as a general MRI reconstruction model.
Magda Duarte, Jonathan Endres, Junaid Rajput, Simon Weinmüller1, Moritz Zaiss
abstract 456

ESMRMB 2023

Book of Abstracts ESMRMB 2023 Online 39th Annual Scientific Meeting 4–7 October 2023. https://doi.org/10.1007/s10334-023-01108-9

Dynamic DREAM MRI: B0, B1 and Tx/Rx-phase mapping for assisting motion tracking systems.
Tim Baum, Simon Weinmüller, Armin Nagel, Martin Vossiek, Moritz Zaiss
LT52

A simple pTx Pulseq extension for pulse-specific B1 shimming.
Martin Freudensprung, Simon Weinmüller1, Jonathan Endres, Armin Nagel, Moritz Zaiss
LT53

A magnetization-prepared DREAM sequence for CEST imaging with an intrinsic dynamic B0 and B1 reference (CEST-MP-DREAM).
Tim Baum, Simon Weinmüller, Jonathan Endres, Patrick Liebig, Moritz Zaiss
LT68

DREAM-zero – Optimized variable flip angles for decreased image blurring in magnetization-prepared DREAM sequences.
Simon Weinmüller, Tim Baum, Hoai Nam Dang, Jonathan Endres, Moritz Zaiss
P198