qmri.dro

Digital Reference Objects for quantitative MRI validation.

Module Overview

The qmri.dro module provides tools for generating synthetic MRI data with known ground truth parameters.

Types

GroundTruth dataclass

GroundTruth(value: T, units: str, description: str)

Bases: Generic[T]

Container for a ground truth parameter value.

Attributes:

Name	Type	Description
value	T	The ground truth parameter value.
units	str	Physical units of the parameter.
description	str	Human-readable description of the parameter.

Example

from qmri.dro import GroundTruth

adc_gt = GroundTruth(
    value=1e-3,
    units="mm²/s",
    description="Apparent Diffusion Coefficient",
)

DWIPhantom dataclass

DWIPhantom(
    signal: NDArray[floating],
    b_values: NDArray[floating],
    ground_truth: dict[
        str, GroundTruth[float | NDArray[floating]]
    ],
    snr: float | None,
    seed: int | None,
)

Digital Reference Object for diffusion-weighted imaging.

Contains synthetic DWI signal with known ground truth ADC values, along with acquisition parameters and optional noise characteristics.

Attributes:

Name	Type	Description
signal	NDArray[floating]	Signal intensities at each b-value. Shape is either (n_bvalues,) for single voxel or (..., n_bvalues) for multi-voxel.
b_values	NDArray[floating]	Diffusion weighting values in s/mm².
ground_truth	dict[str, GroundTruth[float | NDArray[floating]]]	Dictionary of ground truth parameters.
snr	float | None	Signal-to-noise ratio used for noise generation, or None if noiseless.
seed	int | None	Random seed used for reproducibility, or None if not seeded.

Example

from qmri.dro import dwi

phantom = dwi.generate(adc=1e-3, snr=50, seed=42)
print(phantom.ground_truth["adc"].value)  # 0.001
print(phantom.signal.shape)  # (4,) for default b-values

T1Phantom dataclass

T1Phantom(
    signal: NDArray[floating],
    time_points: NDArray[floating],
    method: str,
    model: str | None,
    repetition_time: float | None,
    ground_truth: dict[
        str, GroundTruth[float | NDArray[floating]]
    ],
    snr: float | None,
    seed: int | None,
)

Digital Reference Object for T1 relaxometry.

Contains synthetic T1 signal with known ground truth values, supporting both inversion recovery (IR) and variable TR (VTR) methods.

Attributes:

Name	Type	Description
signal	NDArray[floating]	Signal intensities at each time point. Shape is either (n_timepoints,) for single voxel or (..., n_timepoints) for multi-voxel.
time_points	NDArray[floating]	Time points in seconds (TI for IR, TR for VTR).
method	str	Acquisition method used ("ir" or "vtr").
model	str | None	For IR method, the model used ("general" or "classical").
repetition_time	float | None	TR value(s) for IR method, or None for VTR.
ground_truth	dict[str, GroundTruth[float | NDArray[floating]]]	Dictionary of ground truth parameters.
snr	float | None	Signal-to-noise ratio used for noise generation, or None if noiseless.
seed	int | None	Random seed used for reproducibility, or None if not seeded.

Example

from qmri.dro import relaxometry

phantom = relaxometry.generate_t1_ir(
    t1=1.2,
    inversion_times=[0.1, 0.5, 1.0, 2.0],
    repetition_time=5.0,
    snr=100,
    seed=42,
)
print(phantom.ground_truth["t1"].value)  # 1.2

ASLPhantom dataclass

ASLPhantom(
    control: NDArray[floating],
    label: NDArray[floating],
    m0: NDArray[floating],
    ground_truth: dict[
        str, GroundTruth[float | NDArray[floating]]
    ],
    acquisition_params: dict[str, float],
    snr: float | None,
    seed: int | None,
)

Digital Reference Object for arterial spin labelling.

Contains synthetic ASL control and label images with known ground truth perfusion values, along with acquisition parameters.

Attributes:

Name	Type	Description
control	NDArray[floating]	Control image signal intensity.
label	NDArray[floating]	Label image signal intensity.
m0	NDArray[floating]	Equilibrium magnetisation (M0) image.
ground_truth	dict[str, GroundTruth[float | NDArray[floating]]]	Dictionary of ground truth parameters.
acquisition_params	dict[str, float]	Dictionary of acquisition parameters used.
snr	float | None	Signal-to-noise ratio used for noise generation, or None if noiseless.
seed	int | None	Random seed used for reproducibility, or None if not seeded.

Example

from qmri.dro import perfusion

phantom = perfusion.generate_pcasl(
    perfusion_rate=60.0,
    m0=1000.0,
    snr=50,
    seed=42,
)
delta_m = phantom.control - phantom.label
print(phantom.ground_truth["perfusion_rate"].value)  # 60.0

DWI Phantom Generation

dwi

DWI phantom generation for ADC validation.

This module provides functions for generating synthetic diffusion-weighted imaging (DWI) data with known ground truth ADC values.

Example

from qmri.dro import dwi
from qmri.diffusion import adc

# Generate a single-voxel phantom
phantom = dwi.generate(adc=1e-3, s0=1000, snr=50, seed=42)

# Fit and compare to ground truth
result = adc.fit(phantom.signal, phantom.b_values)
print(f"True ADC: {phantom.ground_truth['adc'].value:.2e}")
print(f"Fitted ADC: {result.adc:.2e}")

generate

generate(
    adc: float,
    s0: float = ...,
    b_values: Sequence[float] = ...,
    *,
    snr: float | None = ...,
    noise_model: Literal["rician", "gaussian"] = ...,
    seed: int | None = ...,
) -> DWIPhantom

generate(
    adc: NDArray[floating],
    s0: NDArray[floating] | float = ...,
    b_values: Sequence[float] = ...,
    *,
    snr: float | None = ...,
    noise_model: Literal["rician", "gaussian"] = ...,
    seed: int | None = ...,
) -> DWIPhantom

