qmri.thermometry¶

thermometry ¶

MR thermometry models.

This module provides functions for:

Proton Resonance Frequency (PRF) thermometry
Multi-echo dual-resonance thermometry (ethylene glycol phantoms)

R_SQUARED_THRESHOLD `module-attribute` ¶

R_SQUARED_THRESHOLD: float = 0.9

Default threshold for acceptable fit quality (R² > 0.9).

RANDOM_SEED `module-attribute` ¶

RANDOM_SEED: int = 840275920

Fixed seed for reproducibility of bootstrap sampling.

DfInitMethod `module-attribute` ¶

DfInitMethod = Literal["multistart", "fixed", "lombscargle"]

Strategy for choosing the frequency-difference starting value of the fit.

"multistart" (default): fit from both the fixed default and the data-driven Lomb-Scargle estimate, and keep the highest-R² result. Most robust against frequency aliasing.
"fixed": a single fit from the fixed default starting value (:data:_DEFAULT_DF_GUESS). Cheapest; can alias on cold phantoms.
"lombscargle": a single fit seeded from the Lomb-Scargle estimate, falling back to the fixed default when no estimate can be made.

RegionAnalysisMethod `module-attribute` ¶

RegionAnalysisMethod = Literal[
    "regionwise", "voxelwise", "regionwise_bootstrap"
]

Analysis methods for image-based (segmentation-driven) thermometry fitting.

MultiEchoResult `dataclass` ¶

MultiEchoResult(
    temperature: float,
    temperature_uncertainty: float,
    df: float,
    r_squared: float,
    fitted_params: NDArray[floating],
    n_bootstrap: int | None = None,
)

Result of multi-echo dual-resonance thermometry fitting.

Attributes:

Name	Type	Description
`temperature`	`float`	Estimated temperature in degrees Celsius.
`temperature_uncertainty`	`float`	Uncertainty in temperature (°C). For single fits, from covariance matrix. For bootstrap fits, from standard deviation of bootstrap samples.
`df`	`float`	Fitted frequency difference in Hz.
`r_squared`	`float`	Coefficient of determination (R²) of the fit. For bootstrap, this is the mean R² across samples.
`fitted_params`	`NDArray[floating]`	Fitted parameters [A1, A2, R21, R22, df, dphi_deg]. For bootstrap, these are the mean parameters.
`n_bootstrap`	`int \| None`	Number of bootstrap samples (None for single fit).

RegionThermometryResult `dataclass` ¶

RegionThermometryResult(
    region_id: int,
    region_size: int,
    temperature: float,
    temperature_uncertainty: float,
    coverage_factor: float,
    temperature_values: NDArray[floating],
    temperature_uncertainty_values: NDArray[floating],
    r_squared: NDArray[floating],
    fitted_params: NDArray[floating],
    mean_fitted_params: NDArray[floating],
    signal_values: NDArray[floating],
)

Per-region results from segmentation-driven multi-echo thermometry.

A "region" is the set of voxels sharing a single non-zero integer label in the segmentation image. The interpretation of the per-fit arrays depends on the analysis method:

regionwise: a single fit of the region-mean signal (arrays length 1).
voxelwise: one fit per voxel (arrays length region_size).
regionwise_bootstrap: one fit per bootstrap sample (arrays length n_bootstrap).

Attributes:

Name	Type	Description
`region_id`	`int`	The integer label of the region in the segmentation.
`region_size`	`int`	Number of voxels in the region.
`temperature`	`float`	Representative region temperature in °C. For `voxelwise` this is the inverse-variance weighted mean of voxel temperatures; for the region methods it is the temperature of the (mean-signal) fit, or the mean across bootstrap samples passing the R² threshold.
`temperature_uncertainty`	`float`	Standard uncertainty (coverage factor `coverage_factor`) of `temperature` in °C.
`coverage_factor`	`float`	Coverage factor k for `temperature_uncertainty` (k=1).
`temperature_values`	`NDArray[floating]`	Temperature estimate from each individual fit (°C).
`temperature_uncertainty_values`	`NDArray[floating]`	Per-fit standard uncertainty in °C, derived from the fitted Δf covariance.
`r_squared`	`NDArray[floating]`	Coefficient of determination R² for each individual fit.
`fitted_params`	`NDArray[floating]`	Fitted parameters per fit, shape `(n_fits, 6)`.
`mean_fitted_params`	`NDArray[floating]`	Mean of `fitted_params` over fits passing the R² threshold, shape `(6,)` (NaN if no fit passes).
`signal_values`	`NDArray[floating]`	Signal values fed to each fit, shape `(n_fits, n_echoes)`.

to_dict ¶

to_dict() -> dict[str, object]

Return a JSON-serialisable dictionary of the region results.

Source code in packages/qmri/src/qmri/thermometry/multiecho.py

def to_dict(self) -> dict[str, object]:
    """Return a JSON-serialisable dictionary of the region results."""
    return {
        "id": self.region_id,
        "temperature": self.temperature,
        "temperature_uncertainty": [
            self.temperature_uncertainty,
            self.coverage_factor,
        ],
        "region_size": self.region_size,
        "mean_fitted_params": self.mean_fitted_params.tolist(),
        "region_temperature_values": self.temperature_values.tolist(),
        "region_temperature_uncertainty_values": (
            self.temperature_uncertainty_values.tolist()
        ),
        "fitted_params": self.fitted_params.tolist(),
        "signal_values": self.signal_values.tolist(),
        "r_squared": self.r_squared.tolist(),
    }

PRFResult `dataclass` ¶

PRFResult(
    temperature_change: NDArray[floating],
    phase_difference: NDArray[floating],
)

Result of PRF thermometry calculation.

Attributes:

Name	Type	Description
`temperature_change`	`NDArray[floating]`	Temperature change in degrees Celsius.
`phase_difference`	`NDArray[floating]`	Phase difference between heated and baseline images in radians.

calculate_df_from_temperature ¶

calculate_df_from_temperature(
    temperature_celsius: float | NDArray[floating],
    magnetic_field_tesla: float,
) -> float | NDArray[floating]

Calculate frequency difference from temperature for ethylene glycol.

Uses the empirical relationship:

\[\Delta f = \frac{(193.35 - T) \cdot \gamma \cdot B_0}{1.02 \times 10^8}\]