generate(
    adc: NDArray[floating] | float,
    s0: NDArray[floating] | float = 1000.0,
    b_values: Sequence[float] = (0, 500, 1000, 2000),
    *,
    snr: float | None = None,
    noise_model: Literal["rician", "gaussian"] = "rician",
    seed: int | None = None,
) -> DWIPhantom

Generate a DWI phantom with known ADC.

Creates synthetic DWI signal using the mono-exponential diffusion model with optional noise for validation and testing purposes.

Parameters:

Name	Type	Description	Default
adc	NDArray[floating] | float	Apparent Diffusion Coefficient in mm²/s. Can be a scalar for single-voxel phantom or an array for multi-voxel phantom.	required
s0	NDArray[floating] | float	Baseline signal intensity (at b=0). Can be scalar or array matching the shape of adc. Default is 1000.0.	1000.0
b_values	Sequence[float]	Diffusion weighting values in s/mm². Default is (0, 500, 1000, 2000).	(0, 500, 1000, 2000)
snr	float | None	Signal-to-noise ratio. If None, no noise is added.	None
noise_model	Literal['rician', 'gaussian']	Type of noise to add ("rician" or "gaussian"). Default is "rician", which is more realistic for magnitude MRI.	'rician'
seed	int | None	Random seed for reproducibility. If None, uses random state.	None

Returns:

Type	Description
DWIPhantom	DWIPhantom containing the synthetic signal, b-values, and ground truth.

Raises:

Type	Description
ValueError	If noise_model is not "rician" or "gaussian".

Example

Single-voxel phantom:

from qmri.dro import dwi

phantom = dwi.generate(adc=1e-3, snr=50, seed=42)
print(phantom.signal)
# array([1000.  , 606.5..., 367.8..., 135.3...])

Multi-voxel phantom:

import numpy as np
from qmri.dro import dwi

adc_map = np.array([[0.5e-3, 1.0e-3],
                   [1.5e-3, 2.0e-3]])
phantom = dwi.generate(adc=adc_map, snr=100, seed=42)
print(phantom.signal.shape)
# (2, 2, 4)

Notes

The signal is generated using the Stejskal-Tanner equation:

\[S(b) = S_0 \\exp(-b \\cdot \\text{ADC})\]

For magnitude MRI, Rician noise is the appropriate model and causes a positive bias at low SNR. Gaussian noise is suitable for complex data or high SNR scenarios.

Source code in packages/qmri-dro/src/qmri/dro/dwi.py

def generate(
    adc: NDArray[np.floating] | float,
    s0: NDArray[np.floating] | float = 1000.0,
    b_values: Sequence[float] = (0, 500, 1000, 2000),
    *,
    snr: float | None = None,
    noise_model: Literal["rician", "gaussian"] = "rician",
    seed: int | None = None,
) -> DWIPhantom:
    r"""Generate a DWI phantom with known ADC.

    Creates synthetic DWI signal using the mono-exponential diffusion model
    with optional noise for validation and testing purposes.

    Args:
        adc: Apparent Diffusion Coefficient in mm²/s. Can be a scalar for
            single-voxel phantom or an array for multi-voxel phantom.
        s0: Baseline signal intensity (at b=0). Can be scalar or array
            matching the shape of `adc`. Default is 1000.0.
        b_values: Diffusion weighting values in s/mm². Default is
            (0, 500, 1000, 2000).
        snr: Signal-to-noise ratio. If None, no noise is added.
        noise_model: Type of noise to add ("rician" or "gaussian").
            Default is "rician", which is more realistic for magnitude MRI.
        seed: Random seed for reproducibility. If None, uses random state.

    Returns:
        DWIPhantom containing the synthetic signal, b-values, and ground truth.

    Raises:
        ValueError: If noise_model is not "rician" or "gaussian".

    Example:
        Single-voxel phantom:

        ```python
        from qmri.dro import dwi

        phantom = dwi.generate(adc=1e-3, snr=50, seed=42)
        print(phantom.signal)
        # array([1000.  , 606.5..., 367.8..., 135.3...])
        ```

        Multi-voxel phantom:

        ```python
        import numpy as np
        from qmri.dro import dwi

        adc_map = np.array([[0.5e-3, 1.0e-3],
                           [1.5e-3, 2.0e-3]])
        phantom = dwi.generate(adc=adc_map, snr=100, seed=42)
        print(phantom.signal.shape)
        # (2, 2, 4)
        ```

    Notes:
        The signal is generated using the Stejskal-Tanner equation:

        $$S(b) = S_0 \\exp(-b \\cdot \\text{ADC})$$

        For magnitude MRI, Rician noise is the appropriate model and causes
        a positive bias at low SNR. Gaussian noise is suitable for complex
        data or high SNR scenarios.
    """
    if noise_model not in ("rician", "gaussian"):
        msg = f"noise_model must be 'rician' or 'gaussian', got '{noise_model}'"
        raise ValueError(msg)

    b_arr = np.asarray(b_values, dtype=np.float64)
    adc_arr = np.asarray(adc, dtype=np.float64)
    s0_arr = np.asarray(s0, dtype=np.float64)

    # Generate clean signal
    # For multi-voxel, we need to handle broadcasting carefully
    if adc_arr.ndim == 0:
        # Single voxel case
        signal = adc_module.signal_model(s0_arr, adc_arr, b_arr)
    else:
        # Multi-voxel case: output shape is (*adc.shape, n_bvalues)
        # Broadcast s0 to match adc shape
        s0_bc = np.broadcast_to(s0_arr, adc_arr.shape)

        # Generate signal for each b-value
        signal = np.zeros((*adc_arr.shape, len(b_arr)), dtype=np.float64)
        for i, b in enumerate(b_arr):
            signal[..., i] = s0_bc * np.exp(-b * adc_arr)

    # Add noise if requested
    if snr is not None:
        rng = np.random.default_rng(seed)
        if noise_model == "rician":
            signal = rician.add_noise(signal, snr=snr, rng=rng)
        else:
            signal = gaussian.add_noise(signal, snr=snr, rng=rng)

    # Build ground truth dictionary
    ground_truth: dict[str, GroundTruth[float | NDArray[np.floating]]] = {
        "adc": GroundTruth(
            value=float(adc) if adc_arr.ndim == 0 else adc_arr,
            units="mm²/s",
            description="Apparent Diffusion Coefficient",
        ),
        "s0": GroundTruth(
            value=float(s0) if np.asarray(s0).ndim == 0 else s0_arr,
            units="a.u.",
            description="Baseline signal intensity (b=0)",
        ),
    }

    return DWIPhantom(
        signal=signal,
        b_values=b_arr,
        ground_truth=ground_truth,
        snr=snr,
        seed=seed,
    )