This relationship is specific to ethylene glycol and should not be used for other substances without recalibration.

Parameters:

Name	Type	Description	Default
`temperature_celsius`	`float \| NDArray[floating]`	Temperature in degrees Celsius.	required
`magnetic_field_tesla`	`float`	Magnetic field strength in Tesla.	required

Returns:

Type	Description
`float \| NDArray[floating]`	Frequency difference in Hz.

Example

from qmri.thermometry.multiecho import calculate_df_from_temperature

# At 37°C and 3T
df = calculate_df_from_temperature(37.0, 3.0)
print(f"Frequency difference: {df:.1f} Hz")

Source code in packages/qmri/src/qmri/thermometry/multiecho.py

def calculate_df_from_temperature(
    temperature_celsius: float | NDArray[np.floating],
    magnetic_field_tesla: float,
) -> float | NDArray[np.floating]:
    r"""Calculate frequency difference from temperature for ethylene glycol.

    Uses the empirical relationship:

    $$\Delta f = \frac{(193.35 - T) \cdot \gamma \cdot B_0}{1.02 \times 10^8}$$

    This relationship is specific to ethylene glycol and should not be
    used for other substances without recalibration.

    Args:
        temperature_celsius: Temperature in degrees Celsius.
        magnetic_field_tesla: Magnetic field strength in Tesla.

    Returns:
        Frequency difference in Hz.

    Example:
        ```python
        from qmri.thermometry.multiecho import calculate_df_from_temperature

        # At 37°C and 3T
        df = calculate_df_from_temperature(37.0, 3.0)
        print(f"Frequency difference: {df:.1f} Hz")
        ```
    """
    return ((193.35 - temperature_celsius) * GAMMA_H * magnetic_field_tesla) / 1.02e8

calculate_temperature_from_df ¶

calculate_temperature_from_df(
    df: float | NDArray[floating],
    magnetic_field_tesla: float,
) -> float | NDArray[floating]

Calculate temperature from frequency difference for ethylene glycol.

Uses the empirical relationship:

\[T[°C] = 193.35 - \frac{1.02 \times 10^8 \cdot |\Delta f|}{\gamma \cdot B_0}\]

This relationship is specific to ethylene glycol and should not be used for other substances without recalibration.

Parameters:

Name	Type	Description	Default
`df`	`float \| NDArray[floating]`	Frequency difference in Hz.	required
`magnetic_field_tesla`	`float`	Magnetic field strength in Tesla.	required

Returns:

Type	Description
`float \| NDArray[floating]`	Temperature in degrees Celsius.

Example

from qmri.thermometry.multiecho import calculate_temperature_from_df

# Convert frequency difference to temperature at 3T
temperature = calculate_temperature_from_df(200.0, 3.0)
print(f"Temperature: {temperature:.1f} °C")

Source code in packages/qmri/src/qmri/thermometry/multiecho.py

def calculate_temperature_from_df(
    df: float | NDArray[np.floating],
    magnetic_field_tesla: float,
) -> float | NDArray[np.floating]:
    r"""Calculate temperature from frequency difference for ethylene glycol.

    Uses the empirical relationship:

    $$T[°C] = 193.35 - \frac{1.02 \times 10^8 \cdot |\Delta f|}{\gamma \cdot B_0}$$

    This relationship is specific to ethylene glycol and should not be
    used for other substances without recalibration.

    Args:
        df: Frequency difference in Hz.
        magnetic_field_tesla: Magnetic field strength in Tesla.

    Returns:
        Temperature in degrees Celsius.

    Example:
        ```python
        from qmri.thermometry.multiecho import calculate_temperature_from_df

        # Convert frequency difference to temperature at 3T
        temperature = calculate_temperature_from_df(200.0, 3.0)
        print(f"Temperature: {temperature:.1f} °C")
        ```
    """
    return 193.35 - (1.02e8 * np.abs(df)) / (GAMMA_H * magnetic_field_tesla)

calculate_temperature_uncertainty ¶

calculate_temperature_uncertainty(
    df_uncertainty: float, magnetic_field_tesla: float
) -> float

Calculate temperature uncertainty from frequency difference uncertainty.

Uses uncertainty propagation:

\[u(T) = \frac{1.02 \times 10^8 \cdot u(\Delta f)}{\gamma \cdot B_0}\]

Parameters:

Name	Type	Description	Default
`df_uncertainty`	`float`	Uncertainty in frequency difference (Hz).	required
`magnetic_field_tesla`	`float`	Magnetic field strength in Tesla.	required

Returns:

Type	Description
`float`	Uncertainty in temperature (°C).

Example

from qmri.thermometry.multiecho import calculate_temperature_uncertainty

# 1 Hz uncertainty in frequency at 3T
temp_uncertainty = calculate_temperature_uncertainty(1.0, 3.0)
print(f"Temperature uncertainty: {temp_uncertainty:.2f} °C")

Source code in packages/qmri/src/qmri/thermometry/multiecho.py

def calculate_temperature_uncertainty(
    df_uncertainty: float,
    magnetic_field_tesla: float,
) -> float:
    r"""Calculate temperature uncertainty from frequency difference uncertainty.

    Uses uncertainty propagation:

    $$u(T) = \frac{1.02 \times 10^8 \cdot u(\Delta f)}{\gamma \cdot B_0}$$

    Args:
        df_uncertainty: Uncertainty in frequency difference (Hz).
        magnetic_field_tesla: Magnetic field strength in Tesla.

    Returns:
        Uncertainty in temperature (°C).

    Example:
        ```python
        from qmri.thermometry.multiecho import calculate_temperature_uncertainty

        # 1 Hz uncertainty in frequency at 3T
        temp_uncertainty = calculate_temperature_uncertainty(1.0, 3.0)
        print(f"Temperature uncertainty: {temp_uncertainty:.2f} °C")
        ```
    """
    return (1.02e8 * df_uncertainty) / (GAMMA_H * magnetic_field_tesla)

fit_multiecho_thermometry ¶

fit_multiecho_thermometry(
    signal: NDArray[floating],
    echo_times: NDArray[floating],
    magnetic_field_tesla: float,
    method: Literal["single", "bootstrap"] = "single",
    n_bootstrap: int = 100,
    r_squared_threshold: float = R_SQUARED_THRESHOLD,
    df_init: DfInitMethod = "multistart",
) -> MultiEchoResult