generate_calibration_phantom

generate_calibration_phantom(
    adc_values: Sequence[float] = (
        0.0003,
        0.0007,
        0.001,
        0.0015,
        0.002,
        0.003,
    ),
    b_values: Sequence[float] = (
        0,
        50,
        100,
        200,
        400,
        600,
        800,
        1000,
    ),
    *,
    s0: float = 1000.0,
    snr: float | None = 50.0,
    noise_model: Literal["rician", "gaussian"] = "rician",
    seed: int | None = None,
) -> DWIPhantom

Generate a calibration phantom with multiple ADC values.

Creates a phantom with multiple voxels at different ADC values, suitable for method validation and benchmarking.

Parameters:

Name	Type	Description	Default
adc_values	Sequence[float]	Sequence of ADC values in mm²/s. Default covers the range from restricted diffusion to free water.	(0.0003, 0.0007, 0.001, 0.0015, 0.002, 0.003)
b_values	Sequence[float]	Diffusion weighting values in s/mm². Default provides good sampling for ADC range.	(0, 50, 100, 200, 400, 600, 800, 1000)
s0	float	Baseline signal intensity for all voxels. Default is 1000.0.	1000.0
snr	float | None	Signal-to-noise ratio. Default is 50.0.	50.0
noise_model	Literal['rician', 'gaussian']	Type of noise to add. Default is "rician".	'rician'
seed	int | None	Random seed for reproducibility.	None

Returns:

Type	Description
DWIPhantom	DWIPhantom with shape (n_adc_values, n_bvalues).

Example

from qmri.dro import dwi
from qmri.diffusion import adc

# Generate calibration phantom
phantom = dwi.generate_calibration_phantom(seed=42)
print(f"Shape: {phantom.signal.shape}")
# Shape: (6, 8)

# Fit each "voxel"
for i, true_adc in enumerate(phantom.ground_truth['adc'].value):
    result = adc.fit(phantom.signal[i], phantom.b_values)
    print(f"True: {true_adc:.2e}, Fitted: {result.adc:.2e}")

Notes

Default ADC values represent:

• 0.3e-3: Highly restricted (e.g., tumour)
• 0.7e-3: White matter
• 1.0e-3: Grey matter
• 1.5e-3: Less restricted tissue
• 2.0e-3: CSF-like
• 3.0e-3: Free water

Source code in packages/qmri-dro/src/qmri/dro/dwi.py

def generate_calibration_phantom(
    adc_values: Sequence[float] = (0.3e-3, 0.7e-3, 1.0e-3, 1.5e-3, 2.0e-3, 3.0e-3),
    b_values: Sequence[float] = (0, 50, 100, 200, 400, 600, 800, 1000),
    *,
    s0: float = 1000.0,
    snr: float | None = 50.0,
    noise_model: Literal["rician", "gaussian"] = "rician",
    seed: int | None = None,
) -> DWIPhantom:
    """Generate a calibration phantom with multiple ADC values.

    Creates a phantom with multiple voxels at different ADC values,
    suitable for method validation and benchmarking.

    Args:
        adc_values: Sequence of ADC values in mm²/s. Default covers the
            range from restricted diffusion to free water.
        b_values: Diffusion weighting values in s/mm². Default provides
            good sampling for ADC range.
        s0: Baseline signal intensity for all voxels. Default is 1000.0.
        snr: Signal-to-noise ratio. Default is 50.0.
        noise_model: Type of noise to add. Default is "rician".
        seed: Random seed for reproducibility.

    Returns:
        DWIPhantom with shape (n_adc_values, n_bvalues).

    Example:
        ```python
        from qmri.dro import dwi
        from qmri.diffusion import adc

        # Generate calibration phantom
        phantom = dwi.generate_calibration_phantom(seed=42)
        print(f"Shape: {phantom.signal.shape}")
        # Shape: (6, 8)

        # Fit each "voxel"
        for i, true_adc in enumerate(phantom.ground_truth['adc'].value):
            result = adc.fit(phantom.signal[i], phantom.b_values)
            print(f"True: {true_adc:.2e}, Fitted: {result.adc:.2e}")
        ```

    Notes:
        Default ADC values represent:

        - 0.3e-3: Highly restricted (e.g., tumour)
        - 0.7e-3: White matter
        - 1.0e-3: Grey matter
        - 1.5e-3: Less restricted tissue
        - 2.0e-3: CSF-like
        - 3.0e-3: Free water
    """
    adc_arr = np.asarray(adc_values, dtype=np.float64)

    return generate(
        adc=adc_arr,
        s0=s0,
        b_values=b_values,
        snr=snr,
        noise_model=noise_model,
        seed=seed,
    )

Relaxometry Phantom Generation