Fit multi-echo signal to dual-resonance model for thermometry.

This function fits the dual-resonance signal model to multi-echo magnitude data and converts the fitted frequency difference to temperature using the ethylene glycol calibration.

Parameters:

Name	Type	Description	Default
`signal`	`NDArray[floating]`	Multi-echo magnitude signal array. Shape should be (n_echoes,).	required
`echo_times`	`NDArray[floating]`	Array of echo times in seconds.	required
`magnetic_field_tesla`	`float`	Magnetic field strength in Tesla.	required
`method`	`Literal['single', 'bootstrap']`	Fitting method. Options: - "single": Single least-squares fit (default). - "bootstrap": Bootstrap resampling for uncertainty estimation.	`'single'`
`n_bootstrap`	`int`	Number of bootstrap samples (default 100). Only used when method="bootstrap".	`100`
`r_squared_threshold`	`float`	Minimum R² for accepting a fit (default 0.9). For bootstrap, samples below threshold are excluded.	`R_SQUARED_THRESHOLD`
`df_init`	`DfInitMethod`	Frequency starting-value strategy — `"multistart"` (default), `"fixed"` or `"lombscargle"`. See :data:`DfInitMethod`.	`'multistart'`

Returns:

Type	Description
`MultiEchoResult`	MultiEchoResult containing temperature, uncertainty, and fit parameters.

Raises:

Type	Description
`ValueError`	If signal and echo_times have different lengths.

Example

import numpy as np
from qmri.thermometry.multiecho import (
    thermometry_signal_model,
    calculate_df_from_temperature,
    fit_multiecho_thermometry,
)

# Generate synthetic data at 25°C
magnetic_field = 3.0
echo_times = np.linspace(0.001, 0.024, 24)
df_true = calculate_df_from_temperature(25.0, magnetic_field)
signal = thermometry_signal_model(
    echo_times, 1.0, 0.5, 50.0, 100.0, df_true, 45.0
)

# Fit the model
result = fit_multiecho_thermometry(
    signal, echo_times, magnetic_field, method="single"
)
temp = result.temperature
uncert = result.temperature_uncertainty
print(f"Temperature: {temp:.1f} ± {uncert:.2f} °C")
print(f"R²: {result.r_squared:.4f}")

# With bootstrap uncertainty
result_bs = fit_multiecho_thermometry(
    signal, echo_times, magnetic_field,
    method="bootstrap", n_bootstrap=50
)
print(f"Bootstrap uncertainty: {result_bs.temperature_uncertainty:.2f} °C")

Note

The temperature-frequency calibration is specific to ethylene glycol. For other substances, the fitted frequency difference (df) can still be used, but the temperature conversion will not be valid.

Source code in packages/qmri/src/qmri/thermometry/multiecho.py

def fit_multiecho_thermometry(
    signal: NDArray[np.floating],
    echo_times: NDArray[np.floating],
    magnetic_field_tesla: float,
    method: Literal["single", "bootstrap"] = "single",
    n_bootstrap: int = 100,
    r_squared_threshold: float = R_SQUARED_THRESHOLD,
    df_init: DfInitMethod = "multistart",
) -> MultiEchoResult:
    """Fit multi-echo signal to dual-resonance model for thermometry.

    This function fits the dual-resonance signal model to multi-echo
    magnitude data and converts the fitted frequency difference to
    temperature using the ethylene glycol calibration.

    Args:
        signal: Multi-echo magnitude signal array. Shape should be (n_echoes,).
        echo_times: Array of echo times in seconds.
        magnetic_field_tesla: Magnetic field strength in Tesla.
        method: Fitting method. Options:
            - "single": Single least-squares fit (default).
            - "bootstrap": Bootstrap resampling for uncertainty estimation.
        n_bootstrap: Number of bootstrap samples (default 100).
            Only used when method="bootstrap".
        r_squared_threshold: Minimum R² for accepting a fit (default 0.9).
            For bootstrap, samples below threshold are excluded.
        df_init: Frequency starting-value strategy — ``"multistart"`` (default),
            ``"fixed"`` or ``"lombscargle"``. See :data:`DfInitMethod`.

    Returns:
        MultiEchoResult containing temperature, uncertainty, and fit parameters.

    Raises:
        ValueError: If signal and echo_times have different lengths.

    Example:
        ```python
        import numpy as np
        from qmri.thermometry.multiecho import (
            thermometry_signal_model,
            calculate_df_from_temperature,
            fit_multiecho_thermometry,
        )

        # Generate synthetic data at 25°C
        magnetic_field = 3.0
        echo_times = np.linspace(0.001, 0.024, 24)
        df_true = calculate_df_from_temperature(25.0, magnetic_field)
        signal = thermometry_signal_model(
            echo_times, 1.0, 0.5, 50.0, 100.0, df_true, 45.0
        )

        # Fit the model
        result = fit_multiecho_thermometry(
            signal, echo_times, magnetic_field, method="single"
        )
        temp = result.temperature
        uncert = result.temperature_uncertainty
        print(f"Temperature: {temp:.1f} ± {uncert:.2f} °C")
        print(f"R²: {result.r_squared:.4f}")

        # With bootstrap uncertainty
        result_bs = fit_multiecho_thermometry(
            signal, echo_times, magnetic_field,
            method="bootstrap", n_bootstrap=50
        )
        print(f"Bootstrap uncertainty: {result_bs.temperature_uncertainty:.2f} °C")
        ```

    Note:
        The temperature-frequency calibration is specific to ethylene glycol.
        For other substances, the fitted frequency difference (df) can still
        be used, but the temperature conversion will not be valid.
    """
    signal = np.asarray(signal)
    echo_times = np.asarray(echo_times)

    if len(signal) != len(echo_times):
        msg = (
            f"Signal length ({len(signal)}) must match "
            f"echo_times length ({len(echo_times)})"
        )
        raise ValueError(msg)

    if method == "single":
        popt, pcov, r_squared_val = _fit_thermometry(echo_times, signal, df_init)
        df = popt[4]
        temperature = float(calculate_temperature_from_df(df, magnetic_field_tesla))

        # Calculate uncertainty from covariance matrix
        param_uncertainties = np.sqrt(np.diag(pcov))
        df_uncert = param_uncertainties[4]
        df_uncertainty = df_uncert if not np.isnan(df_uncert) else np.nan
        temp_uncertainty = calculate_temperature_uncertainty(
            df_uncertainty, magnetic_field_tesla
        )

        return MultiEchoResult(
            temperature=temperature,
            temperature_uncertainty=temp_uncertainty,
            df=float(df),
            r_squared=r_squared_val,
            fitted_params=popt,
            n_bootstrap=None,
        )

    elif method == "bootstrap":
        rng = np.random.default_rng(seed=RANDOM_SEED)
        n_points = len(signal)

        temperatures_list: list[float] = []
        r_squared_list: list[float] = []
        fitted_params_list: list[NDArray[np.floating]] = []

        for _ in range(n_bootstrap):
            # Resample with replacement
            indices = rng.choice(n_points, size=n_points, replace=True)
            resampled_echo_times = echo_times[indices]
            resampled_signal = signal[indices]

            # Sort by echo time for fitting
            sort_idx = np.argsort(resampled_echo_times)
            resampled_echo_times = resampled_echo_times[sort_idx]
            resampled_signal = resampled_signal[sort_idx]

            # Fit
            popt, _, r_squared_val = _fit_thermometry(
                resampled_echo_times, resampled_signal, df_init
            )

            df = popt[4]
            temp = float(calculate_temperature_from_df(df, magnetic_field_tesla))

            temperatures_list.append(temp)
            r_squared_list.append(r_squared_val)
            fitted_params_list.append(popt)

        temperatures = np.array(temperatures_list)
        r_squared_values = np.array(r_squared_list)
        fitted_params_arr = np.array(fitted_params_list)

        # Filter by R² threshold
        good_fits = r_squared_values >= r_squared_threshold

        if np.any(good_fits):
            mean_temperature = float(np.mean(temperatures[good_fits]))
            temp_uncertainty = float(np.std(temperatures[good_fits]))
            mean_r_squared = float(np.mean(r_squared_values[good_fits]))
            mean_params = np.mean(fitted_params_arr[good_fits], axis=0)
            mean_df = float(mean_params[4])
        else:
            mean_temperature = np.nan
            temp_uncertainty = np.nan
            mean_r_squared = np.nan
            mean_params = np.array([np.nan] * 6)
            mean_df = np.nan

        return MultiEchoResult(
            temperature=mean_temperature,
            temperature_uncertainty=temp_uncertainty,
            df=mean_df,
            r_squared=mean_r_squared,
            fitted_params=mean_params,
            n_bootstrap=n_bootstrap,
        )