relaxometry

Relaxometry phantom generation for T1/T2 validation.

This module provides functions for generating synthetic T1 relaxometry data with known ground truth values using inversion recovery (IR) and variable TR (VTR) methods.

Example

from qmri.dro import relaxometry
from qmri.relaxometry import t1

# Generate IR data with known T1
phantom = relaxometry.generate_t1_ir(
    t1=1.2,
    inversion_times=[0.1, 0.5, 1.0, 2.0, 3.0],
    repetition_time=5.0,
    snr=100,
    seed=42,
)

# Fit and validate
result = t1.fit_ir(phantom.signal, phantom.time_points, repetition_times=5.0)
print(f"True T1: {phantom.ground_truth['t1'].value:.2f} s")
print(f"Fitted T1: {float(result.t1):.2f} s")

generate_t1_ir

generate_t1_ir(
    t1: float,
    inversion_times: Sequence[float],
    *,
    s0: float = ...,
    repetition_time: float = ...,
    inversion_efficiency: float = ...,
    model: Literal["general", "classical"] = ...,
    snr: float | None = ...,
    noise_model: Literal["rician", "gaussian"] = ...,
    seed: int | None = ...,
) -> T1Phantom

generate_t1_ir(
    t1: NDArray[floating],
    inversion_times: Sequence[float],
    *,
    s0: NDArray[floating] | float = ...,
    repetition_time: float = ...,
    inversion_efficiency: NDArray[floating] | float = ...,
    model: Literal["general", "classical"] = ...,
    snr: float | None = ...,
    noise_model: Literal["rician", "gaussian"] = ...,
    seed: int | None = ...,
) -> T1Phantom

generate_t1_ir(
    t1: NDArray[floating] | float,
    inversion_times: Sequence[float],
    *,
    s0: NDArray[floating] | float = 1000.0,
    repetition_time: float = 5.0,
    inversion_efficiency: NDArray[floating] | float = 1.0,
    model: Literal["general", "classical"] = "general",
    snr: float | None = None,
    noise_model: Literal["rician", "gaussian"] = "rician",
    seed: int | None = None,
) -> T1Phantom

Generate a T1 phantom using inversion recovery.

Creates synthetic IR signal using either the general or classical model with optional noise for validation and testing.

Parameters:

Name	Type	Description	Default
t1	NDArray[floating] | float	T1 relaxation time in seconds. Can be scalar or array.	required
inversion_times	Sequence[float]	Inversion times (TI) in seconds.	required
s0	NDArray[floating] | float	Signal amplitude at equilibrium. Default is 1000.0.	1000.0
repetition_time	float	Repetition time (TR) in seconds. Default is 5.0.	5.0
inversion_efficiency	NDArray[floating] | float	Inversion efficiency (0-1). Default is 1.0.	1.0
model	Literal['general', 'classical']	IR model to use. Default is "general". - "general": Full model including TR recovery term - "classical": Assumes TR >> T1	'general'
snr	float | None	Signal-to-noise ratio. If None, no noise is added.	None
noise_model	Literal['rician', 'gaussian']	Type of noise to add. Default is "rician".	'rician'
seed	int | None	Random seed for reproducibility.	None

Returns:

Type	Description
T1Phantom	T1Phantom containing the synthetic signal, time points, and ground truth.

Raises:

Type	Description
ValueError	If model is not "general" or "classical".
ValueError	If noise_model is not "rician" or "gaussian".

Example

Single-voxel phantom:

from qmri.dro import relaxometry

phantom = relaxometry.generate_t1_ir(
    t1=1.2,
    inversion_times=[0.1, 0.5, 1.0, 2.0, 3.0],
    snr=100,
    seed=42,
)
print(phantom.signal)

Multi-voxel phantom:

import numpy as np
from qmri.dro import relaxometry

t1_map = np.array([[0.5, 1.0],
                  [1.5, 2.0]])
phantom = relaxometry.generate_t1_ir(
    t1=t1_map,
    inversion_times=[0.1, 0.5, 1.0, 2.0],
    snr=50,
    seed=42,
)
print(phantom.signal.shape)
# (2, 2, 4)

Notes

General model:

\[S = S_0 \\left(1 - 2 \\alpha \\exp\\left(-\\frac{TI}{T_1}\\right) + \\exp\\left(-\\frac{TR}{T_1}\\right)\\right)\]

Classical model (assumes TR >> T1):

\[S = S_0 \\left(1 - 2 \\alpha \\exp\\left(-\\frac{TI}{T_1}\\right)\\right)\]

Source code in packages/qmri-dro/src/qmri/dro/relaxometry.py