fit_multiecho_thermometry_image ¶

fit_multiecho_thermometry_image(
    signal: NDArray[floating],
    segmentation: NDArray[floating],
    echo_times: NDArray[floating],
    magnetic_field_tesla: float,
    method: RegionAnalysisMethod = "regionwise",
    n_bootstrap: int = 100,
    r_squared_threshold: float = R_SQUARED_THRESHOLD,
    df_init: DfInitMethod = "multistart",
) -> tuple[
    NDArray[floating], list[RegionThermometryResult]
]

Fit multi-echo thermometry over a segmented image volume.

The segmentation defines discrete regions by integer label; label 0 is treated as background and ignored. Each non-zero region is fitted with the dual-resonance model and the fitted frequency difference is converted to temperature with the ethylene-glycol calibration.

The arrays must be co-located in world space: signal is the 4D magnitude volume (nx, ny, nz, n_echoes) and segmentation is the 3D label map (nx, ny, nz). echo_times must have length n_echoes and be in seconds.

Parameters:

Name	Type	Description	Default
`signal`	`NDArray[floating]`	Multi-echo magnitude data, shape `(nx, ny, nz, n_echoes)`.	required
`segmentation`	`NDArray[floating]`	Integer label map, shape `(nx, ny, nz)`.	required
`echo_times`	`NDArray[floating]`	Echo times in seconds, shape `(n_echoes,)`.	required
`magnetic_field_tesla`	`float`	Magnetic field strength \(B_0\) in Tesla.	required
`method`	`RegionAnalysisMethod`	Analysis method: `"regionwise"`: fit the mean signal within each region once and assign the resulting temperature to every voxel in the region. Uncertainty comes from the fitted Δf covariance. `"voxelwise"`: fit each voxel independently; the region summary is an inverse-variance weighted mean of voxel temperatures. `"regionwise_bootstrap"`: resample region voxels with replacement, fit each sample's mean signal, and summarise with the mean and standard deviation over samples with \(R^2 \geq\) `r_squared_threshold`.	`'regionwise'`
`n_bootstrap`	`int`	Number of bootstrap samples (`regionwise_bootstrap` only).	`100`
`r_squared_threshold`	`float`	Minimum R² for a fit to contribute to `mean_fitted_params` and to bootstrap summaries.	`R_SQUARED_THRESHOLD`
`df_init`	`DfInitMethod`	Frequency starting-value strategy — `"multistart"` (default), `"fixed"` or `"lombscargle"`. See :data:`DfInitMethod`.	`'multistart'`

Returns:

Type	Description
`NDArray[floating]`	A tuple `(temperature_map, results)` where `temperature_map` is a 3D
`list[RegionThermometryResult]`	array of temperatures in °C co-located with the segmentation, and
`tuple[NDArray[floating], list[RegionThermometryResult]]`	`results` is a list of :class:`RegionThermometryResult`, one per
`tuple[NDArray[floating], list[RegionThermometryResult]]`	non-empty region (in ascending label order).

Raises:

Type	Description
`ValueError`	If the array dimensions or echo-time length are inconsistent, or if `method` is not recognised.

Source code in packages/qmri/src/qmri/thermometry/multiecho.py

def fit_multiecho_thermometry_image(
    signal: NDArray[np.floating],
    segmentation: NDArray[np.floating],
    echo_times: NDArray[np.floating],
    magnetic_field_tesla: float,
    method: RegionAnalysisMethod = "regionwise",
    n_bootstrap: int = 100,
    r_squared_threshold: float = R_SQUARED_THRESHOLD,
    df_init: DfInitMethod = "multistart",
) -> tuple[NDArray[np.floating], list[RegionThermometryResult]]:
    r"""Fit multi-echo thermometry over a segmented image volume.

    The segmentation defines discrete regions by integer label; label ``0`` is
    treated as background and ignored. Each non-zero region is fitted with the
    dual-resonance model and the fitted frequency difference is converted to
    temperature with the ethylene-glycol calibration.

    The arrays must be co-located in world space: ``signal`` is the 4D
    magnitude volume ``(nx, ny, nz, n_echoes)`` and ``segmentation`` is the 3D
    label map ``(nx, ny, nz)``. ``echo_times`` must have length ``n_echoes`` and
    be in seconds.

    Args:
        signal: Multi-echo magnitude data, shape ``(nx, ny, nz, n_echoes)``.
        segmentation: Integer label map, shape ``(nx, ny, nz)``.
        echo_times: Echo times in seconds, shape ``(n_echoes,)``.
        magnetic_field_tesla: Magnetic field strength $B_0$ in Tesla.
        method: Analysis method:

            - ``"regionwise"``: fit the mean signal within each region once and
              assign the resulting temperature to every voxel in the region.
              Uncertainty comes from the fitted Δf covariance.
            - ``"voxelwise"``: fit each voxel independently; the region summary
              is an inverse-variance weighted mean of voxel temperatures.
            - ``"regionwise_bootstrap"``: resample region voxels with
              replacement, fit each sample's mean signal, and summarise with the
              mean and standard deviation over samples with $R^2 \geq$
              ``r_squared_threshold``.
        n_bootstrap: Number of bootstrap samples (``regionwise_bootstrap`` only).
        r_squared_threshold: Minimum R² for a fit to contribute to
            ``mean_fitted_params`` and to bootstrap summaries.
        df_init: Frequency starting-value strategy — ``"multistart"`` (default),
            ``"fixed"`` or ``"lombscargle"``. See :data:`DfInitMethod`.

    Returns:
        A tuple ``(temperature_map, results)`` where ``temperature_map`` is a 3D
        array of temperatures in °C co-located with the segmentation, and
        ``results`` is a list of :class:`RegionThermometryResult`, one per
        non-empty region (in ascending label order).

    Raises:
        ValueError: If the array dimensions or echo-time length are inconsistent,
            or if ``method`` is not recognised.
    """
    signal = np.asarray(signal, dtype=np.float64)
    segmentation = np.asarray(segmentation)
    echo_times = np.asarray(echo_times, dtype=np.float64)

    if signal.ndim != 4:
        msg = f"signal must be 4D (nx, ny, nz, n_echoes), got {signal.ndim}D"
        raise ValueError(msg)
    if segmentation.ndim != 3:
        msg = f"segmentation must be 3D (nx, ny, nz), got {segmentation.ndim}D"
        raise ValueError(msg)
    if segmentation.shape != signal.shape[:3]:
        msg = (
            "segmentation shape must match the spatial shape of signal: "
            f"{segmentation.shape} != {signal.shape[:3]}"
        )
        raise ValueError(msg)
    if echo_times.shape[0] != signal.shape[-1]:
        msg = (
            f"echo_times length ({echo_times.shape[0]}) must match the number "
            f"of echoes ({signal.shape[-1]})"
        )
        raise ValueError(msg)
    if method not in ("regionwise", "voxelwise", "regionwise_bootstrap"):
        msg = (
            f"Unknown method: {method!r}. Use 'regionwise', 'voxelwise', or "
            "'regionwise_bootstrap'."
        )
        raise ValueError(msg)

    n_echoes = signal.shape[-1]
    n_params = 6

    temperature_map: NDArray[np.float64] = np.zeros(
        segmentation.shape, dtype=np.float64
    )

    regions = np.unique(segmentation)
    regions = regions[regions != 0]  # exclude background
    results: list[RegionThermometryResult] = []

    for region in regions:
        region_mask = segmentation == region
        region_size = int(np.sum(region_mask))
        if region_size == 0:
            continue

        temperature_values_list: list[float] = []
        uncertainty_values_list: list[float] = []
        fitted_params_list: list[NDArray[np.floating]] = []
        signal_values_list: list[NDArray[np.floating]] = []
        r_squared_list: list[float] = []

        region_temperature: float = float("nan")
        region_uncertainty: float = float("nan")

        if method == "regionwise":
            region_signal = np.array(
                [
                    float(np.mean(signal[..., echo][region_mask]))
                    for echo in range(n_echoes)
                ]
            )
            popt, pcov, r_squared_value = _fit_thermometry(
                echo_times, region_signal, df_init
            )
            df = float(popt[4])
            df_uncertainty = float(np.sqrt(np.diag(pcov))[4])
            region_temperature = float(
                calculate_temperature_from_df(df, magnetic_field_tesla)
            )
            region_uncertainty = calculate_temperature_uncertainty(
                df_uncertainty, magnetic_field_tesla
            )
            temperature_map[region_mask] = region_temperature

            temperature_values_list.append(region_temperature)
            uncertainty_values_list.append(region_uncertainty)
            fitted_params_list.append(popt)
            signal_values_list.append(region_signal)
            r_squared_list.append(r_squared_value)

        elif method == "voxelwise":
            for index in np.argwhere(region_mask):
                i, j, k = (int(index[0]), int(index[1]), int(index[2]))
                voxel_signal = signal[i, j, k, :]
                popt, pcov, r_squared_value = _fit_thermometry(
                    echo_times, voxel_signal, df_init
                )
                df = float(popt[4])
                df_uncertainty = float(np.sqrt(np.diag(pcov))[4])
                voxel_temperature = float(
                    calculate_temperature_from_df(df, magnetic_field_tesla)
                )
                voxel_uncertainty = calculate_temperature_uncertainty(
                    df_uncertainty, magnetic_field_tesla
                )
                temperature_map[i, j, k] = voxel_temperature

                temperature_values_list.append(voxel_temperature)
                uncertainty_values_list.append(voxel_uncertainty)
                fitted_params_list.append(popt)
                signal_values_list.append(voxel_signal)
                r_squared_list.append(r_squared_value)

            uncertainty_values = np.array(uncertainty_values_list)
            weights = np.divide(
                1.0,
                uncertainty_values**2,
                where=uncertainty_values != 0,
                out=np.zeros_like(uncertainty_values),
            )
            region_temperature = float(
                np.average(np.array(temperature_values_list), weights=weights)
            )
            weight_sum = float(np.sum(weights))
            region_uncertainty = (
                float(np.sqrt(1.0 / weight_sum)) if weight_sum > 0 else float("nan")
            )

        else:  # regionwise_bootstrap
            rng = np.random.default_rng(seed=RANDOM_SEED)
            flat_signal = signal.reshape(-1, n_echoes)
            region_flat_indices = np.where(region_mask.flatten())[0]
            for _ in range(n_bootstrap):
                sampled_indices = rng.choice(
                    region_flat_indices, size=region_size, replace=True
                )
                region_signal = np.mean(flat_signal[sampled_indices], axis=0)
                popt, pcov, r_squared_value = _fit_thermometry(
                    echo_times, region_signal, df_init
                )
                df = float(popt[4])
                df_uncertainty = float(np.sqrt(np.diag(pcov))[4])
                temperature_values_list.append(
                    float(calculate_temperature_from_df(df, magnetic_field_tesla))
                )
                uncertainty_values_list.append(
                    calculate_temperature_uncertainty(
                        df_uncertainty, magnetic_field_tesla
                    )
                )
                fitted_params_list.append(popt)
                signal_values_list.append(region_signal)
                r_squared_list.append(r_squared_value)

            temperature_values = np.array(temperature_values_list)
            r_squared_arr = np.array(r_squared_list)
            passing = r_squared_arr >= r_squared_threshold
            if np.any(passing):
                region_temperature = float(np.mean(temperature_values[passing]))
                region_uncertainty = float(np.std(temperature_values[passing]))
            temperature_map[region_mask] = region_temperature

        r_squared = np.array(r_squared_list)
        fitted_params = np.array(fitted_params_list)
        passing_mask = r_squared >= r_squared_threshold
        if fitted_params.size > 0 and np.any(passing_mask):
            mean_fitted_params = np.mean(fitted_params[passing_mask], axis=0)
        else:
            mean_fitted_params = np.full(n_params, np.nan)

        results.append(
            RegionThermometryResult(
                region_id=int(region),
                region_size=region_size,
                temperature=region_temperature,
                temperature_uncertainty=region_uncertainty,
                coverage_factor=1.0,
                temperature_values=np.array(temperature_values_list),
                temperature_uncertainty_values=np.array(uncertainty_values_list),
                r_squared=r_squared,
                fitted_params=fitted_params,
                mean_fitted_params=mean_fitted_params,
                signal_values=np.array(signal_values_list),
            )
        )

    return temperature_map, results