def generate_t1_ir(
    t1: NDArray[np.floating] | float,
    inversion_times: Sequence[float],
    *,
    s0: NDArray[np.floating] | float = 1000.0,
    repetition_time: float = 5.0,
    inversion_efficiency: NDArray[np.floating] | float = 1.0,
    model: Literal["general", "classical"] = "general",
    snr: float | None = None,
    noise_model: Literal["rician", "gaussian"] = "rician",
    seed: int | None = None,
) -> T1Phantom:
    r"""Generate a T1 phantom using inversion recovery.

    Creates synthetic IR signal using either the general or classical
    model with optional noise for validation and testing.

    Args:
        t1: T1 relaxation time in seconds. Can be scalar or array.
        inversion_times: Inversion times (TI) in seconds.
        s0: Signal amplitude at equilibrium. Default is 1000.0.
        repetition_time: Repetition time (TR) in seconds. Default is 5.0.
        inversion_efficiency: Inversion efficiency (0-1). Default is 1.0.
        model: IR model to use. Default is "general".
            - "general": Full model including TR recovery term
            - "classical": Assumes TR >> T1
        snr: Signal-to-noise ratio. If None, no noise is added.
        noise_model: Type of noise to add. Default is "rician".
        seed: Random seed for reproducibility.

    Returns:
        T1Phantom containing the synthetic signal, time points, and ground truth.

    Raises:
        ValueError: If model is not "general" or "classical".
        ValueError: If noise_model is not "rician" or "gaussian".

    Example:
        Single-voxel phantom:

        ```python
        from qmri.dro import relaxometry

        phantom = relaxometry.generate_t1_ir(
            t1=1.2,
            inversion_times=[0.1, 0.5, 1.0, 2.0, 3.0],
            snr=100,
            seed=42,
        )
        print(phantom.signal)
        ```

        Multi-voxel phantom:

        ```python
        import numpy as np
        from qmri.dro import relaxometry

        t1_map = np.array([[0.5, 1.0],
                          [1.5, 2.0]])
        phantom = relaxometry.generate_t1_ir(
            t1=t1_map,
            inversion_times=[0.1, 0.5, 1.0, 2.0],
            snr=50,
            seed=42,
        )
        print(phantom.signal.shape)
        # (2, 2, 4)
        ```

    Notes:
        **General model**:

        $$S = S_0 \\left(1 - 2 \\alpha \\exp\\left(-\\frac{TI}{T_1}\\right)
        + \\exp\\left(-\\frac{TR}{T_1}\\right)\\right)$$

        **Classical model** (assumes TR >> T1):

        $$S = S_0 \\left(1 - 2 \\alpha \\exp\\left(-\\frac{TI}{T_1}\\right)\\right)$$
    """
    if model not in ("general", "classical"):
        msg = f"model must be 'general' or 'classical', got '{model}'"
        raise ValueError(msg)
    if noise_model not in ("rician", "gaussian"):
        msg = f"noise_model must be 'rician' or 'gaussian', got '{noise_model}'"
        raise ValueError(msg)

    ti_arr = np.asarray(inversion_times, dtype=np.float64)
    t1_arr = np.asarray(t1, dtype=np.float64)
    s0_arr = np.asarray(s0, dtype=np.float64)
    alpha_arr = np.asarray(inversion_efficiency, dtype=np.float64)

    # Generate clean signal
    if t1_arr.ndim == 0:
        # Single voxel case
        if model == "general":
            signal = t1_module.signal_ir(
                s0=s0_arr,
                t1=t1_arr,
                inversion_times=ti_arr,
                repetition_times=repetition_time,
                inversion_efficiency=alpha_arr,
            )
        else:
            signal = t1_module.signal_ir_classical(
                s0=s0_arr,
                t1=t1_arr,
                inversion_times=ti_arr,
                inversion_efficiency=alpha_arr,
            )
    else:
        # Multi-voxel case: output shape is (*t1.shape, n_timepoints)
        s0_bc = np.broadcast_to(s0_arr, t1_arr.shape)
        alpha_bc = np.broadcast_to(alpha_arr, t1_arr.shape)

        signal = np.zeros((*t1_arr.shape, len(ti_arr)), dtype=np.float64)
        for i, ti in enumerate(ti_arr):
            if model == "general":
                signal[..., i] = s0_bc * (
                    1
                    - 2 * alpha_bc * np.exp(-ti / t1_arr)
                    + np.exp(-repetition_time / t1_arr)
                )
            else:
                signal[..., i] = s0_bc * (1 - 2 * alpha_bc * np.exp(-ti / t1_arr))

    # Add noise if requested
    if snr is not None:
        rng = np.random.default_rng(seed)
        if noise_model == "rician":
            signal = rician.add_noise(signal, snr=snr, rng=rng)
        else:
            signal = gaussian.add_noise(signal, snr=snr, rng=rng)

    # Build ground truth dictionary
    ground_truth: dict[str, GroundTruth[float | NDArray[np.floating]]] = {
        "t1": GroundTruth(
            value=float(t1) if t1_arr.ndim == 0 else t1_arr,
            units="s",
            description="T1 relaxation time",
        ),
        "s0": GroundTruth(
            value=float(s0) if np.asarray(s0).ndim == 0 else s0_arr,
            units="a.u.",
            description="Signal amplitude at equilibrium",
        ),
        "inversion_efficiency": GroundTruth(
            value=(
                float(inversion_efficiency)
                if np.asarray(inversion_efficiency).ndim == 0
                else alpha_arr
            ),
            units="",
            description="Inversion efficiency",
        ),
    }

    return T1Phantom(
        signal=signal,
        time_points=ti_arr,
        method="ir",
        model=model,
        repetition_time=repetition_time,
        ground_truth=ground_truth,
        snr=snr,
        seed=seed,
    )

generate_t1_vtr

generate_t1_vtr(
    t1: float,
    repetition_times: Sequence[float],
    *,
    m: float = ...,
    snr: float | None = ...,
    noise_model: Literal["rician", "gaussian"] = ...,
    seed: int | None = ...,
) -> T1Phantom

generate_t1_vtr(
    t1: NDArray[floating],
    repetition_times: Sequence[float],
    *,
    m: NDArray[floating] | float = ...,
    snr: float | None = ...,
    noise_model: Literal["rician", "gaussian"] = ...,
    seed: int | None = ...,
) -> T1Phantom

generate_t1_vtr(
    t1: NDArray[floating] | float,
    repetition_times: Sequence[float],
    *,
    m: NDArray[floating] | float = 1000.0,
    snr: float | None = None,
    noise_model: Literal["rician", "gaussian"] = "rician",
    seed: int | None = None,
) -> T1Phantom

Generate a T1 phantom using variable TR method.

Creates synthetic VTR signal for saturation recovery T1 mapping with optional noise for validation and testing.

Parameters:

Name	Type	Description	Default
t1	NDArray[floating] | float	T1 relaxation time in seconds. Can be scalar or array.	required
repetition_times	Sequence[float]	Repetition times (TR) in seconds.	required
m	NDArray[floating] | float	Equilibrium magnetisation. Default is 1000.0.	1000.0
snr	float | None	Signal-to-noise ratio. If None, no noise is added.	None
noise_model	Literal['rician', 'gaussian']	Type of noise to add. Default is "rician".	'rician'
seed	int | None	Random seed for reproducibility.	None

Returns:

Type	Description
T1Phantom	T1Phantom containing the synthetic signal, time points, and ground truth.

Raises:

Type	Description
ValueError	If noise_model is not "rician" or "gaussian".

Example

from qmri.dro import relaxometry
from qmri.relaxometry import t1

phantom = relaxometry.generate_t1_vtr(
    t1=1.2,
    repetition_times=[0.5, 1.0, 2.0, 4.0, 8.0],
    snr=100,
    seed=42,
)

result = t1.fit_vtr(phantom.signal, phantom.time_points)
print(f"True T1: {phantom.ground_truth['t1'].value:.2f} s")
print(f"Fitted T1: {float(result.t1):.2f} s")

Notes

The VTR signal model is:

\[S = M \\left(1 - \\exp\\left(-\\frac{TR}{T_1}\\right)\\right)\]

where M is the equilibrium magnetisation.

Source code in packages/qmri-dro/src/qmri/dro/relaxometry.py

def generate_t1_vtr(
    t1: NDArray[np.floating] | float,
    repetition_times: Sequence[float],
    *,
    m: NDArray[np.floating] | float = 1000.0,
    snr: float | None = None,
    noise_model: Literal["rician", "gaussian"] = "rician",
    seed: int | None = None,
) -> T1Phantom:
    r"""Generate a T1 phantom using variable TR method.

    Creates synthetic VTR signal for saturation recovery T1 mapping
    with optional noise for validation and testing.

    Args:
        t1: T1 relaxation time in seconds. Can be scalar or array.
        repetition_times: Repetition times (TR) in seconds.
        m: Equilibrium magnetisation. Default is 1000.0.
        snr: Signal-to-noise ratio. If None, no noise is added.
        noise_model: Type of noise to add. Default is "rician".
        seed: Random seed for reproducibility.

    Returns:
        T1Phantom containing the synthetic signal, time points, and ground truth.

    Raises:
        ValueError: If noise_model is not "rician" or "gaussian".

    Example:
        ```python
        from qmri.dro import relaxometry
        from qmri.relaxometry import t1

        phantom = relaxometry.generate_t1_vtr(
            t1=1.2,
            repetition_times=[0.5, 1.0, 2.0, 4.0, 8.0],
            snr=100,
            seed=42,
        )

        result = t1.fit_vtr(phantom.signal, phantom.time_points)
        print(f"True T1: {phantom.ground_truth['t1'].value:.2f} s")
        print(f"Fitted T1: {float(result.t1):.2f} s")
        ```

    Notes:
        The VTR signal model is:

        $$S = M \\left(1 - \\exp\\left(-\\frac{TR}{T_1}\\right)\\right)$$

        where M is the equilibrium magnetisation.
    """
    if noise_model not in ("rician", "gaussian"):
        msg = f"noise_model must be 'rician' or 'gaussian', got '{noise_model}'"
        raise ValueError(msg)

    tr_arr = np.asarray(repetition_times, dtype=np.float64)
    t1_arr = np.asarray(t1, dtype=np.float64)
    m_arr = np.asarray(m, dtype=np.float64)

    # Generate clean signal
    if t1_arr.ndim == 0:
        # Single voxel case
        signal = t1_module.signal_vtr(m=m_arr, t1=t1_arr, repetition_times=tr_arr)
    else:
        # Multi-voxel case: output shape is (*t1.shape, n_timepoints)
        m_bc = np.broadcast_to(m_arr, t1_arr.shape)

        signal = np.zeros((*t1_arr.shape, len(tr_arr)), dtype=np.float64)
        for i, tr in enumerate(tr_arr):
            signal[..., i] = m_bc * (1 - np.exp(-tr / t1_arr))

    # Add noise if requested
    if snr is not None:
        rng = np.random.default_rng(seed)
        if noise_model == "rician":
            signal = rician.add_noise(signal, snr=snr, rng=rng)
        else:
            signal = gaussian.add_noise(signal, snr=snr, rng=rng)

    # Build ground truth dictionary
    ground_truth: dict[str, GroundTruth[float | NDArray[np.floating]]] = {
        "t1": GroundTruth(
            value=float(t1) if t1_arr.ndim == 0 else t1_arr,
            units="s",
            description="T1 relaxation time",
        ),
        "m": GroundTruth(
            value=float(m) if np.asarray(m).ndim == 0 else m_arr,
            units="a.u.",
            description="Equilibrium magnetisation",
        ),
    }

    return T1Phantom(
        signal=signal,
        time_points=tr_arr,
        method="vtr",
        model=None,
        repetition_time=None,
        ground_truth=ground_truth,
        snr=snr,
        seed=seed,
    )

Perfusion Phantom Generation

perfusion

Perfusion phantom generation for ASL validation.

This module provides functions for generating synthetic arterial spin labelling (ASL) data with known ground truth perfusion values.

Example

from qmri.dro import perfusion

# Generate pCASL data with known perfusion
phantom = perfusion.generate_pcasl(
    perfusion_rate=60.0,  # ml/100g/min
    m0=1000.0,
    transit_time=1.0,
    snr=50,
    seed=42,
)

# Access the difference signal
delta_m = phantom.control - phantom.label
print(f"True CBF: {phantom.ground_truth['perfusion_rate'].value} ml/100g/min")

generate_pcasl