lsq_fit_thermometry_signal_model ¶

lsq_fit_thermometry_signal_model(
    echo_times: NDArray[floating],
    signal_values: NDArray[floating],
    initial_guess: list[float],
) -> tuple[NDArray[floating], NDArray[floating], float]

Perform least squares fit of the dual-resonance signal model.

Fits the signal data to the thermometry signal model using scipy.optimize.curve_fit with bounded parameters.

Parameters:

Name	Type	Description	Default
`echo_times`	`NDArray[floating]`	Array of echo times in seconds.	required
`signal_values`	`NDArray[floating]`	Array of signal values at each echo time.	required
`initial_guess`	`list[float]`	Initial parameter estimates [amplitude_1, amplitude_2, r2star_1, r2star_2, df, dphi_deg].	required

Returns:

Type	Description
`tuple[NDArray[floating], NDArray[floating], float]`	Tuple containing: - popt: Optimal parameters [A1, A2, R21, R22, df, dphi]. - pcov: Covariance matrix of the parameters. - r_squared: Coefficient of determination (R²) of the fit.

If the fit fails to converge, returns arrays of NaN values.

Example

import numpy as np
from qmri.thermometry.multiecho import (
    thermometry_signal_model,
    lsq_fit_thermometry_signal_model,
)

# Generate synthetic data
echo_times = np.linspace(0.001, 0.024, 24)
true_params = [1.0, 0.5, 50.0, 100.0, 200.0, 45.0]
signal = thermometry_signal_model(echo_times, *true_params)

# Fit the model
initial_guess = [0.8, 0.4, 40.0, 80.0, 180.0, 30.0]
popt, pcov, r2 = lsq_fit_thermometry_signal_model(
    echo_times, signal, initial_guess
)
print(f"R²: {r2:.4f}")

Source code in packages/qmri/src/qmri/thermometry/multiecho.py

def lsq_fit_thermometry_signal_model(
    echo_times: NDArray[np.floating],
    signal_values: NDArray[np.floating],
    initial_guess: list[float],
) -> tuple[NDArray[np.floating], NDArray[np.floating], float]:
    """Perform least squares fit of the dual-resonance signal model.

    Fits the signal data to the thermometry signal model using
    scipy.optimize.curve_fit with bounded parameters.

    Args:
        echo_times: Array of echo times in seconds.
        signal_values: Array of signal values at each echo time.
        initial_guess: Initial parameter estimates
            [amplitude_1, amplitude_2, r2star_1, r2star_2, df, dphi_deg].

    Returns:
        Tuple containing:
            - popt: Optimal parameters [A1, A2, R2*1, R2*2, df, dphi].
            - pcov: Covariance matrix of the parameters.
            - r_squared: Coefficient of determination (R²) of the fit.

    If the fit fails to converge, returns arrays of NaN values.

    Example:
        ```python
        import numpy as np
        from qmri.thermometry.multiecho import (
            thermometry_signal_model,
            lsq_fit_thermometry_signal_model,
        )

        # Generate synthetic data
        echo_times = np.linspace(0.001, 0.024, 24)
        true_params = [1.0, 0.5, 50.0, 100.0, 200.0, 45.0]
        signal = thermometry_signal_model(echo_times, *true_params)

        # Fit the model
        initial_guess = [0.8, 0.4, 40.0, 80.0, 180.0, 30.0]
        popt, pcov, r2 = lsq_fit_thermometry_signal_model(
            echo_times, signal, initial_guess
        )
        print(f"R²: {r2:.4f}")
        ```
    """
    max_amplitude = np.max(signal_values)
    bounds = (
        [0, 0, 1e-3, 1e-3, 0, -360],
        [10 * max_amplitude, 10 * max_amplitude, 1000, 1000, _DF_UPPER_BOUND, 360],
    )
    try:
        popt, pcov, *_ = curve_fit(
            thermometry_signal_model,
            echo_times,
            signal_values,
            p0=initial_guess,
            bounds=bounds,
            maxfev=10000,
        )
    except RuntimeError:
        popt = np.array([np.nan] * len(initial_guess), dtype=np.float64)
        pcov = np.full((len(initial_guess), len(initial_guess)), np.nan)

    r_squared_val = calculate_r_squared(
        signal_values, thermometry_signal_model(echo_times, *popt)
    )
    return popt, pcov, float(r_squared_val)