generate_pcasl(
    perfusion_rate: float,
    m0: float,
    *,
    transit_time: float = ...,
    label_duration: float = ...,
    post_label_delay: float = ...,
    label_efficiency: float = ...,
    partition_coefficient: float = ...,
    t1_blood: float = ...,
    t1_tissue: float = ...,
    snr: float | None = ...,
    noise_model: Literal["rician", "gaussian"] = ...,
    seed: int | None = ...,
) -> ASLPhantom

generate_pcasl(
    perfusion_rate: NDArray[floating],
    m0: NDArray[floating] | float,
    *,
    transit_time: NDArray[floating] | float = ...,
    label_duration: float = ...,
    post_label_delay: float = ...,
    label_efficiency: float = ...,
    partition_coefficient: NDArray[floating] | float = ...,
    t1_blood: float = ...,
    t1_tissue: NDArray[floating] | float = ...,
    snr: float | None = ...,
    noise_model: Literal["rician", "gaussian"] = ...,
    seed: int | None = ...,
) -> ASLPhantom

generate_pcasl(
    perfusion_rate: NDArray[floating] | float,
    m0: NDArray[floating] | float,
    *,
    transit_time: NDArray[floating] | float = 1.0,
    label_duration: float = 1.8,
    post_label_delay: float = 1.8,
    label_efficiency: float = 0.85,
    partition_coefficient: NDArray[floating] | float = 0.9,
    t1_blood: float = 1.65,
    t1_tissue: NDArray[floating] | float = 1.3,
    snr: float | None = None,
    noise_model: Literal["rician", "gaussian"] = "rician",
    seed: int | None = None,
) -> ASLPhantom

Generate a pCASL phantom with known perfusion.

Creates synthetic pseudo-continuous ASL (pCASL) control and label images using the General Kinetic Model with optional noise.

Parameters:

Name	Type	Description	Default
perfusion_rate	NDArray[floating] | float	Cerebral blood flow (CBF) in ml/100g/min. Can be scalar or array for multi-voxel phantom.	required
m0	NDArray[floating] | float	Equilibrium magnetisation (M0) of tissue.	required
transit_time	NDArray[floating] | float	Arterial transit time (ATT) in seconds. Default is 1.0.	1.0
label_duration	float	Duration of labelling pulse (tau) in seconds. Default is 1.8.	1.8
post_label_delay	float	Post-label delay (PLD) in seconds. Default is 1.8.	1.8
label_efficiency	float	Labelling efficiency (0-1). Default is 0.85.	0.85
partition_coefficient	NDArray[floating] | float	Blood-brain partition coefficient (lambda) in ml/g. Default is 0.9.	0.9
t1_blood	float	T1 of arterial blood in seconds. Default is 1.65.	1.65
t1_tissue	NDArray[floating] | float	T1 of tissue in seconds. Default is 1.3.	1.3
snr	float | None	Signal-to-noise ratio. If None, no noise is added.	None
noise_model	Literal['rician', 'gaussian']	Type of noise to add. Default is "rician".	'rician'
seed	int | None	Random seed for reproducibility.	None

Returns:

Type	Description
ASLPhantom	ASLPhantom containing control, label, and M0 images with ground truth.

Raises:

Type	Description
ValueError	If noise_model is not "rician" or "gaussian".

Example

Single-voxel phantom:

from qmri.dro import perfusion

phantom = perfusion.generate_pcasl(
    perfusion_rate=60.0,
    m0=1000.0,
    snr=50,
    seed=42,
)
print(f"Control: {phantom.control}")
print(f"Label: {phantom.label}")
print(f"Delta M: {phantom.control - phantom.label}")

Multi-voxel phantom:

import numpy as np
from qmri.dro import perfusion

cbf_map = np.array([[40, 60],
                   [80, 100]])
phantom = perfusion.generate_pcasl(
    perfusion_rate=cbf_map,
    m0=1000.0,
    snr=50,
    seed=42,
)
print(phantom.control.shape)  # (2, 2)

Notes

The signal is calculated using the General Kinetic Model (GKM):