thermometry_signal_model ¶

thermometry_signal_model(
    t: NDArray[floating],
    amplitude_1: float,
    amplitude_2: float,
    r2star_1: float,
    r2star_2: float,
    df: float | NDArray[floating],
    dphi_deg: float,
) -> NDArray[floating]

Calculate the dual-resonance signal at time t.

Implements the signal model equation:

\[S(t) = \sqrt{A_1^2 e^{-2R_{2,1}^* t} + A_2^2 e^{-2R_{2,2}^* t} + 2 A_1 A_2 e^{-(R_{2,1}^* + R_{2,2}^*) t} \cos(2\pi \Delta f t + \Delta\phi)}\]

Parameters:

Name	Type	Description	Default
`t`	`NDArray[floating]`	Echo times in seconds.	required
`amplitude_1`	`float`	Amplitude of the first signal component.	required
`amplitude_2`	`float`	Amplitude of the second signal component.	required
`r2star_1`	`float`	R2* relaxation rate of the first component (1/s).	required
`r2star_2`	`float`	R2* relaxation rate of the second component (1/s).	required
`df`	`float \| NDArray[floating]`	Frequency difference between the two components (Hz).	required
`dphi_deg`	`float`	Phase difference between the components (degrees).	required

Returns:

Type	Description
`NDArray[floating]`	Signal magnitude at each echo time.

Example

import numpy as np
from qmri.thermometry.multiecho import thermometry_signal_model

echo_times = np.linspace(0.001, 0.024, 24)
signal = thermometry_signal_model(
    t=echo_times,
    amplitude_1=1.0,
    amplitude_2=0.5,
    r2star_1=50.0,
    r2star_2=100.0,
    df=200.0,
    dphi_deg=45.0,
)

Source code in packages/qmri/src/qmri/thermometry/multiecho.py

def thermometry_signal_model(
    t: NDArray[np.floating],
    amplitude_1: float,
    amplitude_2: float,
    r2star_1: float,
    r2star_2: float,
    df: float | NDArray[np.floating],
    dphi_deg: float,
) -> NDArray[np.floating]:
    r"""Calculate the dual-resonance signal at time t.

    Implements the signal model equation:

    $$S(t) = \sqrt{A_1^2 e^{-2R_{2,1}^* t} + A_2^2 e^{-2R_{2,2}^* t}
           + 2 A_1 A_2 e^{-(R_{2,1}^* + R_{2,2}^*) t}
           \cos(2\pi \Delta f t + \Delta\phi)}$$

    Args:
        t: Echo times in seconds.
        amplitude_1: Amplitude of the first signal component.
        amplitude_2: Amplitude of the second signal component.
        r2star_1: R2* relaxation rate of the first component (1/s).
        r2star_2: R2* relaxation rate of the second component (1/s).
        df: Frequency difference between the two components (Hz).
        dphi_deg: Phase difference between the components (degrees).

    Returns:
        Signal magnitude at each echo time.

    Example:
        ```python
        import numpy as np
        from qmri.thermometry.multiecho import thermometry_signal_model

        echo_times = np.linspace(0.001, 0.024, 24)
        signal = thermometry_signal_model(
            t=echo_times,
            amplitude_1=1.0,
            amplitude_2=0.5,
            r2star_1=50.0,
            r2star_2=100.0,
            df=200.0,
            dphi_deg=45.0,
        )
        ```
    """
    dphi_rad = np.deg2rad(dphi_deg)
    radicand: NDArray[np.float64] = (
        amplitude_1**2 * np.exp(-2 * r2star_1 * t)
        + amplitude_2**2 * np.exp(-2 * r2star_2 * t)
        + 2
        * amplitude_1
        * amplitude_2
        * np.exp(-(r2star_1 + r2star_2) * t)
        * np.cos(2 * np.pi * df * t + dphi_rad)
    )
    # Prevent negative values under the square root
    radicand = np.maximum(radicand, 0)
    return np.sqrt(radicand)

calculate_temperature ¶

calculate_temperature(
    phase_difference: NDArray[floating] | float,
    echo_time: NDArray[floating] | float,
    magnetic_field: float,
    baseline_temperature: NDArray[floating] | float = 37.0,
    prf_coefficient: float = PRF_THERMAL_COEFFICIENT,
) -> PRFResult

Calculate temperature change from phase difference.

Implements the inverse PRF equation:

\[\Delta T = \frac{\Delta\phi}{\gamma \alpha B_0 \cdot TE}\]

where the phase difference is \(\phi_{heated} - \phi_{baseline}\).

Parameters:

Name	Type	Description	Default
`phase_difference`	`NDArray[floating] \| float`	Phase difference between heated and baseline images in radians. Calculated as phase_heated - phase_baseline.	required
`echo_time`	`NDArray[floating] \| float`	Echo time (TE) in seconds.	required
`magnetic_field`	`float`	Main magnetic field strength (B0) in Tesla.	required
`baseline_temperature`	`NDArray[floating] \| float`	Baseline temperature in degrees Celsius (default 37.0). This is the reference temperature at which the baseline phase image was acquired. The absolute temperature is calculated as baseline_temperature + temperature_change.	`37.0`
`prf_coefficient`	`float`	PRF thermal coefficient in per degree Celsius (default -0.01e-6). The standard value for water is approximately -0.01 ppm/°C.	`PRF_THERMAL_COEFFICIENT`

Returns:

Type	Description
`PRFResult`	Result containing temperature change and phase difference.