\[\\Delta M(t) = 2 M_{0,b} f T_1' \\alpha e^{-\\delta/T_{1,b}} (1 - e^{-\\tau/T_1'}) e^{-(t-\\tau-\\delta)/T_1'}\]

Default parameters are based on ASL White Paper recommendations for adult brain imaging at 3T.

Typical perfusion values:

• Grey matter: 50-80 ml/100g/min
• White matter: 20-30 ml/100g/min
• Tumour: Variable, often elevated

Source code in packages/qmri-dro/src/qmri/dro/perfusion.py

def generate_pcasl(
    perfusion_rate: NDArray[np.floating] | float,
    m0: NDArray[np.floating] | float,
    *,
    transit_time: NDArray[np.floating] | float = 1.0,
    label_duration: float = 1.8,
    post_label_delay: float = 1.8,
    label_efficiency: float = 0.85,
    partition_coefficient: NDArray[np.floating] | float = 0.9,
    t1_blood: float = 1.65,
    t1_tissue: NDArray[np.floating] | float = 1.3,
    snr: float | None = None,
    noise_model: Literal["rician", "gaussian"] = "rician",
    seed: int | None = None,
) -> ASLPhantom:
    r"""Generate a pCASL phantom with known perfusion.

    Creates synthetic pseudo-continuous ASL (pCASL) control and label
    images using the General Kinetic Model with optional noise.

    Args:
        perfusion_rate: Cerebral blood flow (CBF) in ml/100g/min.
            Can be scalar or array for multi-voxel phantom.
        m0: Equilibrium magnetisation (M0) of tissue.
        transit_time: Arterial transit time (ATT) in seconds. Default is 1.0.
        label_duration: Duration of labelling pulse (tau) in seconds.
            Default is 1.8.
        post_label_delay: Post-label delay (PLD) in seconds. Default is 1.8.
        label_efficiency: Labelling efficiency (0-1). Default is 0.85.
        partition_coefficient: Blood-brain partition coefficient (lambda)
            in ml/g. Default is 0.9.
        t1_blood: T1 of arterial blood in seconds. Default is 1.65.
        t1_tissue: T1 of tissue in seconds. Default is 1.3.
        snr: Signal-to-noise ratio. If None, no noise is added.
        noise_model: Type of noise to add. Default is "rician".
        seed: Random seed for reproducibility.

    Returns:
        ASLPhantom containing control, label, and M0 images with ground truth.

    Raises:
        ValueError: If noise_model is not "rician" or "gaussian".

    Example:
        Single-voxel phantom:

        ```python
        from qmri.dro import perfusion

        phantom = perfusion.generate_pcasl(
            perfusion_rate=60.0,
            m0=1000.0,
            snr=50,
            seed=42,
        )
        print(f"Control: {phantom.control}")
        print(f"Label: {phantom.label}")
        print(f"Delta M: {phantom.control - phantom.label}")
        ```

        Multi-voxel phantom:

        ```python
        import numpy as np
        from qmri.dro import perfusion

        cbf_map = np.array([[40, 60],
                           [80, 100]])
        phantom = perfusion.generate_pcasl(
            perfusion_rate=cbf_map,
            m0=1000.0,
            snr=50,
            seed=42,
        )
        print(phantom.control.shape)  # (2, 2)
        ```

    Notes:
        The signal is calculated using the General Kinetic Model (GKM):

        $$\\Delta M(t) = 2 M_{0,b} f T_1' \\alpha e^{-\\delta/T_{1,b}}
        (1 - e^{-\\tau/T_1'}) e^{-(t-\\tau-\\delta)/T_1'}$$

        Default parameters are based on ASL White Paper recommendations
        for adult brain imaging at 3T.

        Typical perfusion values:

        - Grey matter: 50-80 ml/100g/min
        - White matter: 20-30 ml/100g/min
        - Tumour: Variable, often elevated
    """
    if noise_model not in ("rician", "gaussian"):
        msg = f"noise_model must be 'rician' or 'gaussian', got '{noise_model}'"
        raise ValueError(msg)

    cbf_arr = np.asarray(perfusion_rate, dtype=np.float64)
    m0_arr = np.asarray(m0, dtype=np.float64)
    att_arr = np.asarray(transit_time, dtype=np.float64)
    lam_arr = np.asarray(partition_coefficient, dtype=np.float64)
    t1_t_arr = np.asarray(t1_tissue, dtype=np.float64)

    # Calculate signal time (time after labelling start)
    signal_time = label_duration + post_label_delay

    # Calculate delta_m using GKM
    gkm_result = gkm.signal_gkm(
        perfusion_rate=cbf_arr,
        transit_time=att_arr,
        m0_tissue=m0_arr,
        label_duration=label_duration,
        signal_time=signal_time,
        label_efficiency=label_efficiency,
        partition_coefficient=lam_arr,
        t1_blood=t1_blood,
        t1_tissue=t1_t_arr,
        label_type="pcasl",
    )

    delta_m = gkm_result.delta_m

    # Control signal is just M0 (equilibrium magnetisation)
    # In reality it would be M0 * signal decay, but for simplicity
    # we use the M0 value as the control baseline
    control: NDArray[np.floating]
    label: NDArray[np.floating]
    if cbf_arr.ndim == 0:
        control = np.asarray(m0_arr, dtype=np.float64)
        label = np.asarray(m0_arr - delta_m, dtype=np.float64)
        m0_out = np.asarray(m0_arr, dtype=np.float64)
    else:
        control = np.broadcast_to(m0_arr, cbf_arr.shape).copy().astype(np.float64)
        label = (control - delta_m).astype(np.float64)
        m0_out = control.copy()

    # Add noise if requested
    if snr is not None:
        rng = np.random.default_rng(seed)
        if noise_model == "rician":
            control = rician.add_noise(control, snr=snr, rng=rng)
            label = rician.add_noise(label, snr=snr, rng=rng)
        else:
            control = gaussian.add_noise(control, snr=snr, rng=rng)
            label = gaussian.add_noise(label, snr=snr, rng=rng)

    # Build ground truth dictionary
    ground_truth: dict[str, GroundTruth[float | NDArray[np.floating]]] = {
        "perfusion_rate": GroundTruth(
            value=float(perfusion_rate) if cbf_arr.ndim == 0 else cbf_arr,
            units="ml/100g/min",
            description="Cerebral blood flow (CBF)",
        ),
        "transit_time": GroundTruth(
            value=float(transit_time) if att_arr.ndim == 0 else att_arr,
            units="s",
            description="Arterial transit time (ATT)",
        ),
        "m0": GroundTruth(
            value=float(m0) if np.asarray(m0).ndim == 0 else m0_arr,
            units="a.u.",
            description="Equilibrium magnetisation",
        ),
        "t1_tissue": GroundTruth(
            value=float(t1_tissue) if t1_t_arr.ndim == 0 else t1_t_arr,
            units="s",
            description="T1 of tissue",
        ),
    }

    # Store acquisition parameters
    acquisition_params = {
        "label_duration": label_duration,
        "post_label_delay": post_label_delay,
        "label_efficiency": label_efficiency,
        "partition_coefficient": (
            float(partition_coefficient)
            if np.asarray(partition_coefficient).ndim == 0
            else float(np.mean(lam_arr))
        ),
        "t1_blood": t1_blood,
    }

    return ASLPhantom(
        control=control,
        label=label,
        m0=m0_out,
        ground_truth=ground_truth,
        acquisition_params=acquisition_params,
        snr=snr,
        seed=seed,
    )