Example

import numpy as np
from qmri.thermometry import prf

# Calculate temperature from measured phase difference
result = prf.calculate_temperature(
    phase_difference=-0.16,  # radians (negative = heating)
    echo_time=0.020,  # 20 ms
    magnetic_field=3.0,
)
print(f"Temperature change: {result.temperature_change:.1f} °C")

Important considerations for PRF thermometry:

Phase wrapping: Phase values are typically wrapped to [-pi, pi]. For large temperature changes, phase unwrapping may be required.
Reference tissue: PRF thermometry provides relative temperature changes. A reference tissue (e.g., subcutaneous fat) with known temperature can be used for drift correction.
Fat signal: Fat does not exhibit PRF shift with temperature. Fat suppression or water-fat separation is often used.
Motion: Motion between baseline and heated acquisitions causes artefacts. Multi-baseline or referenceless methods can help.
Field drift: B0 drift over time causes apparent temperature changes. Drift correction using reference regions is recommended.

signal_phase_shift ¶

signal_phase_shift(
    temperature_change: NDArray[floating] | float,
    echo_time: NDArray[floating] | float,
    magnetic_field: float,
    prf_coefficient: float = PRF_THERMAL_COEFFICIENT,
) -> NDArray[floating]

Calculate temperature-induced phase shift.

Implements the PRF phase shift equation:

\[\Delta\phi = -\gamma \alpha B_0 \Delta T \cdot TE\]

Note that the negative sign arises from the negative PRF thermal coefficient, which means increasing temperature causes a negative phase shift (frequency decrease).

Parameters:

Name	Type	Description	Default
`temperature_change`	`NDArray[floating] \| float`	Temperature change in degrees Celsius.	required
`echo_time`	`NDArray[floating] \| float`	Echo time (TE) in seconds.	required
`magnetic_field`	`float`	Main magnetic field strength (B0) in Tesla.	required
`prf_coefficient`	`float`	PRF thermal coefficient in per degree Celsius (default -0.01e-6). The standard value for water is approximately -0.01 ppm/°C.	`PRF_THERMAL_COEFFICIENT`

Returns:

Type	Description
`NDArray[floating]`	Phase shift in radians.

Example

import numpy as np
from qmri.thermometry import prf

# Calculate phase shift for 10°C temperature increase at 3T
phase_shift = prf.signal_phase_shift(
    temperature_change=10.0,
    echo_time=0.020,  # 20 ms
    magnetic_field=3.0,
)
print(f"Phase shift: {np.degrees(phase_shift):.2f} degrees")

The phase shift is proportional to:

Temperature change (linear)
Echo time (longer TE gives larger phase shift but lower SNR)
Field strength (higher field gives larger phase shift)

Typical phase shifts at 3T with TE=20ms are approximately 5-10 degrees per degree Celsius temperature change.

Source code in packages/qmri/src/qmri/thermometry/prf.py

def signal_phase_shift(
    temperature_change: NDArray[np.floating] | float,
    echo_time: NDArray[np.floating] | float,
    magnetic_field: float,
    prf_coefficient: float = PRF_THERMAL_COEFFICIENT,
) -> NDArray[np.floating]:
    r"""Calculate temperature-induced phase shift.

    Implements the PRF phase shift equation:

    $$\Delta\phi = -\gamma \alpha B_0 \Delta T \cdot TE$$

    Note that the negative sign arises from the negative PRF thermal
    coefficient, which means increasing temperature causes a negative
    phase shift (frequency decrease).

    Args:
        temperature_change: Temperature change in degrees Celsius.
        echo_time: Echo time (TE) in seconds.
        magnetic_field: Main magnetic field strength (B0) in Tesla.
        prf_coefficient: PRF thermal coefficient in per degree Celsius
            (default -0.01e-6). The standard value for water is
            approximately -0.01 ppm/°C.

    Returns:
        Phase shift in radians.

    Example:
        ```python
        import numpy as np
        from qmri.thermometry import prf

        # Calculate phase shift for 10°C temperature increase at 3T
        phase_shift = prf.signal_phase_shift(
            temperature_change=10.0,
            echo_time=0.020,  # 20 ms
            magnetic_field=3.0,
        )
        print(f"Phase shift: {np.degrees(phase_shift):.2f} degrees")
        ```

    The phase shift is proportional to:

    - Temperature change (linear)
    - Echo time (longer TE gives larger phase shift but lower SNR)
    - Field strength (higher field gives larger phase shift)

    Typical phase shifts at 3T with TE=20ms are approximately 5-10 degrees
    per degree Celsius temperature change.
    """
    delta_t = np.asarray(temperature_change)
    te = np.asarray(echo_time)
    b0 = magnetic_field
    alpha = prf_coefficient

    # Convert gyromagnetic ratio to rad/s/T for phase calculation
    # gamma in Hz/T * 2*pi gives rad/s/T
    gamma_rad = 2 * np.pi * GAMMA_H

    # Phase shift: delta_phi = gamma * alpha * B0 * delta_T * TE
    # Note: alpha is already negative for water, so increasing T gives
    # negative phase shift
    result: NDArray[np.floating] = gamma_rad * alpha * b0 * delta_t * te
    return result

qmri.thermometry¶

thermometry ¶

R_SQUARED_THRESHOLD module-attribute ¶

RANDOM_SEED module-attribute ¶

DfInitMethod module-attribute ¶

RegionAnalysisMethod module-attribute ¶

MultiEchoResult dataclass ¶

RegionThermometryResult dataclass ¶

to_dict ¶

PRFResult dataclass ¶

calculate_df_from_temperature ¶

calculate_temperature_from_df ¶

calculate_temperature_uncertainty ¶

fit_multiecho_thermometry ¶

fit_multiecho_thermometry_image ¶

lsq_fit_thermometry_signal_model ¶

thermometry_signal_model ¶

calculate_temperature ¶

signal_phase_shift ¶

R_SQUARED_THRESHOLD `module-attribute` ¶

RANDOM_SEED `module-attribute` ¶

DfInitMethod `module-attribute` ¶

RegionAnalysisMethod `module-attribute` ¶

MultiEchoResult `dataclass` ¶

RegionThermometryResult `dataclass` ¶

PRFResult `dataclass` ¶