Metrics API¶

dosemetrics.metrics exposes domain modules for reference-free computations and direct compare_* functions for the named reference-based plan metrics.

Public module	Responsibility	Dose-reference convention
`conformity`	Single-plan coverage and conformity quantities	Reference-free
`dose_comparison`	General image and voxel comparisons plus image descriptors	`compare_` is reference-based; `compute_` is reference-free
`dvh`	DVH construction, queries, summaries, and comparisons	`compute_` is reference-free; `compare_` is reference-based
`gamma`	Gamma maps, passing-rate summaries, and statistics	Gamma-map creation is reference-based
`geometric`	Structure overlap, volume, and surface metrics	Structure geometry; no dose reference
`homogeneity`	Single-plan target homogeneity and dose falloff	Reference-free

The generated sections below are the authoritative public signatures and defaults from the implementation. For clinical interpretation and reference semantics, see the Metric Framework.

Direct plan-comparison functions¶

Import these functions without an additional module prefix:

from dosemetrics.metrics import compare_ptv_dose

distance_gy = compare_ptv_dose(reference, evaluated, ptv)

compare_ptv_dose ¶

compare_ptv_dose(reference: Dose, evaluated: Dose, ptv: Structure) -> float

Compute the absolute difference in mean PTV dose.

.. math::

|\bar D_{\mathrm{evaluated}} - \bar D_{\mathrm{reference}}|

Source code in src/dosemetrics/metrics/comparison.py

def compare_ptv_dose(
    reference: Dose,
    evaluated: Dose,
    ptv: Structure,
) -> float:
    """Compute the absolute difference in mean PTV dose.

    .. math::

        |\\bar D_{\\mathrm{evaluated}} - \\bar D_{\\mathrm{reference}}|
    """

    _validate_dose_geometry(reference, evaluated)
    reference_values = reference.get_dose_in_structure(ptv)
    evaluated_values = evaluated.get_dose_in_structure(ptv)
    if len(reference_values) == 0:
        return float("nan")
    return float(abs(np.mean(evaluated_values) - np.mean(reference_values)))

compare_paddick_conformity_index ¶

compare_paddick_conformity_index(reference: Dose, evaluated: Dose, ptv: Structure, prescription_dose: float) -> float

Compute absolute distance between two Paddick conformity indices.

The underlying index is

.. math::

CI = \frac{V_{PTV,PIV}}{V_{PTV}}
     \frac{V_{PTV,PIV}}{V_{PIV}}.

Source code in src/dosemetrics/metrics/comparison.py

def compare_paddick_conformity_index(
    reference: Dose,
    evaluated: Dose,
    ptv: Structure,
    prescription_dose: float,
) -> float:
    """Compute absolute distance between two Paddick conformity indices.

    The underlying index is

    .. math::

        CI = \\frac{V_{PTV,PIV}}{V_{PTV}}
             \\frac{V_{PTV,PIV}}{V_{PIV}}.
    """

    _validate_dose_geometry(reference, evaluated)
    reference_ci = compute_paddick_conformity_index(reference, ptv, prescription_dose)
    evaluated_ci = compute_paddick_conformity_index(evaluated, ptv, prescription_dose)
    return float(abs(evaluated_ci - reference_ci))

compare_paddick_gradient_index ¶

compare_paddick_gradient_index(reference: Dose, evaluated: Dose, prescription_dose: float) -> float

Compute absolute distance between Paddick gradient indices.

The target argument accepted by :func:compute_gradient_index is not used by its volume-ratio definition, so this implementation calculates the specified V_50% / V_100% ratio directly.

Source code in src/dosemetrics/metrics/comparison.py

def compare_paddick_gradient_index(
    reference: Dose,
    evaluated: Dose,
    prescription_dose: float,
) -> float:
    """Compute absolute distance between Paddick gradient indices.

    The target argument accepted by :func:`compute_gradient_index` is not used
    by its volume-ratio definition, so this implementation calculates the
    specified ``V_50% / V_100%`` ratio directly.
    """

    _validate_dose_geometry(reference, evaluated)

    def gradient_index(dose: Dose) -> float:
        v_100 = np.count_nonzero(dose.dose_array >= prescription_dose)
        v_50 = np.count_nonzero(dose.dose_array >= 0.5 * prescription_dose)
        if v_100 == 0:
            return float("inf")
        return float(v_50 / v_100)

    reference_gi = gradient_index(reference)
    evaluated_gi = gradient_index(evaluated)
    return float(abs(evaluated_gi - reference_gi))

compare_homogeneity_index ¶

compare_homogeneity_index(reference: Dose, evaluated: Dose, ptv: Structure) -> float

Compute distance between HI = (D2 - D98) / D50 values.

Source code in src/dosemetrics/metrics/comparison.py

def compare_homogeneity_index(
    reference: Dose,
    evaluated: Dose,
    ptv: Structure,
) -> float:
    """Compute distance between ``HI = (D2 - D98) / D50`` values."""

    _validate_dose_geometry(reference, evaluated)
    reference_hi = compute_homogeneity_index(reference, ptv)
    evaluated_hi = compute_homogeneity_index(evaluated, ptv)
    return float(abs(evaluated_hi - reference_hi))

compare_body_rmse ¶

compare_body_rmse(reference: Dose, evaluated: Dose, body: Structure) -> float

Compute root mean squared error over body-mask voxels, in Gy.

Source code in src/dosemetrics/metrics/comparison.py

def compare_body_rmse(
    reference: Dose,
    evaluated: Dose,
    body: Structure,
) -> float:
    """Compute root mean squared error over body-mask voxels, in Gy."""

    _validate_dose_geometry(reference, evaluated)
    reference_values = reference.get_dose_in_structure(body)
    evaluated_values = evaluated.get_dose_in_structure(body)
    if len(reference_values) == 0:
        return float("nan")
    return float(np.sqrt(np.mean((evaluated_values - reference_values) ** 2)))

compare_gamma ¶

compare_gamma(reference: Dose, evaluated: Dose, body: Optional[Structure] = None, dose_criterion_percent: float = 3.0, distance_criterion_mm: float = 3.0, dose_threshold_percent: float = 0.0, global_normalization: bool = True, max_search_distance_mm: Optional[float] = None) -> float

Compute the percentage of reference voxels with gamma <= 1.

Defaults implement 3%/3 mm global gamma and evaluate every reference voxel. Pass body to restrict the reported passing rate to the body mask. A low-dose threshold can be requested explicitly, but is disabled by default because it is not part of this metric definition.

Source code in src/dosemetrics/metrics/comparison.py

def compare_gamma(
    reference: Dose,
    evaluated: Dose,
    body: Optional[Structure] = None,
    dose_criterion_percent: float = 3.0,
    distance_criterion_mm: float = 3.0,
    dose_threshold_percent: float = 0.0,
    global_normalization: bool = True,
    max_search_distance_mm: Optional[float] = None,
) -> float:
    """Compute the percentage of reference voxels with ``gamma <= 1``.

    Defaults implement 3%/3 mm global gamma and evaluate every reference
    voxel. Pass ``body`` to restrict the reported passing rate to the body
    mask. A low-dose threshold can be requested explicitly, but is disabled by
    default because it is not part of this metric definition.
    """

    _validate_dose_geometry(reference, evaluated)
    gamma = compare_gamma_index(
        reference,
        evaluated,
        dose_criterion_percent=dose_criterion_percent,
        distance_criterion_mm=distance_criterion_mm,
        dose_threshold_percent=dose_threshold_percent,
        global_normalization=global_normalization,
        max_search_distance_mm=max_search_distance_mm,
    )
    if body is not None:
        # This also validates compatibility between the body and dose grids.
        reference.get_dose_in_structure(body)
        gamma = gamma[body.mask]
    return compute_gamma_passing_rate(gamma, threshold=1.0)

compare_oar_constraints ¶

compare_oar_constraints(reference_satisfaction: Union[Sequence[bool], Mapping[str, bool]], evaluated_satisfaction: Union[Sequence[bool], Mapping[str, bool]], expected_count: Optional[int] = 38) -> float

Compare OAR-constraint states and return their disagreement fraction.

Mappings must contain identical keys. Sequences are compared by position. By default, the inputs must each contain the 38 CORSAIR-derived states in the head-and-neck protocol. Set expected_count=None for a deliberately different constraint protocol.

Source code in src/dosemetrics/metrics/comparison.py

def compare_oar_constraints(
    reference_satisfaction: Union[Sequence[bool], Mapping[str, bool]],
    evaluated_satisfaction: Union[Sequence[bool], Mapping[str, bool]],
    expected_count: Optional[int] = 38,
) -> float:
    """Compare OAR-constraint states and return their disagreement fraction.

    Mappings must contain identical keys. Sequences are compared by position.
    By default, the inputs must each contain the 38 CORSAIR-derived states in
    the head-and-neck protocol. Set ``expected_count=None`` for a deliberately
    different constraint protocol.
    """

    if isinstance(reference_satisfaction, Mapping):
        if not isinstance(evaluated_satisfaction, Mapping):
            raise TypeError("Both satisfaction inputs must use the same container type")
        if set(reference_satisfaction) != set(evaluated_satisfaction):
            raise ValueError("Constraint mappings must contain identical keys")
        keys = tuple(reference_satisfaction)
        reference = np.asarray(
            [reference_satisfaction[key] for key in keys], dtype=bool
        )
        evaluated = np.asarray(
            [evaluated_satisfaction[key] for key in keys], dtype=bool
        )
    else:
        if isinstance(evaluated_satisfaction, Mapping):
            raise TypeError("Both satisfaction inputs must use the same container type")
        reference = np.asarray(reference_satisfaction, dtype=bool)
        evaluated = np.asarray(evaluated_satisfaction, dtype=bool)

    if reference.ndim != 1 or evaluated.ndim != 1:
        raise ValueError("Constraint satisfaction inputs must be one-dimensional")
    if len(reference) == 0:
        raise ValueError("At least one constraint is required")
    if len(reference) != len(evaluated):
        raise ValueError(
            "Constraint satisfaction inputs must have the same number of entries"
        )
    if expected_count is not None and len(reference) != expected_count:
        raise ValueError(
            f"Expected {expected_count} constraint states, got {len(reference)}"
        )
    return float(np.mean(reference != evaluated))

compare_oar_dvh_auc ¶

compare_oar_dvh_auc(reference: Dose, evaluated: Dose, oar: Structure, num_bins: int = 100) -> float

Return |AUC_evaluated - AUC_reference| for one OAR, in Gy.

Source code in src/dosemetrics/metrics/dvh.py

def compare_oar_dvh_auc(
    reference: Dose,
    evaluated: Dose,
    oar: Structure,
    num_bins: int = 100,
) -> float:
    """Return ``|AUC_evaluated - AUC_reference|`` for one OAR, in Gy."""
    _validate_comparable_doses(reference, evaluated)
    if num_bins < 2:
        raise ValueError("num_bins must be at least 2")
    reference_values = reference.get_dose_in_structure(oar)
    evaluated_values = evaluated.get_dose_in_structure(oar)
    if len(reference_values) == 0:
        return float("nan")
    min_dose = float(min(np.min(reference_values), np.min(evaluated_values)))
    max_dose = float(max(np.max(reference_values), np.max(evaluated_values)))
    if np.isclose(min_dose, max_dose):
        return 0.0
    bins = np.linspace(min_dose, max_dose, num_bins)
    reference_dvh = np.asarray([np.mean(reference_values >= d) for d in bins])
    evaluated_dvh = np.asarray([np.mean(evaluated_values >= d) for d in bins])
    return float(abs(np.trapz(evaluated_dvh, bins) - np.trapz(reference_dvh, bins)))

compare_mean_oar_dvh_auc ¶

compare_mean_oar_dvh_auc(reference: Dose, evaluated: Dose, oars: Sequence[Structure], num_bins: int = 100) -> float

Average the per-OAR absolute AUC differences.

Source code in src/dosemetrics/metrics/dvh.py

def compare_mean_oar_dvh_auc(
    reference: Dose,
    evaluated: Dose,
    oars: Sequence[Structure],
    num_bins: int = 100,
) -> float:
    """Average the per-OAR absolute AUC differences."""
    if not oars:
        raise ValueError("At least one OAR is required")
    return float(
        np.mean(
            [
                compare_oar_dvh_auc(reference, evaluated, oar, num_bins=num_bins)
                for oar in oars
            ]
        )
    )

compare_dvh_score ¶

compare_dvh_score(reference: Dose, evaluated: Dose, targets: Union[Structure, Sequence[Structure]], oars: Sequence[Structure] = ()) -> float

Compare plans using the complete OpenKBP DVH score.

The score is the unweighted mean absolute error over D1, D95, and D99 for every target and mean dose and D0.1cc for every OAR. It is reported in Gy.

Source code in src/dosemetrics/metrics/dvh.py

def compare_dvh_score(
    reference: Dose,
    evaluated: Dose,
    targets: Union[Structure, Sequence[Structure]],
    oars: Sequence[Structure] = (),
) -> float:
    """Compare plans using the complete OpenKBP DVH score.

    The score is the unweighted mean absolute error over D1, D95, and D99 for
    every target and mean dose and D0.1cc for every OAR. It is reported in Gy.
    """
    _validate_comparable_doses(reference, evaluated)
    target_structures = (targets,) if isinstance(targets, Structure) else tuple(targets)
    errors = []

    for target in target_structures:
        for volume_percent in (1.0, 95.0, 99.0):
            reference_value = compute_dose_at_volume(reference, target, volume_percent)
            evaluated_value = compute_dose_at_volume(evaluated, target, volume_percent)
            errors.append(abs(evaluated_value - reference_value))

    for oar in oars:
        errors.append(
            abs(compute_mean_dose(evaluated, oar) - compute_mean_dose(reference, oar))
        )
        errors.append(
            abs(
                compute_dose_at_volume_cc(evaluated, oar, 0.1)
                - compute_dose_at_volume_cc(reference, oar, 0.1)
            )
        )

    if not errors:
        raise ValueError("At least one target or OAR criterion is required")
    return float(np.mean(errors))

The dosemetrics.metrics.comparison module remains importable for backward compatibility, but direct imports are the documented invocation style.

Comparison metadata¶

comparison ¶

Metrics that compare an evaluated radiotherapy plan with a reference plan.

Every public function in this module accepts reference before evaluated. Single-plan quantities remain in their clinical domain modules (dvh, conformity, and homogeneity); this module contains only reference-based comparisons.

Attributes¶

COMPARISON_METRICS `module-attribute` ¶

COMPARISON_METRICS: Tuple[MetricDefinition, ...] = (MetricDefinition('DVH Score', 'compare_dvh_score', MetricCategory.GLOBAL, (EvaluationTask.DOSE_PREDICTION,), 'Gy'), MetricDefinition('Root Mean Squared Error', 'compare_body_rmse', MetricCategory.VOXEL_BASED, (EvaluationTask.DOSE_CALCULATION,), 'Gy'), MetricDefinition('Gamma Index Passing Rate', 'compare_gamma', MetricCategory.VOXEL_BASED, (EvaluationTask.DOSE_CALCULATION,), '%', lower_is_better=False), MetricDefinition('PTV Dose Distance', 'compare_ptv_dose', MetricCategory.PTV_COVERAGE, _BOTH_TASKS, 'Gy'), MetricDefinition('Paddick Conformity Index Distance', 'compare_paddick_conformity_index', MetricCategory.PTV_COVERAGE, _BOTH_TASKS, 'dimensionless'), MetricDefinition('Homogeneity Index Distance', 'compare_homogeneity_index', MetricCategory.PTV_HOMOGENEITY, (EvaluationTask.DOSE_CALCULATION,), 'dimensionless'), MetricDefinition('Paddick Gradient Index Distance', 'compare_paddick_gradient_index', MetricCategory.OAR_SPARING, _BOTH_TASKS, 'dimensionless'), MetricDefinition('OAR Constraint Disagreement', 'compare_oar_constraints', MetricCategory.OAR_SPARING, _BOTH_TASKS, 'fraction'), MetricDefinition('OAR DVH Area Between Curves', 'compare_oar_dvh_auc', MetricCategory.OAR_SPARING, _BOTH_TASKS, 'Gy'))

Classes¶

EvaluationTask ¶

Bases: str, Enum

Dose-estimation task for which a metric is intended.

MetricCategory ¶

Bases: str, Enum

Clinical category used by the plan-comparison metric catalogue.

MetricDefinition `dataclass` ¶

MetricDefinition(name: str, function: str, category: MetricCategory, tasks: Tuple[EvaluationTask, ...], unit: str, lower_is_better: bool = True)

Metadata for one implemented plan-comparison metric.

`conformity`¶

conformity ¶

Conformity indices for target coverage evaluation.

This module provides various conformity indices used to evaluate how well the prescription isodose conforms to the target volume. These metrics are critical for assessing treatment plan quality.

Classes¶

Functions:¶

compute_conformity_index ¶

compute_conformity_index(dose: Dose, target: Structure, prescription_dose: float) -> float

Compute Conformity Index (CI).

CI = V_target_rx / V_rx

Where: - V_target_rx = volume of target receiving >= prescription dose - V_rx = total volume receiving >= prescription dose

Measures how well the prescription isodose conforms to the target. Ideal value is 1.0. Values < 1.0 indicate dose spillage outside target.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`target`	`Structure`	Target structure (PTV, CTV, etc.)	required
`prescription_dose`	`float`	Prescription dose in Gy	required

Returns:

Type	Description
`float`	Conformity index (dimensionless, typically 0-1)

References

ICRU Report 62 (1999)

Examples:

>>> ci = compute_conformity_index(dose, ptv, prescription_dose=60.0)
>>> print(f"Conformity Index: {ci:.3f}")

Source code in src/dosemetrics/metrics/conformity.py

def compute_conformity_index(
    dose: Dose, target: Structure, prescription_dose: float
) -> float:
    """
    Compute Conformity Index (CI).

    CI = V_target_rx / V_rx

    Where:
    - V_target_rx = volume of target receiving >= prescription dose
    - V_rx = total volume receiving >= prescription dose

    Measures how well the prescription isodose conforms to the target.
    Ideal value is 1.0. Values < 1.0 indicate dose spillage outside target.

    Args:
        dose: Dose distribution object
        target: Target structure (PTV, CTV, etc.)
        prescription_dose: Prescription dose in Gy

    Returns:
        Conformity index (dimensionless, typically 0-1)

    References:
        ICRU Report 62 (1999)

    Examples:
        >>> ci = compute_conformity_index(dose, ptv, prescription_dose=60.0)
        >>> print(f"Conformity Index: {ci:.3f}")
    """
    # Volume of target receiving >= prescription dose
    target_dose_values = dose.get_dose_in_structure(target)
    v_target_rx = np.sum(target_dose_values >= prescription_dose)

    # Total volume receiving >= prescription dose
    v_rx = np.sum(dose.dose_array >= prescription_dose)

    if v_rx == 0:
        return 0.0

    return float(v_target_rx / v_rx)

compute_conformity_number ¶

compute_conformity_number(dose: Dose, target: Structure, prescription_dose: float) -> float

Compute Conformity Number (CN) or Conformation Number.

CN = (V_target_rx / V_target) * (V_target_rx / V_rx)

Combines target coverage and dose spillage into a single metric. Ideal value is 1.0.

The first factor (V_target_rx / V_target) represents target coverage. The second factor (V_target_rx / V_rx) represents conformity.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`target`	`Structure`	Target structure	required
`prescription_dose`	`float`	Prescription dose in Gy	required

Returns:

Type	Description
`float`	Conformity number (0-1)

References

van't Riet et al., Int J Radiat Oncol Biol Phys 1997

Examples:

>>> cn = compute_conformity_number(dose, ptv, prescription_dose=60.0)
>>> print(f"Conformity Number: {cn:.3f}")

Source code in src/dosemetrics/metrics/conformity.py

def compute_conformity_number(
    dose: Dose, target: Structure, prescription_dose: float
) -> float:
    """
    Compute Conformity Number (CN) or Conformation Number.

    CN = (V_target_rx / V_target) * (V_target_rx / V_rx)

    Combines target coverage and dose spillage into a single metric.
    Ideal value is 1.0.

    The first factor (V_target_rx / V_target) represents target coverage.
    The second factor (V_target_rx / V_rx) represents conformity.

    Args:
        dose: Dose distribution object
        target: Target structure
        prescription_dose: Prescription dose in Gy

    Returns:
        Conformity number (0-1)

    References:
        van't Riet et al., Int J Radiat Oncol Biol Phys 1997

    Examples:
        >>> cn = compute_conformity_number(dose, ptv, prescription_dose=60.0)
        >>> print(f"Conformity Number: {cn:.3f}")
    """
    target_dose_values = dose.get_dose_in_structure(target)

    v_target = len(target_dose_values)
    if v_target == 0:
        return 0.0

    v_target_rx = np.sum(target_dose_values >= prescription_dose)
    v_rx = np.sum(dose.dose_array >= prescription_dose)

    if v_rx == 0:
        return 0.0

    coverage = v_target_rx / v_target
    conformity = v_target_rx / v_rx

    return float(coverage * conformity)

compute_paddick_conformity_index ¶

compute_paddick_conformity_index(dose: Dose, target: Structure, prescription_dose: float) -> float

Compute Paddick Conformity Index (CI_Paddick).

CI_Paddick = (V_target_rx)^2 / (V_target * V_rx)

This index is commonly used for radiosurgery and SBRT plans. Ideal value is 1.0.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`target`	`Structure`	Target structure	required
`prescription_dose`	`float`	Prescription dose in Gy	required

Returns:

Type	Description
`float`	Paddick conformity index (0-1)

References

Paddick, J Neurosurg 2000

Examples:

>>> # Often used for stereotactic radiosurgery
>>> ci_paddick = compute_paddick_conformity_index(dose, gtv, prescription_dose=18.0)
>>> print(f"Paddick CI: {ci_paddick:.3f}")

Source code in src/dosemetrics/metrics/conformity.py

def compute_paddick_conformity_index(
    dose: Dose, target: Structure, prescription_dose: float
) -> float:
    """
    Compute Paddick Conformity Index (CI_Paddick).

    CI_Paddick = (V_target_rx)^2 / (V_target * V_rx)

    This index is commonly used for radiosurgery and SBRT plans.
    Ideal value is 1.0.

    Args:
        dose: Dose distribution object
        target: Target structure
        prescription_dose: Prescription dose in Gy

    Returns:
        Paddick conformity index (0-1)

    References:
        Paddick, J Neurosurg 2000

    Examples:
        >>> # Often used for stereotactic radiosurgery
        >>> ci_paddick = compute_paddick_conformity_index(dose, gtv, prescription_dose=18.0)
        >>> print(f"Paddick CI: {ci_paddick:.3f}")
    """
    target_dose_values = dose.get_dose_in_structure(target)

    v_target = len(target_dose_values)
    if v_target == 0:
        return 0.0

    v_target_rx = np.sum(target_dose_values >= prescription_dose)
    v_rx = np.sum(dose.dose_array >= prescription_dose)

    if v_rx == 0 or v_target == 0:
        return 0.0

    return float((v_target_rx**2) / (v_target * v_rx))

compute_coverage ¶

compute_coverage(dose: Dose, target: Structure, prescription_dose: float) -> float

Compute target coverage.

Coverage = V_target_rx / V_target

Percentage of target volume receiving at least the prescription dose.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`target`	`Structure`	Target structure	required
`prescription_dose`	`float`	Prescription dose in Gy	required

Returns:

Type	Description
`float`	Coverage as fraction (0-1) or percentage if multiplied by 100

Examples:

>>> coverage = compute_coverage(dose, ptv, prescription_dose=60.0)
>>> print(f"Target coverage: {coverage*100:.1f}%")

Source code in src/dosemetrics/metrics/conformity.py

def compute_coverage(dose: Dose, target: Structure, prescription_dose: float) -> float:
    """
    Compute target coverage.

    Coverage = V_target_rx / V_target

    Percentage of target volume receiving at least the prescription dose.

    Args:
        dose: Dose distribution object
        target: Target structure
        prescription_dose: Prescription dose in Gy

    Returns:
        Coverage as fraction (0-1) or percentage if multiplied by 100

    Examples:
        >>> coverage = compute_coverage(dose, ptv, prescription_dose=60.0)
        >>> print(f"Target coverage: {coverage*100:.1f}%")
    """
    target_dose_values = dose.get_dose_in_structure(target)

    v_target = len(target_dose_values)
    if v_target == 0:
        return 0.0

    v_target_rx = np.sum(target_dose_values >= prescription_dose)

    return float(v_target_rx / v_target)

compute_spillage ¶

compute_spillage(dose: Dose, target: Structure, prescription_dose: float) -> float

Compute dose spillage outside target.

Spillage = (V_rx - V_target_rx) / V_rx

Fraction of prescription isodose volume that is outside the target. Lower values indicate better conformity.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`target`	`Structure`	Target structure	required
`prescription_dose`	`float`	Prescription dose in Gy	required

Returns:

Type	Description
`float`	Spillage as fraction (0-1)

Examples:

>>> spillage = compute_spillage(dose, ptv, prescription_dose=60.0)
>>> print(f"Dose spillage: {spillage*100:.1f}%")

Source code in src/dosemetrics/metrics/conformity.py

def compute_spillage(dose: Dose, target: Structure, prescription_dose: float) -> float:
    """
    Compute dose spillage outside target.

    Spillage = (V_rx - V_target_rx) / V_rx

    Fraction of prescription isodose volume that is outside the target.
    Lower values indicate better conformity.

    Args:
        dose: Dose distribution object
        target: Target structure
        prescription_dose: Prescription dose in Gy

    Returns:
        Spillage as fraction (0-1)

    Examples:
        >>> spillage = compute_spillage(dose, ptv, prescription_dose=60.0)
        >>> print(f"Dose spillage: {spillage*100:.1f}%")
    """
    target_dose_values = dose.get_dose_in_structure(target)
    v_target_rx = np.sum(target_dose_values >= prescription_dose)
    v_rx = np.sum(dose.dose_array >= prescription_dose)

    if v_rx == 0:
        return 0.0

    return float((v_rx - v_target_rx) / v_rx)

compute_rtog_conformity_index ¶

compute_rtog_conformity_index(dose: Dose, target: Structure, prescription_dose: float) -> float

Compute the RTOG Conformity Index (RTOG CI).

RTOG CI = V_Rx / V_target

Where: - V_Rx = total volume receiving >= prescription dose (prescription isodose volume) - V_target = target structure volume

The RTOG CI measures how well the prescription isodose conforms to the target. Values close to 1.0 are ideal. Values > 1.0 indicate over-coverage (dose spillage); values < 1.0 indicate under-coverage.

This differs from the ICRU-based CI in this library (V_target_rx / V_rx), which measures how much of the prescription isodose overlaps the target.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`target`	`Structure`	Target structure (PTV, CTV, etc.)	required
`prescription_dose`	`float`	Prescription dose in Gy	required

Returns:

Type	Description
`float`	RTOG Conformity Index (dimensionless). Ideal value: 1.0.

References

Shaw E, et al. Int J Radiat Oncol Biol Phys. 1993;27(5):1231-9. RTOG 90-05 stereotactic radiosurgery protocol.

Examples:

>>> rtog_ci = compute_rtog_conformity_index(dose, ptv, prescription_dose=60.0)
>>> if 0.9 <= rtog_ci <= 1.1:
...     print("Excellent conformity (RTOG criteria)")
>>> elif 0.7 <= rtog_ci <= 1.5:
...     print("Acceptable conformity (RTOG criteria)")

Source code in src/dosemetrics/metrics/conformity.py

def compute_rtog_conformity_index(
    dose: Dose,
    target: Structure,
    prescription_dose: float,
) -> float:
    """
    Compute the RTOG Conformity Index (RTOG CI).

    RTOG CI = V_Rx / V_target

    Where:
    - V_Rx = total volume receiving >= prescription dose (prescription isodose volume)
    - V_target = target structure volume

    The RTOG CI measures how well the prescription isodose conforms to the target.
    Values close to 1.0 are ideal. Values > 1.0 indicate over-coverage (dose spillage);
    values < 1.0 indicate under-coverage.

    This differs from the ICRU-based CI in this library (V_target_rx / V_rx), which
    measures how much of the prescription isodose overlaps the target.

    Args:
        dose: Dose distribution object
        target: Target structure (PTV, CTV, etc.)
        prescription_dose: Prescription dose in Gy

    Returns:
        RTOG Conformity Index (dimensionless). Ideal value: 1.0.

    References:
        Shaw E, et al. Int J Radiat Oncol Biol Phys. 1993;27(5):1231-9.
        RTOG 90-05 stereotactic radiosurgery protocol.

    Examples:
        >>> rtog_ci = compute_rtog_conformity_index(dose, ptv, prescription_dose=60.0)
        >>> if 0.9 <= rtog_ci <= 1.1:
        ...     print("Excellent conformity (RTOG criteria)")
        >>> elif 0.7 <= rtog_ci <= 1.5:
        ...     print("Acceptable conformity (RTOG criteria)")
    """
    v_rx = int(np.sum(dose.dose_array >= prescription_dose))
    v_target = int(np.sum(target.mask))

    if v_target == 0:
        return float("nan")

    return float(v_rx / v_target)

compute_prescription_mae ¶

compute_prescription_mae(dose: Dose, target: Structure, prescription_dose: float) -> float

Compute the Mean Absolute Error (MAE) between actual dose and prescription dose within target.

Prescription MAE = mean(|dose_in_target - prescription_dose|)

This metric measures how well the dose within the target matches the prescription. A value of 0.0 means every voxel in the target received exactly the prescription dose. Useful for quantifying underdosing and overdosing within the target volume.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`target`	`Structure`	Target structure (PTV, CTV, etc.)	required
`prescription_dose`	`float`	Prescription dose in Gy	required

Returns:

Type	Description
`float`	Mean absolute error from prescription dose in Gy

References

Adapted from PTVPrescriptionMAE in GDP-HMM AAPM Challenge evaluation.

Examples:

>>> mae = compute_prescription_mae(dose, ptv, prescription_dose=60.0)
>>> print(f"Prescription MAE: {mae:.2f} Gy ({mae/60.0*100:.1f}% of prescription)")

Source code in src/dosemetrics/metrics/conformity.py

def compute_prescription_mae(
    dose: Dose,
    target: Structure,
    prescription_dose: float,
) -> float:
    """
    Compute the Mean Absolute Error (MAE) between actual dose and prescription dose within target.

    Prescription MAE = mean(|dose_in_target - prescription_dose|)

    This metric measures how well the dose within the target matches the prescription.
    A value of 0.0 means every voxel in the target received exactly the prescription dose.
    Useful for quantifying underdosing and overdosing within the target volume.

    Args:
        dose: Dose distribution object
        target: Target structure (PTV, CTV, etc.)
        prescription_dose: Prescription dose in Gy

    Returns:
        Mean absolute error from prescription dose in Gy

    References:
        Adapted from PTVPrescriptionMAE in GDP-HMM AAPM Challenge evaluation.

    Examples:
        >>> mae = compute_prescription_mae(dose, ptv, prescription_dose=60.0)
        >>> print(f"Prescription MAE: {mae:.2f} Gy ({mae/60.0*100:.1f}% of prescription)")
    """
    dose_values = dose.get_dose_in_structure(target)

    if len(dose_values) == 0:
        return float("nan")

    return float(np.mean(np.abs(dose_values - prescription_dose)))

`dose_comparison`¶

dose_comparison ¶

Dose distribution comparison metrics beyond DVH.

This module provides image-based metrics for comparing 3D dose distributions. All compare_* functions accept reference before evaluated; compute_* functions characterize one dose distribution.

Future Implementation TODOs

Structural Similarity Index (SSIM) for dose volumes
Mean Squared Error (MSE) and variants
Peak Signal-to-Noise Ratio (PSNR)
Mutual Information
Normalized Cross-Correlation
Dose-volume histogram difference maps

Classes¶

Functions:¶

compare_ssim ¶

compare_ssim(reference: Dose, evaluated: Dose, structure: Optional[Structure] = None, window_size: int = 11, k1: float = 0.01, k2: float = 0.03) -> float

Compute Structural Similarity Index (SSIM) between two dose distributions.

SSIM is a perceptual metric that quantifies image quality degradation based on luminance, contrast, and structure. Originally developed for image comparison, it's applicable to dose distributions.

Parameters¶

reference : Dose Reference dose distribution. evaluated : Dose Evaluated dose distribution. structure : Structure, optional If provided, compute SSIM only within structure volume. If None, compute for entire dose grid. window_size : int, optional Size of sliding window for local SSIM computation (default: 11). k1 : float, optional Algorithm parameter (default: 0.01). k2 : float, optional Algorithm parameter (default: 0.03).

Returns¶

ssim : float Mean SSIM value (0-1, where 1 is perfect similarity).

Notes¶

SSIM ranges from -1 to 1: - 1: Perfect similarity - 0: No structural similarity - -1: Perfect anti-correlation

SSIM considers three components

Luminance: Compares mean intensities
Contrast: Compares standard deviations
Structure: Compares correlation

References¶

Wang Z, Bovik AC, Sheikh HR, Simoncelli EP. "Image quality assessment: from error visibility to structural similarity." IEEE Trans Image Process. 2004;13(4):600-12.

Examples¶

ssim = compare_ssim(planned_dose, delivered_dose, ptv) print(f"Dose SSIM: {ssim:.3f}") if ssim > 0.95: ... print("Excellent agreement")

Raises¶

NotImplementedError This function is a stub for future implementation. ValueError If dose distributions have incompatible geometry.

Source code in src/dosemetrics/metrics/dose_comparison.py

def compare_ssim(
    reference: Dose,
    evaluated: Dose,
    structure: Optional[Structure] = None,
    window_size: int = 11,
    k1: float = 0.01,
    k2: float = 0.03,
) -> float:
    """
    Compute Structural Similarity Index (SSIM) between two dose distributions.

    SSIM is a perceptual metric that quantifies image quality degradation
    based on luminance, contrast, and structure. Originally developed for
    image comparison, it's applicable to dose distributions.

    Parameters
    ----------
    reference : Dose
        Reference dose distribution.
    evaluated : Dose
        Evaluated dose distribution.
    structure : Structure, optional
        If provided, compute SSIM only within structure volume.
        If None, compute for entire dose grid.
    window_size : int, optional
        Size of sliding window for local SSIM computation (default: 11).
    k1 : float, optional
        Algorithm parameter (default: 0.01).
    k2 : float, optional
        Algorithm parameter (default: 0.03).

    Returns
    -------
    ssim : float
        Mean SSIM value (0-1, where 1 is perfect similarity).

    Notes
    -----
    SSIM ranges from -1 to 1:
        - 1: Perfect similarity
        - 0: No structural similarity
        - -1: Perfect anti-correlation

    SSIM considers three components:
        - Luminance: Compares mean intensities
        - Contrast: Compares standard deviations
        - Structure: Compares correlation

    References
    ----------
    - Wang Z, Bovik AC, Sheikh HR, Simoncelli EP. "Image quality assessment:
      from error visibility to structural similarity." IEEE Trans Image Process.
      2004;13(4):600-12.

    Examples
    --------
    >>> ssim = compare_ssim(planned_dose, delivered_dose, ptv)
    >>> print(f"Dose SSIM: {ssim:.3f}")
    >>> if ssim > 0.95:
    ...     print("Excellent agreement")

    Raises
    ------
    NotImplementedError
        This function is a stub for future implementation.
    ValueError
        If dose distributions have incompatible geometry.
    """
    # Get dose arrays
    arr1 = reference.dose_array
    arr2 = evaluated.dose_array

    # Check shapes match
    if arr1.shape != arr2.shape:
        raise ValueError(f"Dose shapes must match: {arr1.shape} vs {arr2.shape}")

    # Apply structure mask if provided
    if structure is not None:
        mask = structure.mask
        # For 3D SSIM, we need to work with the full volume
        # but we'll compute SSIM and then weight by the mask
        arr1_masked = np.where(mask, arr1, 0)
        arr2_masked = np.where(mask, arr2, 0)
    else:
        arr1_masked = arr1
        arr2_masked = arr2

    # Compute SSIM for 3D volume
    # Use smaller window for medical images
    win_size = min(window_size, min(arr1.shape) - 1)
    if win_size % 2 == 0:
        win_size -= 1  # Must be odd
    win_size = max(3, win_size)  # At least 3

    data_range = max(np.max(arr1), np.max(arr2))

    try:
        ssim_value = structural_similarity(
            arr1_masked,
            arr2_masked,
            data_range=data_range,
            win_size=win_size,
            K1=k1,
            K2=k2,
        )
    except ValueError:
        # If window size is too large, reduce it
        win_size = 3
        ssim_value = structural_similarity(
            arr1_masked,
            arr2_masked,
            data_range=data_range,
            win_size=win_size,
            K1=k1,
            K2=k2,
        )

    return float(ssim_value)

compare_mse ¶

compare_mse(reference: Dose, evaluated: Dose, structure: Optional[Structure] = None) -> float

Compute Mean Squared Error between two dose distributions.

Parameters¶

reference : Dose Reference dose. evaluated : Dose Evaluated dose. structure : Structure, optional If provided, compute MSE only within structure.

Returns¶

mse : float Mean squared error in Gy^2.

Raises¶

ValueError If dose distributions have incompatible shapes.

Source code in src/dosemetrics/metrics/dose_comparison.py

def compare_mse(
    reference: Dose, evaluated: Dose, structure: Optional[Structure] = None
) -> float:
    """
    Compute Mean Squared Error between two dose distributions.

    Parameters
    ----------
    reference : Dose
        Reference dose.
    evaluated : Dose
        Evaluated dose.
    structure : Structure, optional
        If provided, compute MSE only within structure.

    Returns
    -------
    mse : float
        Mean squared error in Gy^2.

    Raises
    ------
    ValueError
        If dose distributions have incompatible shapes.
    """
    # Get dose arrays
    arr1 = reference.dose_array
    arr2 = evaluated.dose_array

    # Check shapes match
    if arr1.shape != arr2.shape:
        raise ValueError(f"Dose shapes must match: {arr1.shape} vs {arr2.shape}")

    # Apply structure mask if provided
    if structure is not None:
        mask = structure.mask
        arr1 = arr1[mask]
        arr2 = arr2[mask]

    # Compute MSE
    mse = np.mean((arr1 - arr2) ** 2)
    return float(mse)

compare_mae ¶

compare_mae(reference: Dose, evaluated: Dose, structure: Optional[Structure] = None) -> float

Compute Mean Absolute Error between two dose distributions.

Parameters¶

reference : Dose Reference dose. evaluated : Dose Evaluated dose. structure : Structure, optional If provided, compute MAE only within structure.

Returns¶

mae : float Mean absolute error in Gy.

Notes¶

MAE is often more interpretable than MSE for dose comparison as it's in the same units as dose (Gy).

Raises¶

ValueError If dose distributions have incompatible shapes.

Source code in src/dosemetrics/metrics/dose_comparison.py

def compare_mae(
    reference: Dose, evaluated: Dose, structure: Optional[Structure] = None
) -> float:
    """
    Compute Mean Absolute Error between two dose distributions.

    Parameters
    ----------
    reference : Dose
        Reference dose.
    evaluated : Dose
        Evaluated dose.
    structure : Structure, optional
        If provided, compute MAE only within structure.

    Returns
    -------
    mae : float
        Mean absolute error in Gy.

    Notes
    -----
    MAE is often more interpretable than MSE for dose comparison as it's
    in the same units as dose (Gy).

    Raises
    ------
    ValueError
        If dose distributions have incompatible shapes.
    """
    # Get dose arrays
    arr1 = reference.dose_array
    arr2 = evaluated.dose_array

    # Check shapes match
    if arr1.shape != arr2.shape:
        raise ValueError(f"Dose shapes must match: {arr1.shape} vs {arr2.shape}")

    # Apply structure mask if provided
    if structure is not None:
        mask = structure.mask
        arr1 = arr1[mask]
        arr2 = arr2[mask]

    # Compute MAE
    mae = np.mean(np.abs(arr1 - arr2))
    return float(mae)

compare_psnr ¶

compare_psnr(reference: Dose, evaluated: Dose, structure: Optional[Structure] = None, data_range: Optional[float] = None) -> float

Compute Peak Signal-to-Noise Ratio between two dose distributions.

Parameters¶

reference : Dose Reference dose. evaluated : Dose Evaluated dose. structure : Structure, optional If provided, compute PSNR only within structure. data_range : float, optional Data range (max - min). If None, computed from doses.

Returns¶

psnr : float Peak signal-to-noise ratio in dB.

Notes¶

PSNR is defined as: PSNR = 10 * log10((MAX^2) / MSE) Higher values indicate better similarity.

Raises¶

ValueError If dose distributions have incompatible shapes or MSE is zero.

Source code in src/dosemetrics/metrics/dose_comparison.py

def compare_psnr(
    reference: Dose,
    evaluated: Dose,
    structure: Optional[Structure] = None,
    data_range: Optional[float] = None,
) -> float:
    """
    Compute Peak Signal-to-Noise Ratio between two dose distributions.

    Parameters
    ----------
    reference : Dose
        Reference dose.
    evaluated : Dose
        Evaluated dose.
    structure : Structure, optional
        If provided, compute PSNR only within structure.
    data_range : float, optional
        Data range (max - min). If None, computed from doses.

    Returns
    -------
    psnr : float
        Peak signal-to-noise ratio in dB.

    Notes
    -----
    PSNR is defined as: PSNR = 10 * log10((MAX^2) / MSE)
    Higher values indicate better similarity.

    Raises
    ------
    ValueError
        If dose distributions have incompatible shapes or MSE is zero.
    """
    # Compute MSE
    mse = compare_mse(reference, evaluated, structure)

    if mse == 0:
        return float("inf")  # Perfect match

    # Determine data range
    if data_range is None:
        arr1 = reference.dose_array
        arr2 = evaluated.dose_array
        if structure is not None:
            mask = structure.mask
            arr1 = arr1[mask]
            arr2 = arr2[mask]
        data_range = max(np.max(arr1), np.max(arr2))

    # Compute PSNR
    psnr = 10 * np.log10((data_range**2) / mse)
    return float(psnr)

compare_mutual_information ¶

compare_mutual_information(reference: Dose, evaluated: Dose, structure: Optional[Structure] = None, bins: int = 256) -> float

Compute Mutual Information between two dose distributions.

Parameters¶

reference : Dose Reference dose distribution. evaluated : Dose Evaluated dose distribution. structure : Structure, optional If provided, compute MI only within structure. bins : int, optional Number of histogram bins (default: 256).

Returns¶

mi : float Mutual information value (higher indicates more similarity).

Notes¶

Mutual Information quantifies the information shared between two distributions. It's particularly useful for multimodal comparison.

Raises¶

ValueError If dose distributions have incompatible shapes.

Source code in src/dosemetrics/metrics/dose_comparison.py

def compare_mutual_information(
    reference: Dose,
    evaluated: Dose,
    structure: Optional[Structure] = None,
    bins: int = 256,
) -> float:
    """
    Compute Mutual Information between two dose distributions.

    Parameters
    ----------
    reference : Dose
        Reference dose distribution.
    evaluated : Dose
        Evaluated dose distribution.
    structure : Structure, optional
        If provided, compute MI only within structure.
    bins : int, optional
        Number of histogram bins (default: 256).

    Returns
    -------
    mi : float
        Mutual information value (higher indicates more similarity).

    Notes
    -----
    Mutual Information quantifies the information shared between two
    distributions. It's particularly useful for multimodal comparison.

    Raises
    ------
    ValueError
        If dose distributions have incompatible shapes.
    """
    # Get dose arrays
    arr1 = reference.dose_array.flatten()
    arr2 = evaluated.dose_array.flatten()

    # Check shapes match
    if arr1.shape != arr2.shape:
        raise ValueError("Dose shapes must match")

    # Apply structure mask if provided
    if structure is not None:
        mask = structure.mask.flatten()
        arr1 = arr1[mask]
        arr2 = arr2[mask]

    # Compute 2D histogram
    hist_2d, x_edges, y_edges = np.histogram2d(arr1, arr2, bins=bins)

    # Add small epsilon to avoid log(0)
    hist_2d = hist_2d + np.finfo(float).eps

    # Normalize to get joint probability
    pxy = hist_2d / np.sum(hist_2d)

    # Compute marginal probabilities
    px = np.sum(pxy, axis=1)
    py = np.sum(pxy, axis=0)

    # Compute mutual information
    # MI = sum(p(x,y) * log(p(x,y) / (p(x) * p(y))))
    px_py = px[:, None] * py[None, :]

    # Only compute where both are non-zero
    nonzero = (pxy > 0) & (px_py > 0)
    mi = np.sum(pxy[nonzero] * np.log(pxy[nonzero] / px_py[nonzero]))

    return float(mi)

compare_normalized_cross_correlation ¶

compare_normalized_cross_correlation(reference: Dose, evaluated: Dose, structure: Optional[Structure] = None) -> float

Compute Normalized Cross-Correlation between two dose distributions.

Parameters¶

reference : Dose Reference dose distribution. evaluated : Dose Evaluated dose distribution. structure : Structure, optional If provided, compute NCC only within structure.

Returns¶

ncc : float Normalized cross-correlation (-1 to 1).

Notes¶

NCC is Pearson correlation coefficient for images/volumes. Values close to 1 indicate high positive correlation.

Raises¶

ValueError If dose distributions have incompatible shapes.

Source code in src/dosemetrics/metrics/dose_comparison.py

def compare_normalized_cross_correlation(
    reference: Dose, evaluated: Dose, structure: Optional[Structure] = None
) -> float:
    """
    Compute Normalized Cross-Correlation between two dose distributions.

    Parameters
    ----------
    reference : Dose
        Reference dose distribution.
    evaluated : Dose
        Evaluated dose distribution.
    structure : Structure, optional
        If provided, compute NCC only within structure.

    Returns
    -------
    ncc : float
        Normalized cross-correlation (-1 to 1).

    Notes
    -----
    NCC is Pearson correlation coefficient for images/volumes.
    Values close to 1 indicate high positive correlation.

    Raises
    ------
    ValueError
        If dose distributions have incompatible shapes.
    """
    # Get dose arrays
    arr1 = reference.dose_array.flatten()
    arr2 = evaluated.dose_array.flatten()

    # Check shapes match
    if arr1.shape != arr2.shape:
        raise ValueError("Dose shapes must match")

    # Apply structure mask if provided
    if structure is not None:
        mask = structure.mask.flatten()
        arr1 = arr1[mask]
        arr2 = arr2[mask]

    # Compute NCC (Pearson correlation)
    # NCC = sum((x - mean_x) * (y - mean_y)) / (std_x * std_y * N)
    mean1 = np.mean(arr1)
    mean2 = np.mean(arr2)

    numerator = np.sum((arr1 - mean1) * (arr2 - mean2))
    denominator = np.sqrt(np.sum((arr1 - mean1) ** 2) * np.sum((arr2 - mean2) ** 2))

    if denominator == 0:
        return 0.0  # No variation in one or both images

    ncc = numerator / denominator
    return float(ncc)

compare_dose_difference_map ¶

compare_dose_difference_map(reference: Dose, evaluated: Dose, absolute: bool = False) -> Dose

Compute voxel-wise dose difference map.

Parameters¶

reference : Dose Reference dose. evaluated : Dose Evaluated dose. absolute : bool, optional If True, return absolute differences (default: False).

Returns¶

diff_dose : Dose Dose object containing difference map.

Notes¶

Useful for visualizing spatial dose discrepancies.

Raises¶

ValueError If dose distributions have incompatible shapes.

Source code in src/dosemetrics/metrics/dose_comparison.py

def compare_dose_difference_map(
    reference: Dose, evaluated: Dose, absolute: bool = False
) -> Dose:
    """
    Compute voxel-wise dose difference map.

    Parameters
    ----------
    reference : Dose
        Reference dose.
    evaluated : Dose
        Evaluated dose.
    absolute : bool, optional
        If True, return absolute differences (default: False).

    Returns
    -------
    diff_dose : Dose
        Dose object containing difference map.

    Notes
    -----
    Useful for visualizing spatial dose discrepancies.

    Raises
    ------
    ValueError
        If dose distributions have incompatible shapes.
    """
    # Check shapes match
    if reference.dose_array.shape != evaluated.dose_array.shape:
        raise ValueError(
            f"Dose shapes must match: {reference.dose_array.shape} vs "
            f"{evaluated.dose_array.shape}"
        )

    # Compute difference
    if absolute:
        diff_grid = np.abs(reference.dose_array - evaluated.dose_array)
    else:
        diff_grid = reference.dose_array - evaluated.dose_array

    # Create new Dose object with difference
    diff_dose = Dose(
        dose_array=diff_grid,
        spacing=reference.spacing,
        origin=reference.origin,
        name=f"{reference.name}_diff",
    )

    return diff_dose

compare_dose ¶

compare_dose(reference: Dose, evaluated: Dose, structure: Optional[Structure] = None) -> Dict[str, float]

Compute comprehensive set of dose comparison metrics.

Parameters¶

reference : Dose Reference dose. evaluated : Dose Evaluated dose. structure : Structure, optional If provided, compute metrics only within structure.

Returns¶

metrics : dict Dictionary containing: - 'ssim': Structural similarity index - 'mse': Mean squared error - 'mae': Mean absolute error - 'psnr': Peak signal-to-noise ratio - 'ncc': Normalized cross-correlation - 'mi': Mutual information

Examples¶

metrics = compare_dose(reference, evaluated, ptv) print(f"SSIM: {metrics['ssim']:.3f}") print(f"MAE: {metrics['mae']:.2f} Gy")

Raises¶

ValueError If dose distributions have incompatible shapes.

Source code in src/dosemetrics/metrics/dose_comparison.py

def compare_dose(
    reference: Dose, evaluated: Dose, structure: Optional[Structure] = None
) -> Dict[str, float]:
    """
    Compute comprehensive set of dose comparison metrics.

    Parameters
    ----------
    reference : Dose
        Reference dose.
    evaluated : Dose
        Evaluated dose.
    structure : Structure, optional
        If provided, compute metrics only within structure.

    Returns
    -------
    metrics : dict
        Dictionary containing:
            - 'ssim': Structural similarity index
            - 'mse': Mean squared error
            - 'mae': Mean absolute error
            - 'psnr': Peak signal-to-noise ratio
            - 'ncc': Normalized cross-correlation
            - 'mi': Mutual information

    Examples
    --------
    >>> metrics = compare_dose(reference, evaluated, ptv)
    >>> print(f"SSIM: {metrics['ssim']:.3f}")
    >>> print(f"MAE: {metrics['mae']:.2f} Gy")

    Raises
    ------
    ValueError
        If dose distributions have incompatible shapes.
    """
    metrics = {}

    try:
        metrics["mse"] = compare_mse(reference, evaluated, structure)
    except Exception as e:
        warnings.warn(f"MSE computation failed: {e}")
        metrics["mse"] = np.nan

    try:
        metrics["mae"] = compare_mae(reference, evaluated, structure)
    except Exception as e:
        warnings.warn(f"MAE computation failed: {e}")
        metrics["mae"] = np.nan

    try:
        metrics["psnr"] = compare_psnr(reference, evaluated, structure)
    except Exception as e:
        warnings.warn(f"PSNR computation failed: {e}")
        metrics["psnr"] = np.nan

    try:
        metrics["ssim"] = compare_ssim(reference, evaluated, structure)
    except Exception as e:
        warnings.warn(f"SSIM computation failed: {e}")
        metrics["ssim"] = np.nan

    try:
        metrics["ncc"] = compare_normalized_cross_correlation(
            reference, evaluated, structure
        )
    except Exception as e:
        warnings.warn(f"NCC computation failed: {e}")
        metrics["ncc"] = np.nan

    try:
        metrics["mi"] = compare_mutual_information(reference, evaluated, structure)
    except Exception as e:
        warnings.warn(f"MI computation failed: {e}")
        metrics["mi"] = np.nan

    return metrics

compute_3d_dose_gradient ¶

compute_3d_dose_gradient(dose: Dose) -> Tuple[np.ndarray, np.ndarray, np.ndarray]

Compute 3D dose gradient (useful for dose falloff analysis).

Parameters¶

dose : Dose Dose distribution.

Returns¶

grad_x : np.ndarray Gradient in x direction. grad_y : np.ndarray Gradient in y direction. grad_z : np.ndarray Gradient in z direction.

Notes¶

Uses numpy gradient function which computes central differences in the interior and first differences at the boundaries.

The gradient is useful for analyzing dose falloff regions and identifying high-gradient areas.

Source code in src/dosemetrics/metrics/dose_comparison.py

def compute_3d_dose_gradient(dose: Dose) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Compute 3D dose gradient (useful for dose falloff analysis).

    Parameters
    ----------
    dose : Dose
        Dose distribution.

    Returns
    -------
    grad_x : np.ndarray
        Gradient in x direction.
    grad_y : np.ndarray
        Gradient in y direction.
    grad_z : np.ndarray
        Gradient in z direction.

    Notes
    -----
    Uses numpy gradient function which computes central differences
    in the interior and first differences at the boundaries.

    The gradient is useful for analyzing dose falloff regions and
    identifying high-gradient areas.
    """
    dose_array = dose.dose_array

    # Get voxel spacing from dose object
    spacing = dose.spacing

    # Compute gradients in each direction
    # Note: numpy.gradient returns gradients in the order of axes
    grad_z, grad_y, grad_x = np.gradient(dose_array, spacing[2], spacing[1], spacing[0])

    return grad_x, grad_y, grad_z

compute_variance_of_laplacian ¶

compute_variance_of_laplacian(dose: Dose, structure: Optional[Structure] = None) -> float

Compute the Variance of Laplacian (VoL) as a measure of dose distribution sharpness.

A higher variance indicates sharper, more spatially complex dose gradients (common in modern IMRT/VMAT plans). A lower variance indicates smoother, more homogeneous dose distributions.

The Laplacian operator highlights regions of rapid dose change (edges/interfaces). For 3D volumes, the Laplacian is applied slice-by-slice along the first axis and the variance is averaged across all slices.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure`	`Optional[Structure]`	If provided, only consider voxels within this structure for variance computation. If None, uses the full dose volume.	`None`

Returns:

Type	Description
`float`	Average variance of the Laplacian (dimensionless). Higher = sharper dose edges.

References

Adapted from VarianceOfLaplacian metric in GDP-HMM AAPM Challenge. Laplacian kernel: [[0,1,0],[1,-4,1],[0,1,0]]

Examples:

>>> vol = compute_variance_of_laplacian(dose)
>>> print(f"Dose sharpness (VoL): {vol:.4f}")
>>>
>>> # Compare sharpness within target vs globally
>>> vol_ptv = compute_variance_of_laplacian(dose, ptv)
>>> vol_global = compute_variance_of_laplacian(dose)

Source code in src/dosemetrics/metrics/dose_comparison.py

def compute_variance_of_laplacian(
    dose: Dose,
    structure: Optional[Structure] = None,
) -> float:
    """
    Compute the Variance of Laplacian (VoL) as a measure of dose distribution sharpness.

    A higher variance indicates sharper, more spatially complex dose gradients
    (common in modern IMRT/VMAT plans). A lower variance indicates smoother,
    more homogeneous dose distributions.

    The Laplacian operator highlights regions of rapid dose change (edges/interfaces).
    For 3D volumes, the Laplacian is applied slice-by-slice along the first axis
    and the variance is averaged across all slices.

    Args:
        dose: Dose distribution object
        structure: If provided, only consider voxels within this structure
            for variance computation. If None, uses the full dose volume.

    Returns:
        Average variance of the Laplacian (dimensionless). Higher = sharper dose edges.

    References:
        Adapted from VarianceOfLaplacian metric in GDP-HMM AAPM Challenge.
        Laplacian kernel: [[0,1,0],[1,-4,1],[0,1,0]]

    Examples:
        >>> vol = compute_variance_of_laplacian(dose)
        >>> print(f"Dose sharpness (VoL): {vol:.4f}")
        >>>
        >>> # Compare sharpness within target vs globally
        >>> vol_ptv = compute_variance_of_laplacian(dose, ptv)
        >>> vol_global = compute_variance_of_laplacian(dose)
    """
    from scipy.ndimage import laplace

    dose_array = dose.dose_array.astype(float)

    if structure is not None:
        # Compute VoL only within the bounding box of the structure mask
        mask = structure.mask
        if not np.any(mask):
            return float("nan")
        # Apply mask: set non-structure voxels to mean value to avoid edge artifacts
        mean_val = float(np.mean(dose_array[mask]))
        masked = np.where(mask, dose_array, mean_val)
        laplacian = laplace(masked)
        return float(np.var(laplacian[mask]))

    # Global: apply Laplacian to each axial slice and average variance
    slice_variances = []
    for z in range(dose_array.shape[2]):
        lap_slice = laplace(dose_array[:, :, z])
        slice_variances.append(float(np.var(lap_slice)))

    return float(np.mean(slice_variances)) if slice_variances else float("nan")

compare_normalized_mae ¶

compare_normalized_mae(reference: Dose, evaluated: Dose, structure: Optional[Structure] = None, normalization_value: Optional[float] = None, dose_threshold_gy: Optional[float] = None) -> float

Compute Normalized MAE with optional threshold masking.

Normalized MAE = mean(|dose_ref - dose_eval|) / normalization_value

Optionally restricts computation to voxels where the reference dose exceeds a threshold, focusing the metric on clinically relevant dose regions.

This is a generalization of the GDP-HMM Challenge MAE metric, adapted for use with arbitrary structures and normalization values.

Parameters:

Name	Type	Description	Default
`reference`	`Dose`	Reference dose distribution	required
`evaluated`	`Dose`	Evaluated dose distribution to compare	required
`structure`	`Optional[Structure]`	If provided, restrict computation to this structure. Uses the full dose volume if None.	`None`
`normalization_value`	`Optional[float]`	Value to normalize the MAE by (e.g., prescription dose). If None, returns un-normalized MAE (equivalent to compare_mae).	`None`
`dose_threshold_gy`	`Optional[float]`	If provided, only include voxels where the reference dose exceeds this threshold in Gy. Useful for focusing on high-dose regions and ignoring low-dose areas outside the treatment field.	`None`

Returns:

Type	Description
`float`	Normalized MAE (dimensionless if normalization_value provided, else Gy).
`float`	Returns NaN if no voxels remain after applying the threshold mask.

References

Adapted from ChallengeMAE in GDP-HMM AAPM Challenge evaluation.

Examples:

>>> # Normalized by prescription dose (60 Gy), only high-dose region
>>> nMAE = compare_normalized_mae(
...     reference_dose, predicted_dose,
...     structure=body,
...     normalization_value=60.0,
...     dose_threshold_gy=5.0
... )
>>> print(f"Normalized MAE: {nMAE:.4f}")

Source code in src/dosemetrics/metrics/dose_comparison.py

def compare_normalized_mae(
    reference: Dose,
    evaluated: Dose,
    structure: Optional[Structure] = None,
    normalization_value: Optional[float] = None,
    dose_threshold_gy: Optional[float] = None,
) -> float:
    """
    Compute Normalized MAE with optional threshold masking.

    Normalized MAE = mean(|dose_ref - dose_eval|) / normalization_value

    Optionally restricts computation to voxels where the reference dose exceeds a
    threshold, focusing the metric on clinically relevant dose regions.

    This is a generalization of the GDP-HMM Challenge MAE metric, adapted for
    use with arbitrary structures and normalization values.

    Args:
        reference: Reference dose distribution
        evaluated: Evaluated dose distribution to compare
        structure: If provided, restrict computation to this structure. Uses the
            full dose volume if None.
        normalization_value: Value to normalize the MAE by (e.g., prescription dose).
            If None, returns un-normalized MAE (equivalent to compare_mae).
        dose_threshold_gy: If provided, only include voxels where the reference
            dose exceeds this threshold in Gy. Useful for focusing on high-dose
            regions and ignoring low-dose areas outside the treatment field.

    Returns:
        Normalized MAE (dimensionless if normalization_value provided, else Gy).
        Returns NaN if no voxels remain after applying the threshold mask.

    References:
        Adapted from ChallengeMAE in GDP-HMM AAPM Challenge evaluation.

    Examples:
        >>> # Normalized by prescription dose (60 Gy), only high-dose region
        >>> nMAE = compare_normalized_mae(
        ...     reference_dose, predicted_dose,
        ...     structure=body,
        ...     normalization_value=60.0,
        ...     dose_threshold_gy=5.0
        ... )
        >>> print(f"Normalized MAE: {nMAE:.4f}")
    """
    if structure is not None:
        ref_arr = reference.get_dose_in_structure(structure)
        eval_arr = evaluated.get_dose_in_structure(structure)
    else:
        ref_arr = reference.dose_array.flatten()
        eval_arr = evaluated.dose_array.flatten()

    if len(ref_arr) == 0:
        return float("nan")

    if dose_threshold_gy is not None:
        mask = ref_arr >= dose_threshold_gy
        ref_arr = ref_arr[mask]
        eval_arr = eval_arr[mask]

    if len(ref_arr) == 0:
        return float("nan")

    mae = float(np.mean(np.abs(ref_arr - eval_arr)))

    if normalization_value is not None and normalization_value > 0:
        return mae / normalization_value

    return mae

`dvh`¶

dvh ¶

Dose-volume histogram computation, metrics, and comparisons.

compute_* functions describe one dose plan or a collection of plans. compare_* functions always accept reference before evaluated. Keeping both in this module makes the DVH data model and all operations on it discoverable in one place.

Classes¶

Functions:¶

compute_dvh ¶

compute_dvh(dose: Dose, structure: Structure, max_dose: Optional[float] = None, step_size: float = 0.1, verbose: bool = False) -> Tuple[np.ndarray, np.ndarray]

Compute dose-volume histogram for a structure.

A DVH shows the percentage of structure volume that receives at least a given dose level.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure`	`Structure`	Structure to compute DVH for	required
`max_dose`	`Optional[float]`	Maximum dose for histogram bins (auto-detect if None)	`None`
`step_size`	`float`	Bin width in Gy	`0.1`
`verbose`	`bool`	Print a compact summary when `True` (default: `False`)	`False`

Returns:

Type	Description
`ndarray`	Tuple of (dose_bins, volume_percentages)
`ndarray`	dose_bins: Array of dose levels (Gy)
`Tuple[ndarray, ndarray]`	volume_percentages: Percentage of volume receiving >= each dose (0-100)

Examples:

>>> from dosemetrics.dose import Dose
>>> from dosemetrics.metrics import dvh
>>>
>>> dose = Dose.from_dicom("rtdose.dcm")
>>> ptv = structures.get_structure("PTV")
>>>
>>> dose_bins, volumes = dvh.compute_dvh(dose, ptv, verbose=True)
>>>
>>> # Plot DVH
>>> import matplotlib.pyplot as plt
>>> plt.plot(dose_bins, volumes)
>>> plt.xlabel("Dose (Gy)")
>>> plt.ylabel("Volume (%)")

Source code in src/dosemetrics/metrics/dvh.py

def compute_dvh(
    dose: Dose,
    structure: Structure,
    max_dose: Optional[float] = None,
    step_size: float = 0.1,
    verbose: bool = False,
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Compute dose-volume histogram for a structure.

    A DVH shows the percentage of structure volume that receives at least
    a given dose level.

    Args:
        dose: Dose distribution object
        structure: Structure to compute DVH for
        max_dose: Maximum dose for histogram bins (auto-detect if None)
        step_size: Bin width in Gy
        verbose: Print a compact summary when ``True`` (default: ``False``)

    Returns:
        Tuple of (dose_bins, volume_percentages)
        - dose_bins: Array of dose levels (Gy)
        - volume_percentages: Percentage of volume receiving >= each dose (0-100)

    Examples:
        >>> from dosemetrics.dose import Dose
        >>> from dosemetrics.metrics import dvh
        >>>
        >>> dose = Dose.from_dicom("rtdose.dcm")
        >>> ptv = structures.get_structure("PTV")
        >>>
        >>> dose_bins, volumes = dvh.compute_dvh(dose, ptv, verbose=True)
        >>>
        >>> # Plot DVH
        >>> import matplotlib.pyplot as plt
        >>> plt.plot(dose_bins, volumes)
        >>> plt.xlabel("Dose (Gy)")
        >>> plt.ylabel("Volume (%)")
    """
    dose_values = dose.get_dose_in_structure(structure)

    if len(dose_values) == 0:
        bins = np.array([0.0])
        volumes = np.array([0.0])
        if verbose:
            print(f"{structure.name} DVH: empty structure")
        return bins, volumes

    if max_dose is None:
        max_dose = float(np.max(dose_values))

    bins = np.arange(0, max_dose + step_size, step_size)
    volumes = np.array(
        [
            100.0 * np.sum(dose_values >= dose_bin) / len(dose_values)
            for dose_bin in bins
        ]
    )

    if verbose:
        print(
            f"{structure.name} DVH: {len(bins)} bins, "
            f"{bins[0]:.2f}-{bins[-1]:.2f} Gy, "
            f"{volumes.min():.1f}-{volumes.max():.1f}% volume"
        )

    return bins, volumes

compute_volume_at_dose ¶

compute_volume_at_dose(dose: Dose, structure: Structure, dose_threshold: float) -> float

Compute percentage of structure receiving at least the dose threshold.

This computes VX where X is the dose threshold (e.g., V20 = % volume >= 20 Gy).

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure`	`Structure`	Structure to analyze	required
`dose_threshold`	`float`	Dose threshold in Gy	required

Returns:

Type	Description
`float`	Percentage of volume (0-100) receiving >= dose_threshold

Examples:

>>> # V20: percentage of lung receiving >= 20 Gy
>>> v20 = compute_volume_at_dose(dose, lung, 20.0)
>>> print(f"V20: {v20:.1f}%")
>>>
>>> # V5: percentage of heart receiving >= 5 Gy
>>> v5 = compute_volume_at_dose(dose, heart, 5.0)

Source code in src/dosemetrics/metrics/dvh.py

def compute_volume_at_dose(
    dose: Dose, structure: Structure, dose_threshold: float
) -> float:
    """
    Compute percentage of structure receiving at least the dose threshold.

    This computes VX where X is the dose threshold (e.g., V20 = % volume >= 20 Gy).

    Args:
        dose: Dose distribution object
        structure: Structure to analyze
        dose_threshold: Dose threshold in Gy

    Returns:
        Percentage of volume (0-100) receiving >= dose_threshold

    Examples:
        >>> # V20: percentage of lung receiving >= 20 Gy
        >>> v20 = compute_volume_at_dose(dose, lung, 20.0)
        >>> print(f"V20: {v20:.1f}%")
        >>>
        >>> # V5: percentage of heart receiving >= 5 Gy
        >>> v5 = compute_volume_at_dose(dose, heart, 5.0)
    """
    dose_values = dose.get_dose_in_structure(structure)

    if len(dose_values) == 0:
        return 0.0

    return float(100.0 * np.sum(dose_values >= dose_threshold) / len(dose_values))

compute_dose_at_volume ¶

compute_dose_at_volume(dose: Dose, structure: Structure, volume_percent: float) -> float

Compute dose received by a given percentage of structure volume.

This computes DX where X is the volume percentage (e.g., D95 = dose to 95% of volume).

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure`	`Structure`	Structure to analyze	required
`volume_percent`	`float`	Volume percentage (0-100)	required

Returns:

Type	Description
`float`	Dose in Gy that the specified volume percentage receives

Raises:

Type	Description
`ValueError`	If volume_percent is not in range 0-100

Examples:

>>> # D95: dose covering 95% of PTV
>>> d95 = compute_dose_at_volume(dose, ptv, 95)
>>> print(f"D95: {d95:.2f} Gy")
>>>
>>> # D_0.1cc for OAR (requires volume in cc conversion)
>>> # For now, use percentile approximation
>>> d_max = compute_dose_at_volume(dose, brainstem, 0.1)

Source code in src/dosemetrics/metrics/dvh.py

def compute_dose_at_volume(
    dose: Dose, structure: Structure, volume_percent: float
) -> float:
    """
    Compute dose received by a given percentage of structure volume.

    This computes DX where X is the volume percentage (e.g., D95 = dose to 95% of volume).

    Args:
        dose: Dose distribution object
        structure: Structure to analyze
        volume_percent: Volume percentage (0-100)

    Returns:
        Dose in Gy that the specified volume percentage receives

    Raises:
        ValueError: If volume_percent is not in range 0-100

    Examples:
        >>> # D95: dose covering 95% of PTV
        >>> d95 = compute_dose_at_volume(dose, ptv, 95)
        >>> print(f"D95: {d95:.2f} Gy")
        >>>
        >>> # D_0.1cc for OAR (requires volume in cc conversion)
        >>> # For now, use percentile approximation
        >>> d_max = compute_dose_at_volume(dose, brainstem, 0.1)
    """
    if not 0 <= volume_percent <= 100:
        raise ValueError(f"Volume percent must be 0-100, got {volume_percent}")

    dose_values = dose.get_dose_in_structure(structure)

    if len(dose_values) == 0:
        return 0.0

    # DX means X% of volume receives AT LEAST this dose
    # This is the (100-X)th percentile of dose distribution
    percentile = 100 - volume_percent
    return float(np.percentile(dose_values, percentile))

compute_dose_at_volume_cc ¶

compute_dose_at_volume_cc(dose: Dose, structure: Structure, volume_cc: float) -> float

Compute dose received by a given absolute volume in cc.

This computes D_Xcc (e.g., D_0.1cc = dose to hottest 0.1 cc).

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure`	`Structure`	Structure to analyze	required
`volume_cc`	`float`	Absolute volume in cubic centimeters	required

Returns:

Type	Description
`float`	Dose in Gy received by the specified volume

Examples:

>>> # D_0.1cc: dose to hottest 0.1 cc (common OAR metric)
>>> d_0_1cc = compute_dose_at_volume_cc(dose, brainstem, 0.1)
>>> print(f"D_0.1cc: {d_0_1cc:.2f} Gy")

Source code in src/dosemetrics/metrics/dvh.py

def compute_dose_at_volume_cc(
    dose: Dose, structure: Structure, volume_cc: float
) -> float:
    """
    Compute dose received by a given absolute volume in cc.

    This computes D_Xcc (e.g., D_0.1cc = dose to hottest 0.1 cc).

    Args:
        dose: Dose distribution object
        structure: Structure to analyze
        volume_cc: Absolute volume in cubic centimeters

    Returns:
        Dose in Gy received by the specified volume

    Examples:
        >>> # D_0.1cc: dose to hottest 0.1 cc (common OAR metric)
        >>> d_0_1cc = compute_dose_at_volume_cc(dose, brainstem, 0.1)
        >>> print(f"D_0.1cc: {d_0_1cc:.2f} Gy")
    """
    dose_values = dose.get_dose_in_structure(structure)

    if len(dose_values) == 0:
        return 0.0

    # Convert cc to number of voxels
    voxel_volume_cc = np.prod(structure.spacing) / 1000.0  # mm³ to cc
    num_voxels = int(np.round(volume_cc / voxel_volume_cc))

    if num_voxels >= len(dose_values):
        # Requested volume exceeds structure volume
        return float(np.min(dose_values))

    if num_voxels <= 0:
        return float(np.max(dose_values))

    # Sort dose values in descending order and take the dose at num_voxels
    sorted_doses = np.sort(dose_values)[::-1]
    return float(sorted_doses[num_voxels - 1])

compute_equivalent_uniform_dose ¶

compute_equivalent_uniform_dose(dose: Dose, structure: Structure, a_parameter: float) -> float

Compute Equivalent Uniform Dose (EUD).

EUD = (mean(D_i^a))^(1/a)

The a-parameter depends on tissue type: - a < 0 for tumors (emphasizes cold spots) - a > 0 for normal tissues (emphasizes hot spots)

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure`	`Structure`	Structure to analyze	required
`a_parameter`	`float`	Tissue-specific parameter	required

Returns:

Type	Description
`float`	Equivalent uniform dose in Gy

References

Niemierko, Med Phys 1997

Examples:

>>> # For tumor (emphasize underdosage)
>>> eud_tumor = compute_equivalent_uniform_dose(dose, ptv, a_parameter=-10)
>>>
>>> # For OAR (emphasize overdosage)
>>> eud_oar = compute_equivalent_uniform_dose(dose, brainstem, a_parameter=5)

Source code in src/dosemetrics/metrics/dvh.py

def compute_equivalent_uniform_dose(
    dose: Dose, structure: Structure, a_parameter: float
) -> float:
    """
    Compute Equivalent Uniform Dose (EUD).

    EUD = (mean(D_i^a))^(1/a)

    The a-parameter depends on tissue type:
    - a < 0 for tumors (emphasizes cold spots)
    - a > 0 for normal tissues (emphasizes hot spots)

    Args:
        dose: Dose distribution object
        structure: Structure to analyze
        a_parameter: Tissue-specific parameter

    Returns:
        Equivalent uniform dose in Gy

    References:
        Niemierko, Med Phys 1997

    Examples:
        >>> # For tumor (emphasize underdosage)
        >>> eud_tumor = compute_equivalent_uniform_dose(dose, ptv, a_parameter=-10)
        >>>
        >>> # For OAR (emphasize overdosage)
        >>> eud_oar = compute_equivalent_uniform_dose(dose, brainstem, a_parameter=5)
    """
    dose_values = dose.get_dose_in_structure(structure)

    if len(dose_values) == 0:
        return 0.0

    if a_parameter == 0:
        # Limit case: geometric mean
        return float(np.exp(np.mean(np.log(dose_values + 1e-10))))

    powered_doses = np.power(dose_values, a_parameter)
    mean_powered = np.mean(powered_doses)
    eud = np.power(mean_powered, 1.0 / a_parameter)

    return float(eud)

create_dvh_table ¶

create_dvh_table(dose: Dose, structure_set: StructureSet, structure_names: Optional[list] = None, max_dose: Optional[float] = None, step_size: float = 0.1) -> pd.DataFrame

Create DVH table for multiple structures in long format.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure_set`	`StructureSet`	StructureSet containing structures	required
`structure_names`	`Optional[list]`	List of structure names to include (optional)	`None`
`max_dose`	`Optional[float]`	Maximum dose for bins	`None`
`step_size`	`float`	Dose bin width in Gy	`0.1`

Returns:

Type	Description
`DataFrame`	DataFrame with columns [Dose, Structure, Volume]

Examples:

>>> dvh_df = create_dvh_table(dose, structures,
...                           structure_names=["PTV", "Brainstem", "SpinalCord"])
>>> dvh_df.to_csv("dvh_data.csv")

Source code in src/dosemetrics/metrics/dvh.py

def create_dvh_table(
    dose: Dose,
    structure_set: StructureSet,
    structure_names: Optional[list] = None,
    max_dose: Optional[float] = None,
    step_size: float = 0.1,
) -> pd.DataFrame:
    """
    Create DVH table for multiple structures in long format.

    Args:
        dose: Dose distribution object
        structure_set: StructureSet containing structures
        structure_names: List of structure names to include (optional)
        max_dose: Maximum dose for bins
        step_size: Dose bin width in Gy

    Returns:
        DataFrame with columns [Dose, Structure, Volume]

    Examples:
        >>> dvh_df = create_dvh_table(dose, structures,
        ...                           structure_names=["PTV", "Brainstem", "SpinalCord"])
        >>> dvh_df.to_csv("dvh_data.csv")
    """
    if structure_names is None:
        structure_names = structure_set.structure_names

    dvh_data = []

    for name in structure_names:
        try:
            structure = structure_set.get_structure(name)
            dose_bins, volumes = compute_dvh(dose, structure, max_dose, step_size)

            for dose_val, vol_val in zip(dose_bins, volumes):
                dvh_data.append(
                    {"Dose": dose_val, "Structure": name, "Volume": vol_val}
                )
        except ValueError:
            # Structure not found
            continue

    return pd.DataFrame(dvh_data)

extract_dvh_metrics ¶

extract_dvh_metrics(dose: Dose, structure: Structure, dose_thresholds: Optional[list] = None, volume_percentages: Optional[list] = None) -> Dict[str, float]

Extract common DVH metrics for a structure.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure`	`Structure`	Structure to analyze	required
`dose_thresholds`	`Optional[list]`	List of dose levels for VX metrics (Gy)	`None`
`volume_percentages`	`Optional[list]`	List of volume percentages for DX metrics	`None`

Returns:

Type	Description
`Dict[str, float]`	Dictionary with DVH metrics

Examples:

>>> metrics = extract_dvh_metrics(
...     dose, ptv,
...     dose_thresholds=[20, 40, 60],
...     volume_percentages=[2, 50, 95, 98]
... )
>>> print(metrics)
{'V20': 98.5, 'V40': 97.2, 'V60': 95.8, 'D2': 63.5, 'D50': 60.2, ...}

Source code in src/dosemetrics/metrics/dvh.py

def extract_dvh_metrics(
    dose: Dose,
    structure: Structure,
    dose_thresholds: Optional[list] = None,
    volume_percentages: Optional[list] = None,
) -> Dict[str, float]:
    """
    Extract common DVH metrics for a structure.

    Args:
        dose: Dose distribution object
        structure: Structure to analyze
        dose_thresholds: List of dose levels for VX metrics (Gy)
        volume_percentages: List of volume percentages for DX metrics

    Returns:
        Dictionary with DVH metrics

    Examples:
        >>> metrics = extract_dvh_metrics(
        ...     dose, ptv,
        ...     dose_thresholds=[20, 40, 60],
        ...     volume_percentages=[2, 50, 95, 98]
        ... )
        >>> print(metrics)
        {'V20': 98.5, 'V40': 97.2, 'V60': 95.8, 'D2': 63.5, 'D50': 60.2, ...}
    """
    metrics = {}

    # Volume at dose metrics (VX)
    if dose_thresholds:
        for threshold in dose_thresholds:
            v_x = compute_volume_at_dose(dose, structure, threshold)
            metrics[f"V{threshold}"] = v_x

    # Dose at volume metrics (DX)
    if volume_percentages:
        for vol_pct in volume_percentages:
            d_x = compute_dose_at_volume(dose, structure, vol_pct)
            metrics[f"D{vol_pct}"] = d_x

    return metrics

compute_dose_statistics ¶

compute_dose_statistics(dose: Dose, structure: Structure, verbose: bool = False) -> Dict[str, float]

Compute comprehensive dose statistics for a structure.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure`	`Structure`	Structure to analyze	required
`verbose`	`bool`	Print the returned statistics when `True` (default: `False`)	`False`

Returns:

Type	Description
`Dict[str, float]`	Dictionary with statistics including:
`Dict[str, float]`	mean_dose, max_dose, min_dose, median_dose, std_dose
`Dict[str, float]`	D95, D50, D05, D02, D98 (dose percentiles)

Examples:

>>> from dosemetrics.dose import Dose
>>> from dosemetrics.structure_set import StructureSet
>>> from dosemetrics.metrics import dvh
>>>
>>> dose = Dose.from_dicom("rtdose.dcm")
>>> structures = StructureSet(...)
>>> ptv = structures.get_structure("PTV")
>>>
>>> stats = dvh.compute_dose_statistics(dose, ptv, verbose=True)

Source code in src/dosemetrics/metrics/dvh.py

def compute_dose_statistics(
    dose: Dose, structure: Structure, verbose: bool = False
) -> Dict[str, float]:
    """
    Compute comprehensive dose statistics for a structure.

    Args:
        dose: Dose distribution object
        structure: Structure to analyze
        verbose: Print the returned statistics when ``True`` (default: ``False``)

    Returns:
        Dictionary with statistics including:
        - mean_dose, max_dose, min_dose, median_dose, std_dose
        - D95, D50, D05, D02, D98 (dose percentiles)

    Examples:
        >>> from dosemetrics.dose import Dose
        >>> from dosemetrics.structure_set import StructureSet
        >>> from dosemetrics.metrics import dvh
        >>>
        >>> dose = Dose.from_dicom("rtdose.dcm")
        >>> structures = StructureSet(...)
        >>> ptv = structures.get_structure("PTV")
        >>>
        >>> stats = dvh.compute_dose_statistics(dose, ptv, verbose=True)
    """
    dose_values = dose.get_dose_in_structure(structure)

    if len(dose_values) == 0:
        statistics = {
            "mean_dose": 0.0,
            "max_dose": 0.0,
            "min_dose": 0.0,
            "median_dose": 0.0,
            "std_dose": 0.0,
            "D95": 0.0,
            "D50": 0.0,
            "D05": 0.0,
            "D02": 0.0,
            "D98": 0.0,
        }
    else:
        statistics = {
            "mean_dose": float(np.mean(dose_values)),
            "max_dose": float(np.max(dose_values)),
            "min_dose": float(np.min(dose_values)),
            "median_dose": float(np.median(dose_values)),
            "std_dose": float(np.std(dose_values)),
            "D95": float(np.percentile(dose_values, 5)),  # 95% receives at least this
            "D50": float(np.percentile(dose_values, 50)),
            "D05": float(np.percentile(dose_values, 95)),  # 5% receives at least this
            "D02": float(np.percentile(dose_values, 98)),  # 2% receives at least this
            "D98": float(np.percentile(dose_values, 2)),  # 98% receives at least this
        }

    if verbose:
        print(f"{structure.name} dose statistics (Gy)")
        for label, key in (
            ("Mean", "mean_dose"),
            ("Minimum", "min_dose"),
            ("Maximum", "max_dose"),
            ("D98", "D98"),
            ("D95", "D95"),
            ("D50", "D50"),
            ("D05", "D05"),
            ("D02", "D02"),
        ):
            print(f"  {label:7s} {statistics[key]:.2f}")

    return statistics

compute_mean_dose ¶

compute_mean_dose(dose: Dose, structure: Structure) -> float

Compute mean dose in structure.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure`	`Structure`	Structure to analyze	required

Returns:

Type	Description
`float`	Mean dose in Gy

Source code in src/dosemetrics/metrics/dvh.py

def compute_mean_dose(dose: Dose, structure: Structure) -> float:
    """
    Compute mean dose in structure.

    Args:
        dose: Dose distribution object
        structure: Structure to analyze

    Returns:
        Mean dose in Gy
    """
    dose_values = dose.get_dose_in_structure(structure)
    return float(np.mean(dose_values)) if len(dose_values) > 0 else 0.0

compute_max_dose ¶

compute_max_dose(dose: Dose, structure: Structure) -> float

Compute maximum dose in structure.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure`	`Structure`	Structure to analyze	required

Returns:

Type	Description
`float`	Maximum dose in Gy

Source code in src/dosemetrics/metrics/dvh.py

def compute_max_dose(dose: Dose, structure: Structure) -> float:
    """
    Compute maximum dose in structure.

    Args:
        dose: Dose distribution object
        structure: Structure to analyze

    Returns:
        Maximum dose in Gy
    """
    dose_values = dose.get_dose_in_structure(structure)
    return float(np.max(dose_values)) if len(dose_values) > 0 else 0.0

compute_min_dose ¶

compute_min_dose(dose: Dose, structure: Structure) -> float

Compute minimum dose in structure.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure`	`Structure`	Structure to analyze	required

Returns:

Type	Description
`float`	Minimum dose in Gy

Source code in src/dosemetrics/metrics/dvh.py

def compute_min_dose(dose: Dose, structure: Structure) -> float:
    """
    Compute minimum dose in structure.

    Args:
        dose: Dose distribution object
        structure: Structure to analyze

    Returns:
        Minimum dose in Gy
    """
    dose_values = dose.get_dose_in_structure(structure)
    return float(np.min(dose_values)) if len(dose_values) > 0 else 0.0

compute_median_dose ¶

compute_median_dose(dose: Dose, structure: Structure) -> float

Compute median dose in structure.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure`	`Structure`	Structure to analyze	required

Returns:

Type	Description
`float`	Median dose in Gy

Source code in src/dosemetrics/metrics/dvh.py

def compute_median_dose(dose: Dose, structure: Structure) -> float:
    """
    Compute median dose in structure.

    Args:
        dose: Dose distribution object
        structure: Structure to analyze

    Returns:
        Median dose in Gy
    """
    dose_values = dose.get_dose_in_structure(structure)
    return float(np.median(dose_values)) if len(dose_values) > 0 else 0.0

compare_dvh_score ¶

compare_dvh_score(reference: Dose, evaluated: Dose, targets: Union[Structure, Sequence[Structure]], oars: Sequence[Structure] = ()) -> float

Compare plans using the complete OpenKBP DVH score.

The score is the unweighted mean absolute error over D1, D95, and D99 for every target and mean dose and D0.1cc for every OAR. It is reported in Gy.

Source code in src/dosemetrics/metrics/dvh.py

def compare_dvh_score(
    reference: Dose,
    evaluated: Dose,
    targets: Union[Structure, Sequence[Structure]],
    oars: Sequence[Structure] = (),
) -> float:
    """Compare plans using the complete OpenKBP DVH score.

    The score is the unweighted mean absolute error over D1, D95, and D99 for
    every target and mean dose and D0.1cc for every OAR. It is reported in Gy.
    """
    _validate_comparable_doses(reference, evaluated)
    target_structures = (targets,) if isinstance(targets, Structure) else tuple(targets)
    errors = []

    for target in target_structures:
        for volume_percent in (1.0, 95.0, 99.0):
            reference_value = compute_dose_at_volume(reference, target, volume_percent)
            evaluated_value = compute_dose_at_volume(evaluated, target, volume_percent)
            errors.append(abs(evaluated_value - reference_value))

    for oar in oars:
        errors.append(
            abs(compute_mean_dose(evaluated, oar) - compute_mean_dose(reference, oar))
        )
        errors.append(
            abs(
                compute_dose_at_volume_cc(evaluated, oar, 0.1)
                - compute_dose_at_volume_cc(reference, oar, 0.1)
            )
        )

    if not errors:
        raise ValueError("At least one target or OAR criterion is required")
    return float(np.mean(errors))

compute_dvh_auc ¶

compute_dvh_auc(dose: Dose, structure: Structure, num_bins: int = 100, normalize: bool = True, dose_range: Optional[Tuple[float, float]] = None) -> float

Compute the Area Under the DVH Curve (DVH-AUC) using the trapezoidal rule.

The DVH-AUC is the integral of volume percentage over the dose range. A higher AUC indicates that more of the structure volume receives higher doses. This is a single-distribution metric (unlike area-between-curves which compares two).

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure`	`Structure`	Structure to compute DVH AUC for	required
`num_bins`	`int`	Number of dose bins for DVH computation (default: 100)	`100`
`normalize`	`bool`	If True, normalize AUC to [0, 1] by dividing by the maximum possible area (100% volume × dose range). Default: True.	`True`
`dose_range`	`Optional[Tuple[float, float]]`	Fixed (min_dose, max_dose) in Gy for binning. Uses per-structure min/max if None.	`None`

Returns:

Type	Description
`float`	DVH AUC value. If normalize=True, returns value in [0, 1].
`float`	If normalize=False, returns the trapezoidal integral in Gy·% for a
`float`	non-degenerate dose range.

References

Adapted from DVHAUC metric in GDP-HMM AAPM Challenge.

Examples:

>>> auc = compute_dvh_auc(dose, ptv, normalize=True)
>>> print(f"DVH AUC (normalized): {auc:.3f}")
>>>
>>> # Compare two structures
>>> ptv_auc = compute_dvh_auc(dose, ptv)
>>> oar_auc = compute_dvh_auc(dose, brainstem)

Source code in src/dosemetrics/metrics/dvh.py

def compute_dvh_auc(
    dose: Dose,
    structure: Structure,
    num_bins: int = 100,
    normalize: bool = True,
    dose_range: Optional[Tuple[float, float]] = None,
) -> float:
    """
    Compute the Area Under the DVH Curve (DVH-AUC) using the trapezoidal rule.

    The DVH-AUC is the integral of volume percentage over the dose range.
    A higher AUC indicates that more of the structure volume receives higher doses.
    This is a single-distribution metric (unlike area-between-curves which compares two).

    Args:
        dose: Dose distribution object
        structure: Structure to compute DVH AUC for
        num_bins: Number of dose bins for DVH computation (default: 100)
        normalize: If True, normalize AUC to [0, 1] by dividing by the maximum
            possible area (100% volume × dose range). Default: True.
        dose_range: Fixed (min_dose, max_dose) in Gy for binning. Uses
            per-structure min/max if None.

    Returns:
        DVH AUC value. If normalize=True, returns value in [0, 1].
        If normalize=False, returns the trapezoidal integral in Gy·% for a
        non-degenerate dose range.

    References:
        Adapted from DVHAUC metric in GDP-HMM AAPM Challenge.

    Examples:
        >>> auc = compute_dvh_auc(dose, ptv, normalize=True)
        >>> print(f"DVH AUC (normalized): {auc:.3f}")
        >>>
        >>> # Compare two structures
        >>> ptv_auc = compute_dvh_auc(dose, ptv)
        >>> oar_auc = compute_dvh_auc(dose, brainstem)
    """
    dose_values = dose.get_dose_in_structure(structure)

    if len(dose_values) == 0:
        return 0.0

    if dose_range is not None:
        min_dose, max_dose = float(dose_range[0]), float(dose_range[1])
    else:
        min_dose = float(np.min(dose_values))
        max_dose = float(np.max(dose_values))

    if max_dose - min_dose < 1e-10:
        # All voxels at same dose: AUC = 1.0 normalized, or max_dose unnormalized
        return 1.0 if normalize else float(max_dose)

    dose_bins = np.linspace(min_dose, max_dose, num_bins)
    dvh_values = np.array(
        [100.0 * np.sum(dose_values >= d) / len(dose_values) for d in dose_bins]
    )

    auc = float(np.trapz(dvh_values, dose_bins))

    if normalize:
        max_area = (max_dose - min_dose) * 100.0
        auc = auc / max_area if max_area > 1e-10 else 0.0

    return auc

compute_dose_percentile ¶

compute_dose_percentile(dose: Dose, structure: Structure, percentile: float) -> float

Compute dose percentile (DX).

D95 means 95% of the volume receives at least this dose. This corresponds to the 5th percentile of the dose distribution.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`structure`	`Structure`	Structure to analyze	required
`percentile`	`float`	Volume percentage (0-100). For D95, use percentile=95	required

Returns:

Type	Description
`float`	Dose in Gy that the specified percentage of volume receives

Raises:

Type	Description
`ValueError`	If percentile is not in range 0-100

Examples:

>>> # D95: dose received by 95% of volume
>>> d95 = compute_dose_percentile(dose, ptv, 95)
>>>
>>> # D50: median dose
>>> d50 = compute_dose_percentile(dose, ptv, 50)
>>>
>>> # D05: near-maximum dose (hot spot)
>>> d05 = compute_dose_percentile(dose, ptv, 5)

Source code in src/dosemetrics/metrics/dvh.py

def compute_dose_percentile(
    dose: Dose, structure: Structure, percentile: float
) -> float:
    """
    Compute dose percentile (DX).

    D95 means 95% of the volume receives at least this dose.
    This corresponds to the 5th percentile of the dose distribution.

    Args:
        dose: Dose distribution object
        structure: Structure to analyze
        percentile: Volume percentage (0-100). For D95, use percentile=95

    Returns:
        Dose in Gy that the specified percentage of volume receives

    Raises:
        ValueError: If percentile is not in range 0-100

    Examples:
        >>> # D95: dose received by 95% of volume
        >>> d95 = compute_dose_percentile(dose, ptv, 95)
        >>>
        >>> # D50: median dose
        >>> d50 = compute_dose_percentile(dose, ptv, 50)
        >>>
        >>> # D05: near-maximum dose (hot spot)
        >>> d05 = compute_dose_percentile(dose, ptv, 5)
    """
    if not 0 <= percentile <= 100:
        raise ValueError(f"Percentile must be 0-100, got {percentile}")

    dose_values = dose.get_dose_in_structure(structure)

    if len(dose_values) == 0:
        return 0.0

    # DX means X% receives AT LEAST this dose
    # This is the (100-X)th percentile of the dose array
    return float(np.percentile(dose_values, 100 - percentile))

compare_dvh_wasserstein ¶

compare_dvh_wasserstein(reference: Dose, evaluated: Dose, structure: Structure) -> float

Return the Wasserstein distance between structure dose samples, in Gy.

Source code in src/dosemetrics/metrics/dvh.py

def compare_dvh_wasserstein(
    reference: Dose,
    evaluated: Dose,
    structure: Structure,
) -> float:
    """Return the Wasserstein distance between structure dose samples, in Gy."""
    from scipy.stats import wasserstein_distance

    _validate_comparable_doses(reference, evaluated)
    reference_values = reference.get_dose_in_structure(structure)
    evaluated_values = evaluated.get_dose_in_structure(structure)
    if len(reference_values) == 0:
        return float("nan")
    return float(wasserstein_distance(reference_values, evaluated_values))

compare_dvh_area ¶

compare_dvh_area(reference: Dose, evaluated: Dose, structure: Structure, norm: str = 'l1', step_size: float = 0.1) -> float

Compare two cumulative DVHs using an integrated L1 or L2 distance.

norm="l1" integrates the pointwise absolute separation. norm="l2" returns the square root of the integrated squared separation.

Source code in src/dosemetrics/metrics/dvh.py

def compare_dvh_area(
    reference: Dose,
    evaluated: Dose,
    structure: Structure,
    norm: str = "l1",
    step_size: float = 0.1,
) -> float:
    """Compare two cumulative DVHs using an integrated L1 or L2 distance.

    ``norm="l1"`` integrates the pointwise absolute separation. ``norm="l2"``
    returns the square root of the integrated squared separation.
    """
    norm = norm.lower()
    if norm not in {"l1", "l2"}:
        raise ValueError("norm must be 'l1' or 'l2'")
    bins, reference_volume, evaluated_volume = _common_dvh_grid(
        reference, evaluated, structure, step_size
    )
    difference = evaluated_volume - reference_volume
    integrand = np.abs(difference) if norm == "l1" else difference**2
    integral = float(np.trapz(integrand, bins))
    return integral if norm == "l1" else float(np.sqrt(integral))

compare_dvh_chi_square ¶

compare_dvh_chi_square(reference: Dose, evaluated: Dose, structure: Structure, step_size: float = 0.1) -> Tuple[float, float]

Compare differential DVHs with a chi-square goodness-of-fit test.

Source code in src/dosemetrics/metrics/dvh.py

def compare_dvh_chi_square(
    reference: Dose,
    evaluated: Dose,
    structure: Structure,
    step_size: float = 0.1,
) -> Tuple[float, float]:
    """Compare differential DVHs with a chi-square goodness-of-fit test."""
    from scipy.stats import chisquare

    _, reference_volume, evaluated_volume = _common_dvh_grid(
        reference, evaluated, structure, step_size
    )
    expected = np.maximum(-np.diff(np.append(reference_volume, 0.0)), 0.0)
    observed = np.maximum(-np.diff(np.append(evaluated_volume, 0.0)), 0.0)
    expected += np.finfo(float).eps
    observed *= np.sum(expected) / max(np.sum(observed), np.finfo(float).eps)
    statistic, p_value = chisquare(observed, expected)
    return float(statistic), float(p_value)

compare_dvh_ks ¶

compare_dvh_ks(reference: Dose, evaluated: Dose, structure: Structure) -> Tuple[float, float]

Compare structure dose samples with a two-sample KS test.

Source code in src/dosemetrics/metrics/dvh.py

def compare_dvh_ks(
    reference: Dose,
    evaluated: Dose,
    structure: Structure,
) -> Tuple[float, float]:
    """Compare structure dose samples with a two-sample KS test."""
    from scipy.stats import ks_2samp

    _validate_comparable_doses(reference, evaluated)
    reference_values = reference.get_dose_in_structure(structure)
    evaluated_values = evaluated.get_dose_in_structure(structure)
    if len(reference_values) == 0:
        return float("nan"), float("nan")
    statistic, p_value = ks_2samp(reference_values, evaluated_values)
    return float(statistic), float(p_value)

compute_dvh_confidence_interval ¶

compute_dvh_confidence_interval(doses: Sequence[Dose], structure: Structure, confidence: float = 0.95, step_size: float = 0.1) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]

Summarize a collection of plans with mean DVH and percentile interval.

Source code in src/dosemetrics/metrics/dvh.py

def compute_dvh_confidence_interval(
    doses: Sequence[Dose],
    structure: Structure,
    confidence: float = 0.95,
    step_size: float = 0.1,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """Summarize a collection of plans with mean DVH and percentile interval."""
    if not 0 < confidence < 1:
        raise ValueError("confidence must be between 0 and 1")
    bins, dvhs = _aligned_dvhs(doses, structure, step_size)
    tail = (1.0 - confidence) * 50.0
    return (
        bins,
        np.mean(dvhs, axis=0),
        np.percentile(dvhs, tail, axis=0),
        np.percentile(dvhs, 100.0 - tail, axis=0),
    )

compute_dvh_bandwidth ¶

compute_dvh_bandwidth(doses: Sequence[Dose], structure: Structure, step_size: float = 0.1) -> Tuple[np.ndarray, np.ndarray]

Return dose bins and the pointwise max-minus-min DVH bandwidth.

Source code in src/dosemetrics/metrics/dvh.py

def compute_dvh_bandwidth(
    doses: Sequence[Dose],
    structure: Structure,
    step_size: float = 0.1,
) -> Tuple[np.ndarray, np.ndarray]:
    """Return dose bins and the pointwise max-minus-min DVH bandwidth."""
    bins, dvhs = _aligned_dvhs(doses, structure, step_size)
    return bins, np.max(dvhs, axis=0) - np.min(dvhs, axis=0)

compare_dvh_similarity ¶

compare_dvh_similarity(reference: Dose, evaluated: Dose, structure: Structure, method: str = 'dice', step_size: float = 0.1) -> float

Compare two DVHs using Dice, Jaccard, correlation, or cosine similarity.

Source code in src/dosemetrics/metrics/dvh.py

def compare_dvh_similarity(
    reference: Dose,
    evaluated: Dose,
    structure: Structure,
    method: str = "dice",
    step_size: float = 0.1,
) -> float:
    """Compare two DVHs using Dice, Jaccard, correlation, or cosine similarity."""
    method = method.lower()
    if method not in {"dice", "jaccard", "correlation", "cosine"}:
        raise ValueError("method must be 'dice', 'jaccard', 'correlation', or 'cosine'")
    _, reference_volume, evaluated_volume = _common_dvh_grid(
        reference, evaluated, structure, step_size
    )
    intersection = np.minimum(reference_volume, evaluated_volume)
    if method == "dice":
        denominator = np.sum(reference_volume + evaluated_volume)
        return float(2.0 * np.sum(intersection) / denominator) if denominator else 0.0
    if method == "jaccard":
        union = np.sum(np.maximum(reference_volume, evaluated_volume))
        return float(np.sum(intersection) / union) if union else 0.0
    if method == "correlation":
        correlation = np.corrcoef(reference_volume, evaluated_volume)[0, 1]
        return float(correlation) if np.isfinite(correlation) else 0.0
    denominator = np.linalg.norm(reference_volume) * np.linalg.norm(evaluated_volume)
    return (
        float(np.dot(reference_volume, evaluated_volume) / denominator)
        if denominator
        else 0.0
    )

compare_oar_dvh_auc ¶

compare_oar_dvh_auc(reference: Dose, evaluated: Dose, oar: Structure, num_bins: int = 100) -> float

Return |AUC_evaluated - AUC_reference| for one OAR, in Gy.

Source code in src/dosemetrics/metrics/dvh.py

def compare_oar_dvh_auc(
    reference: Dose,
    evaluated: Dose,
    oar: Structure,
    num_bins: int = 100,
) -> float:
    """Return ``|AUC_evaluated - AUC_reference|`` for one OAR, in Gy."""
    _validate_comparable_doses(reference, evaluated)
    if num_bins < 2:
        raise ValueError("num_bins must be at least 2")
    reference_values = reference.get_dose_in_structure(oar)
    evaluated_values = evaluated.get_dose_in_structure(oar)
    if len(reference_values) == 0:
        return float("nan")
    min_dose = float(min(np.min(reference_values), np.min(evaluated_values)))
    max_dose = float(max(np.max(reference_values), np.max(evaluated_values)))
    if np.isclose(min_dose, max_dose):
        return 0.0
    bins = np.linspace(min_dose, max_dose, num_bins)
    reference_dvh = np.asarray([np.mean(reference_values >= d) for d in bins])
    evaluated_dvh = np.asarray([np.mean(evaluated_values >= d) for d in bins])
    return float(abs(np.trapz(evaluated_dvh, bins) - np.trapz(reference_dvh, bins)))

compare_mean_oar_dvh_auc ¶

compare_mean_oar_dvh_auc(reference: Dose, evaluated: Dose, oars: Sequence[Structure], num_bins: int = 100) -> float

Average the per-OAR absolute AUC differences.

Source code in src/dosemetrics/metrics/dvh.py

def compare_mean_oar_dvh_auc(
    reference: Dose,
    evaluated: Dose,
    oars: Sequence[Structure],
    num_bins: int = 100,
) -> float:
    """Average the per-OAR absolute AUC differences."""
    if not oars:
        raise ValueError("At least one OAR is required")
    return float(
        np.mean(
            [
                compare_oar_dvh_auc(reference, evaluated, oar, num_bins=num_bins)
                for oar in oars
            ]
        )
    )

`gamma`¶

gamma ¶

Gamma analysis for dose distribution comparison.

This module provides gamma index calculation following the methodology of Low et al. (1998) and subsequent refinements.

References

Low DA, Harms WB, Mutic S, Purdy JA. "A technique for the quantitative evaluation of dose distributions." Med Phys. 1998;25(5):656-61.
Depuydt T, Van Esch A, Huyskens DP. "A quantitative evaluation of IMRT dose distributions: refinement and clinical assessment of the gamma evaluation." Radiother Oncol. 2002;62(3):309-19.

Classes¶

Functions:¶

compare_gamma_index ¶

compare_gamma_index(reference: Dose, evaluated: Dose, dose_criterion_percent: float = 3.0, distance_criterion_mm: float = 3.0, dose_threshold_percent: float = 10.0, global_normalization: bool = True, max_search_distance_mm: Optional[float] = None) -> np.ndarray

Compute 3D gamma index between reference and evaluated dose distributions.

The gamma index quantifies the agreement between two dose distributions by combining dose difference and distance-to-agreement criteria.

Parameters¶

reference : Dose Reference (planned) dose distribution. evaluated : Dose Evaluated (measured/calculated) dose distribution to compare. dose_criterion_percent : float, optional Dose difference criterion as percentage (default: 3.0 for 3%). distance_criterion_mm : float, optional Distance-to-agreement criterion in mm (default: 3.0 for 3mm). dose_threshold_percent : float, optional Low dose threshold below which gamma is not calculated (default: 10%). global_normalization : bool, optional If True, normalize to global maximum dose. If False, use local dose (default: True). max_search_distance_mm : float, optional Maximum search distance for gamma calculation. If None, uses 3 * distance_criterion_mm (default: None).

Returns¶

gamma : np.ndarray 3D array of gamma values. Values < 1 indicate passing points, values >= 1 indicate failing points. NaN for points below threshold.

Notes¶

Common gamma criteria: - Clinical QA: 3%/3mm (dose_criterion=3.0, distance_criterion=3.0) - Stricter QA: 2%/2mm - Research: 1%/1mm

The gamma passing rate is typically calculated as the percentage of points with gamma <= 1.0.

Examples¶

gamma = compare_gamma_index(planned_dose, measured_dose) passing_rate = np.sum(gamma <= 1.0) / np.sum(~np.isnan(gamma)) * 100 print(f"Gamma passing rate: {passing_rate:.1f}%")

Raises¶

ValueError If dose distributions have incompatible geometry.

Source code in src/dosemetrics/metrics/gamma.py

def compare_gamma_index(
    reference: Dose,
    evaluated: Dose,
    dose_criterion_percent: float = 3.0,
    distance_criterion_mm: float = 3.0,
    dose_threshold_percent: float = 10.0,
    global_normalization: bool = True,
    max_search_distance_mm: Optional[float] = None,
) -> np.ndarray:
    """
    Compute 3D gamma index between reference and evaluated dose distributions.

    The gamma index quantifies the agreement between two dose distributions by
    combining dose difference and distance-to-agreement criteria.

    Parameters
    ----------
    reference : Dose
        Reference (planned) dose distribution.
    evaluated : Dose
        Evaluated (measured/calculated) dose distribution to compare.
    dose_criterion_percent : float, optional
        Dose difference criterion as percentage (default: 3.0 for 3%).
    distance_criterion_mm : float, optional
        Distance-to-agreement criterion in mm (default: 3.0 for 3mm).
    dose_threshold_percent : float, optional
        Low dose threshold below which gamma is not calculated (default: 10%).
    global_normalization : bool, optional
        If True, normalize to global maximum dose. If False, use local dose
        (default: True).
    max_search_distance_mm : float, optional
        Maximum search distance for gamma calculation. If None, uses
        3 * distance_criterion_mm (default: None).

    Returns
    -------
    gamma : np.ndarray
        3D array of gamma values. Values < 1 indicate passing points,
        values >= 1 indicate failing points. NaN for points below threshold.

    Notes
    -----
    Common gamma criteria:
        - Clinical QA: 3%/3mm (dose_criterion=3.0, distance_criterion=3.0)
        - Stricter QA: 2%/2mm
        - Research: 1%/1mm

    The gamma passing rate is typically calculated as the percentage of
    points with gamma <= 1.0.

    Examples
    --------
    >>> gamma = compare_gamma_index(planned_dose, measured_dose)
    >>> passing_rate = np.sum(gamma <= 1.0) / np.sum(~np.isnan(gamma)) * 100
    >>> print(f"Gamma passing rate: {passing_rate:.1f}%")

    Raises
    ------
    ValueError
        If dose distributions have incompatible geometry.
    """
    # Validate spatial compatibility
    if reference.dose_array.shape != evaluated.dose_array.shape:
        raise ValueError(
            f"Dose shapes must match: {reference.dose_array.shape} vs "
            f"{evaluated.dose_array.shape}"
        )

    if not np.allclose(reference.spacing, evaluated.spacing):
        raise ValueError(
            f"Dose spacings must match: {reference.spacing} vs {evaluated.spacing}"
        )

    # Get dose arrays and spatial information
    ref_dose = np.ascontiguousarray(reference.dose_array)
    eval_dose = np.ascontiguousarray(evaluated.dose_array)
    spacing = np.array(reference.spacing)

    # Set max search distance
    if max_search_distance_mm is None:
        max_search_distance_mm = 3 * distance_criterion_mm

    # Determine normalization dose
    if global_normalization:
        normalization_dose = np.max(ref_dose)
    else:
        normalization_dose = 0.0  # The kernel uses each reference voxel locally.

    # Calculate absolute dose threshold
    dose_threshold = dose_threshold_percent / 100.0 * np.max(ref_dose)

    # Create search grid offsets (in voxel indices)
    search_radius_voxels = np.ceil(max_search_distance_mm / spacing).astype(int)

    # For efficiency, create a search template
    i_range = np.arange(-search_radius_voxels[0], search_radius_voxels[0] + 1)
    j_range = np.arange(-search_radius_voxels[1], search_radius_voxels[1] + 1)
    k_range = np.arange(-search_radius_voxels[2], search_radius_voxels[2] + 1)

    di, dj, dk = np.meshgrid(i_range, j_range, k_range, indexing="ij")

    # Physical distances for search template
    dx = di * spacing[0]
    dy = dj * spacing[1]
    dz = dk * spacing[2]
    distances_template = np.sqrt(dx**2 + dy**2 + dz**2)

    # Only keep points within search distance
    valid_search = distances_template <= max_search_distance_mm
    di_valid = di[valid_search]
    dj_valid = dj[valid_search]
    dk_valid = dk[valid_search]
    distances_valid = distances_template[valid_search]

    # Sorting enables an exact branch-and-bound search: once an offset's
    # spatial contribution exceeds the best gamma found so far, no later
    # candidate can improve it because the dose contribution is non-negative.
    spatial_gamma_squared = (distances_valid / distance_criterion_mm) ** 2
    order = np.argsort(spatial_gamma_squared, kind="stable")

    return _compute_gamma_grid(
        ref_dose,
        eval_dose,
        np.ascontiguousarray(di_valid[order]),
        np.ascontiguousarray(dj_valid[order]),
        np.ascontiguousarray(dk_valid[order]),
        np.ascontiguousarray(spatial_gamma_squared[order]),
        dose_threshold,
        float(normalization_dose),
        dose_criterion_percent / 100.0,
        global_normalization,
    )

compute_gamma_passing_rate ¶

compute_gamma_passing_rate(gamma: ndarray, threshold: float = 1.0) -> float

Compute gamma passing rate from gamma index array.

Parameters¶

gamma : np.ndarray Gamma index values from compare_gamma_index(). threshold : float, optional Gamma threshold for passing (default: 1.0).

Returns¶

passing_rate : float Percentage of points with gamma <= threshold (0-100).

Source code in src/dosemetrics/metrics/gamma.py

def compute_gamma_passing_rate(gamma: np.ndarray, threshold: float = 1.0) -> float:
    """
    Compute gamma passing rate from gamma index array.

    Parameters
    ----------
    gamma : np.ndarray
        Gamma index values from compare_gamma_index().
    threshold : float, optional
        Gamma threshold for passing (default: 1.0).

    Returns
    -------
    passing_rate : float
        Percentage of points with gamma <= threshold (0-100).
    """
    # Remove NaN values (below threshold points)
    valid_gamma = gamma[~np.isnan(gamma)]

    if len(valid_gamma) == 0:
        return 0.0

    # Calculate passing rate
    passing = np.sum(valid_gamma <= threshold)
    total = len(valid_gamma)
    passing_rate = (passing / total) * 100.0

    return float(passing_rate)

compute_gamma_statistics ¶

compute_gamma_statistics(gamma: ndarray) -> Dict[str, float]

Compute comprehensive statistics from gamma index array.

Parameters¶

gamma : np.ndarray Gamma index values.

Returns¶

stats : dict Dictionary containing: - 'passing_rate_1_0': Passing rate at gamma=1.0 - 'mean_gamma': Mean gamma value - 'max_gamma': Maximum gamma value - 'gamma_50': Median gamma value - 'gamma_95': 95th percentile gamma

Source code in src/dosemetrics/metrics/gamma.py

def compute_gamma_statistics(gamma: np.ndarray) -> Dict[str, float]:
    """
    Compute comprehensive statistics from gamma index array.

    Parameters
    ----------
    gamma : np.ndarray
        Gamma index values.

    Returns
    -------
    stats : dict
        Dictionary containing:
            - 'passing_rate_1_0': Passing rate at gamma=1.0
            - 'mean_gamma': Mean gamma value
            - 'max_gamma': Maximum gamma value
            - 'gamma_50': Median gamma value
            - 'gamma_95': 95th percentile gamma
    """
    # Remove NaN values
    valid_gamma = gamma[~np.isnan(gamma)]

    if len(valid_gamma) == 0:
        return {
            "passing_rate_1_0": 0.0,
            "mean_gamma": np.nan,
            "max_gamma": np.nan,
            "gamma_50": np.nan,
            "gamma_95": np.nan,
        }

    stats = {
        "passing_rate_1_0": compute_gamma_passing_rate(gamma, threshold=1.0),
        "mean_gamma": float(np.mean(valid_gamma)),
        "max_gamma": float(np.max(valid_gamma)),
        "gamma_50": float(np.percentile(valid_gamma, 50)),
        "gamma_95": float(np.percentile(valid_gamma, 95)),
    }

    return stats

compare_2d_gamma ¶

compare_2d_gamma(reference: ndarray, evaluated: ndarray, dose_criterion_percent: float = 3.0, distance_criterion_mm: float = 3.0, pixel_spacing: Tuple[float, float] = (1.0, 1.0)) -> np.ndarray

Compute 2D gamma index for a single slice (faster than 3D).

Parameters¶

reference : np.ndarray 2D reference dose slice. evaluated : np.ndarray 2D evaluated dose slice. dose_criterion_percent : float Dose criterion (%). distance_criterion_mm : float Distance criterion (mm). pixel_spacing : tuple of float Pixel spacing in mm (row_spacing, col_spacing).

Returns¶

gamma : np.ndarray 2D gamma index array.

Raises¶

ValueError If slice shapes don't match or are not 2D.

Source code in src/dosemetrics/metrics/gamma.py

def compare_2d_gamma(
    reference: np.ndarray,
    evaluated: np.ndarray,
    dose_criterion_percent: float = 3.0,
    distance_criterion_mm: float = 3.0,
    pixel_spacing: Tuple[float, float] = (1.0, 1.0),
) -> np.ndarray:
    """
    Compute 2D gamma index for a single slice (faster than 3D).

    Parameters
    ----------
    reference : np.ndarray
        2D reference dose slice.
    evaluated : np.ndarray
        2D evaluated dose slice.
    dose_criterion_percent : float
        Dose criterion (%).
    distance_criterion_mm : float
        Distance criterion (mm).
    pixel_spacing : tuple of float
        Pixel spacing in mm (row_spacing, col_spacing).

    Returns
    -------
    gamma : np.ndarray
        2D gamma index array.

    Raises
    ------
    ValueError
        If slice shapes don't match or are not 2D.
    """
    # Validate input
    if reference.ndim != 2:
        raise ValueError(f"Reference slice must be 2D, got shape {reference.shape}")
    if evaluated.ndim != 2:
        raise ValueError(f"Evaluated slice must be 2D, got shape {evaluated.shape}")
    if reference.shape != evaluated.shape:
        raise ValueError(
            f"Slice shapes must match: {reference.shape} vs {evaluated.shape}"
        )

    # Get shape and spacing
    shape = reference.shape
    spacing = np.array(pixel_spacing)

    # Determine normalization dose (global max)
    normalization_dose = np.max(reference)
    if normalization_dose == 0:
        normalization_dose = 1.0

    # Initialize gamma array
    gamma_result = np.full(shape, np.nan, dtype=np.float32)

    # Create search grid offsets (in voxel indices)
    max_search_distance_mm = 3 * distance_criterion_mm
    search_radius_voxels = np.ceil(max_search_distance_mm / spacing).astype(int)

    # Create search template
    i_range = np.arange(-search_radius_voxels[0], search_radius_voxels[0] + 1)
    j_range = np.arange(-search_radius_voxels[1], search_radius_voxels[1] + 1)

    di, dj = np.meshgrid(i_range, j_range, indexing="ij")

    # Physical distances for search template
    dx = di * spacing[0]
    dy = dj * spacing[1]
    distances_template = np.sqrt(dx**2 + dy**2)

    # Only keep points within search distance
    valid_search = distances_template <= max_search_distance_mm
    di_valid = di[valid_search]
    dj_valid = dj[valid_search]
    distances_valid = distances_template[valid_search]

    # Iterate through reference dose grid
    for i in range(shape[0]):
        for j in range(shape[1]):
            ref_value = reference[i, j]

            # Get search positions in index space
            i_search = i + di_valid
            j_search = j + dj_valid

            # Filter for valid indices
            valid_mask = (
                (i_search >= 0)
                & (i_search < shape[0])
                & (j_search >= 0)
                & (j_search < shape[1])
            )

            i_search = i_search[valid_mask]
            j_search = j_search[valid_mask]
            local_distances = distances_valid[valid_mask]

            if len(i_search) == 0:
                continue

            # Get evaluated dose values at search positions
            eval_values = evaluated[i_search, j_search]

            # Compute dose differences (normalized)
            dose_diff = np.abs(eval_values - ref_value) / normalization_dose

            # Compute gamma values
            gamma_values = np.sqrt(
                (local_distances / distance_criterion_mm) ** 2
                + (dose_diff / (dose_criterion_percent / 100.0)) ** 2
            )

            # Store minimum gamma value
            gamma_result[i, j] = np.min(gamma_values)

    return gamma_result

compare_gamma_index_gpu ¶

compare_gamma_index_gpu(reference: Dose, evaluated: Dose, dose_criterion_percent: float = 3.0, distance_criterion_mm: float = 3.0) -> np.ndarray

GPU-accelerated gamma index calculation (requires CuPy or similar).

Note: This is a placeholder for future GPU acceleration using CuPy or similar.

Parameters¶

reference : Dose Reference dose. evaluated : Dose Evaluated dose. dose_criterion_percent : float Dose criterion (%). distance_criterion_mm : float Distance criterion (mm).

Returns¶

gamma : np.ndarray Gamma index array.

Raises¶

NotImplementedError GPU acceleration not implemented yet.

Source code in src/dosemetrics/metrics/gamma.py

def compare_gamma_index_gpu(
    reference: Dose,
    evaluated: Dose,
    dose_criterion_percent: float = 3.0,
    distance_criterion_mm: float = 3.0,
) -> np.ndarray:
    """
    GPU-accelerated gamma index calculation (requires CuPy or similar).

    Note: This is a placeholder for future GPU acceleration using CuPy or similar.

    Parameters
    ----------
    reference : Dose
        Reference dose.
    evaluated : Dose
        Evaluated dose.
    dose_criterion_percent : float
        Dose criterion (%).
    distance_criterion_mm : float
        Distance criterion (mm).

    Returns
    -------
    gamma : np.ndarray
        Gamma index array.

    Raises
    ------
    NotImplementedError
        GPU acceleration not implemented yet.
    """
    warnings.warn(
        "GPU-accelerated gamma is not implemented. "
        "Use compare_gamma_index() for CPU-based calculation.",
        FutureWarning,
    )
    raise NotImplementedError(
        "GPU-accelerated gamma not implemented. Use compare_gamma_index()."
    )

`geometric`¶

geometric ¶

Geometric similarity and overlap metrics for structure comparison.

This module provides metrics to compare two structure sets, typically used for evaluating auto-segmentation algorithms or inter-observer variability.

Classes¶

Functions:¶

compute_dice_coefficient ¶

compute_dice_coefficient(structure1: Structure, structure2: Structure) -> float

Compute Dice coefficient (Sørensen-Dice index).

Dice = 2 * |A ∩ B| / (|A| + |B|)

Measures overlap between two structures. Range [0, 1], where 1 is perfect overlap.

Parameters:

Name	Type	Description	Default
`structure1`	`Structure`	First structure	required
`structure2`	`Structure`	Second structure	required

Returns:

Type	Description
`float`	Dice coefficient (0-1)

References

Dice, Ecology 1945; Sørensen, Biologiske Skrifter 1948

Examples:

>>> auto_ptv = structures_auto.get_structure("PTV")
>>> manual_ptv = structures_manual.get_structure("PTV")
>>> dice = compute_dice_coefficient(auto_ptv, manual_ptv)
>>> print(f"Dice: {dice:.3f}")

Source code in src/dosemetrics/metrics/geometric.py

def compute_dice_coefficient(structure1: Structure, structure2: Structure) -> float:
    """
    Compute Dice coefficient (Sørensen-Dice index).

    Dice = 2 * |A ∩ B| / (|A| + |B|)

    Measures overlap between two structures. Range [0, 1], where 1 is perfect overlap.

    Args:
        structure1: First structure
        structure2: Second structure

    Returns:
        Dice coefficient (0-1)

    References:
        Dice, Ecology 1945; Sørensen, Biologiske Skrifter 1948

    Examples:
        >>> auto_ptv = structures_auto.get_structure("PTV")
        >>> manual_ptv = structures_manual.get_structure("PTV")
        >>> dice = compute_dice_coefficient(auto_ptv, manual_ptv)
        >>> print(f"Dice: {dice:.3f}")
    """
    if structure1.mask is None or structure2.mask is None:
        return 0.0

    intersection = np.logical_and(structure1.mask, structure2.mask)
    sum_volumes = structure1.volume_voxels() + structure2.volume_voxels()

    if sum_volumes == 0:
        return 0.0

    return float(2.0 * np.sum(intersection) / sum_volumes)

compute_jaccard_index ¶

compute_jaccard_index(structure1: Structure, structure2: Structure) -> float

Compute Jaccard index (Intersection over Union, IoU).

Jaccard = |A ∩ B| / |A ∪ B|

Measures overlap between two structures. Range [0, 1], where 1 is perfect overlap. More conservative than Dice coefficient.

Parameters:

Name	Type	Description	Default
`structure1`	`Structure`	First structure	required
`structure2`	`Structure`	Second structure	required

Returns:

Type	Description
`float`	Jaccard index (0-1)

References

Jaccard, New Phytologist 1912

Examples:

>>> jaccard = compute_jaccard_index(auto_ptv, manual_ptv)
>>> print(f"IoU: {jaccard:.3f}")

Source code in src/dosemetrics/metrics/geometric.py

def compute_jaccard_index(structure1: Structure, structure2: Structure) -> float:
    """
    Compute Jaccard index (Intersection over Union, IoU).

    Jaccard = |A ∩ B| / |A ∪ B|

    Measures overlap between two structures. Range [0, 1], where 1 is perfect overlap.
    More conservative than Dice coefficient.

    Args:
        structure1: First structure
        structure2: Second structure

    Returns:
        Jaccard index (0-1)

    References:
        Jaccard, New Phytologist 1912

    Examples:
        >>> jaccard = compute_jaccard_index(auto_ptv, manual_ptv)
        >>> print(f"IoU: {jaccard:.3f}")
    """
    if structure1.mask is None or structure2.mask is None:
        return 0.0

    intersection = np.logical_and(structure1.mask, structure2.mask)
    union = np.logical_or(structure1.mask, structure2.mask)

    union_sum = np.sum(union)
    if union_sum == 0:
        return 0.0

    return float(np.sum(intersection) / union_sum)

compute_volume_difference ¶

compute_volume_difference(structure1: Structure, structure2: Structure) -> float

Compute absolute volume difference.

Parameters:

Name	Type	Description	Default
`structure1`	`Structure`	First structure	required
`structure2`	`Structure`	Second structure	required

Returns:

Type	Description
`float`	Absolute volume difference in cubic centimeters

Examples:

>>> vol_diff = compute_volume_difference(auto_ptv, manual_ptv)
>>> print(f"Volume difference: {vol_diff:.2f} cc")

Source code in src/dosemetrics/metrics/geometric.py

def compute_volume_difference(structure1: Structure, structure2: Structure) -> float:
    """
    Compute absolute volume difference.

    Args:
        structure1: First structure
        structure2: Second structure

    Returns:
        Absolute volume difference in cubic centimeters

    Examples:
        >>> vol_diff = compute_volume_difference(auto_ptv, manual_ptv)
        >>> print(f"Volume difference: {vol_diff:.2f} cc")
    """
    return abs(structure1.volume_cc() - structure2.volume_cc())

compute_volume_ratio ¶

compute_volume_ratio(structure1: Structure, structure2: Structure) -> float

Compute volume ratio V1/V2.

Parameters:

Name	Type	Description	Default
`structure1`	`Structure`	First structure (numerator)	required
`structure2`	`Structure`	Second structure (denominator)	required

Returns:

Type	Description
`float`	Volume ratio (dimensionless)

Examples:

>>> ratio = compute_volume_ratio(auto_ptv, manual_ptv)
>>> print(f"Volume ratio: {ratio:.3f}")

Source code in src/dosemetrics/metrics/geometric.py

def compute_volume_ratio(structure1: Structure, structure2: Structure) -> float:
    """
    Compute volume ratio V1/V2.

    Args:
        structure1: First structure (numerator)
        structure2: Second structure (denominator)

    Returns:
        Volume ratio (dimensionless)

    Examples:
        >>> ratio = compute_volume_ratio(auto_ptv, manual_ptv)
        >>> print(f"Volume ratio: {ratio:.3f}")
    """
    v2 = structure2.volume_cc()
    if v2 == 0:
        return float('inf') if structure1.volume_cc() > 0 else 1.0

    return structure1.volume_cc() / v2

compute_sensitivity ¶

compute_sensitivity(structure1: Structure, structure2: Structure) -> float

Compute sensitivity (recall, true positive rate).

Sensitivity = TP / (TP + FN) = |A ∩ B| / |B|

Measures how much of structure2 is covered by structure1.

Parameters:

Name	Type	Description	Default
`structure1`	`Structure`	Predicted/test structure	required
`structure2`	`Structure`	Reference/ground truth structure	required

Returns:

Type	Description
`float`	Sensitivity (0-1)

Examples:

>>> sens = compute_sensitivity(auto_structure, manual_structure)
>>> print(f"Sensitivity: {sens:.3f}")

Source code in src/dosemetrics/metrics/geometric.py

def compute_sensitivity(structure1: Structure, structure2: Structure) -> float:
    """
    Compute sensitivity (recall, true positive rate).

    Sensitivity = TP / (TP + FN) = |A ∩ B| / |B|

    Measures how much of structure2 is covered by structure1.

    Args:
        structure1: Predicted/test structure
        structure2: Reference/ground truth structure

    Returns:
        Sensitivity (0-1)

    Examples:
        >>> sens = compute_sensitivity(auto_structure, manual_structure)
        >>> print(f"Sensitivity: {sens:.3f}")
    """
    if structure1.mask is None or structure2.mask is None:
        return 0.0

    intersection = np.logical_and(structure1.mask, structure2.mask)
    v2 = structure2.volume_voxels()

    if v2 == 0:
        return 0.0

    return float(np.sum(intersection) / v2)

compute_specificity ¶

compute_specificity(structure1: Structure, structure2: Structure, background_mask: Optional[ndarray] = None) -> float

Compute specificity (true negative rate).

Specificity = TN / (TN + FP)

Requires definition of background/universe. If not provided, uses the bounding box union of both structures.

Parameters:

Name	Type	Description	Default
`structure1`	`Structure`	Predicted/test structure	required
`structure2`	`Structure`	Reference/ground truth structure	required
`background_mask`	`Optional[ndarray]`	Mask defining the universe (optional)	`None`

Returns:

Type	Description
`float`	Specificity (0-1)

Examples:

>>> spec = compute_specificity(auto_structure, manual_structure)
>>> print(f"Specificity: {spec:.3f}")

Source code in src/dosemetrics/metrics/geometric.py

def compute_specificity(
    structure1: Structure, 
    structure2: Structure,
    background_mask: Optional[np.ndarray] = None
) -> float:
    """
    Compute specificity (true negative rate).

    Specificity = TN / (TN + FP)

    Requires definition of background/universe. If not provided, uses
    the bounding box union of both structures.

    Args:
        structure1: Predicted/test structure
        structure2: Reference/ground truth structure
        background_mask: Mask defining the universe (optional)

    Returns:
        Specificity (0-1)

    Examples:
        >>> spec = compute_specificity(auto_structure, manual_structure)
        >>> print(f"Specificity: {spec:.3f}")
    """
    if structure1.mask is None or structure2.mask is None:
        return 0.0

    # True negatives: voxels outside both structures
    # False positives: in structure1 but not in structure2
    not_s1 = ~structure1.mask
    not_s2 = ~structure2.mask

    true_negatives = np.logical_and(not_s1, not_s2)
    false_positives = np.logical_and(structure1.mask, not_s2)

    denominator = np.sum(true_negatives) + np.sum(false_positives)

    if denominator == 0:
        return 0.0

    return float(np.sum(true_negatives) / denominator)

compute_hausdorff_distance ¶

compute_hausdorff_distance(structure1: Structure, structure2: Structure, percentile: Optional[float] = None) -> float

Compute Hausdorff distance between two structures.

If percentile is specified, computes the percentile Hausdorff distance (e.g., 95th percentile HD95), which is more robust to outliers.

Parameters:

Name	Type	Description	Default
`structure1`	`Structure`	First structure	required
`structure2`	`Structure`	Second structure	required
`percentile`	`Optional[float]`	If specified, compute percentile HD (e.g., 95 for HD95)	`None`

Returns:

Type	Description
`float`	Hausdorff distance in mm

Examples:

>>> hd = compute_hausdorff_distance(auto_structure, manual_structure)
>>> hd95 = compute_hausdorff_distance(auto_structure, manual_structure, percentile=95)

Source code in src/dosemetrics/metrics/geometric.py

def compute_hausdorff_distance(
    structure1: Structure,
    structure2: Structure,
    percentile: Optional[float] = None
) -> float:
    """
    Compute Hausdorff distance between two structures.

    If percentile is specified, computes the percentile Hausdorff distance
    (e.g., 95th percentile HD95), which is more robust to outliers.

    Args:
        structure1: First structure
        structure2: Second structure
        percentile: If specified, compute percentile HD (e.g., 95 for HD95)

    Returns:
        Hausdorff distance in mm

    Examples:
        >>> hd = compute_hausdorff_distance(auto_structure, manual_structure)
        >>> hd95 = compute_hausdorff_distance(auto_structure, manual_structure, percentile=95)
    """
    # Validate percentile
    if percentile is not None:
        if not (0 < percentile <= 100):
            raise ValueError(f"Percentile must be between 0 and 100, got {percentile}")

    if structure1.mask is None or structure2.mask is None:
        return float('inf')

    # Get surface points (boundary voxels)
    # Use binary erosion to get boundary
    eroded1 = ndimage.binary_erosion(structure1.mask)
    eroded2 = ndimage.binary_erosion(structure2.mask)
    surface1 = structure1.mask & ~eroded1
    surface2 = structure2.mask & ~eroded2

    # Get coordinates of surface points
    points1 = np.argwhere(surface1)
    points2 = np.argwhere(surface2)

    if len(points1) == 0 or len(points2) == 0:
        return float('inf')

    # Scale by voxel spacing to get mm
    spacing = np.array(structure1.spacing)
    points1_mm = points1 * spacing
    points2_mm = points2 * spacing

    if percentile is not None:
        # Compute percentile Hausdorff distance
        # Calculate distances from points1 to points2
        from scipy.spatial.distance import cdist
        distances = cdist(points1_mm, points2_mm)

        # For each point in set 1, find min distance to set 2
        min_distances_1_to_2 = np.min(distances, axis=1)
        # For each point in set 2, find min distance to set 1
        min_distances_2_to_1 = np.min(distances, axis=0)

        # Compute percentile
        hd_1_to_2 = np.percentile(min_distances_1_to_2, percentile)
        hd_2_to_1 = np.percentile(min_distances_2_to_1, percentile)

        return float(max(hd_1_to_2, hd_2_to_1))
    else:
        # Standard Hausdorff distance
        hd_1_to_2, _, _ = directed_hausdorff(points1_mm, points2_mm)
        hd_2_to_1, _, _ = directed_hausdorff(points2_mm, points1_mm)

        return float(max(hd_1_to_2, hd_2_to_1))

compute_mean_surface_distance ¶

compute_mean_surface_distance(structure1: Structure, structure2: Structure) -> float

Compute mean surface distance between two structures.

Average of all point-to-surface distances (symmetric).

Parameters:

Name	Type	Description	Default
`structure1`	`Structure`	First structure	required
`structure2`	`Structure`	Second structure	required

Returns:

Type	Description
`float`	Mean surface distance in mm

Source code in src/dosemetrics/metrics/geometric.py

def compute_mean_surface_distance(
    structure1: Structure,
    structure2: Structure
) -> float:
    """
    Compute mean surface distance between two structures.

    Average of all point-to-surface distances (symmetric).

    Args:
        structure1: First structure
        structure2: Second structure

    Returns:
        Mean surface distance in mm
    """
    if structure1.mask is None or structure2.mask is None:
        return float('inf')

    # Get surface points (boundary voxels)
    eroded1 = ndimage.binary_erosion(structure1.mask)
    eroded2 = ndimage.binary_erosion(structure2.mask)
    surface1 = structure1.mask & ~eroded1
    surface2 = structure2.mask & ~eroded2

    # Get coordinates of surface points
    points1 = np.argwhere(surface1)
    points2 = np.argwhere(surface2)

    if len(points1) == 0 or len(points2) == 0:
        return float('inf')

    # Scale by voxel spacing to get mm
    spacing = np.array(structure1.spacing)
    points1_mm = points1 * spacing
    points2_mm = points2 * spacing

    # Compute pairwise distances
    from scipy.spatial.distance import cdist
    distances = cdist(points1_mm, points2_mm)

    # Mean of minimum distances from each point to other surface
    mean_1_to_2 = np.mean(np.min(distances, axis=1))
    mean_2_to_1 = np.mean(np.min(distances, axis=0))

    # Return symmetric average
    return float((mean_1_to_2 + mean_2_to_1) / 2.0)

compare_structure_sets ¶

compare_structure_sets(structure_set1: StructureSet, structure_set2: StructureSet, structure_names: Optional[list] = None) -> pd.DataFrame

Compute geometric metrics between two structure sets.

Parameters:

Name	Type	Description	Default
`structure_set1`	`StructureSet`	First structure set (e.g., auto-segmentation)	required
`structure_set2`	`StructureSet`	Second structure set (e.g., manual segmentation)	required
`structure_names`	`Optional[list]`	List of structure names to compare (optional)	`None`

Returns:

Type	Description
`DataFrame`	DataFrame with geometric metrics for each structure

Examples:

>>> auto_structures = load_structure_set("auto/")
>>> manual_structures = load_structure_set("manual/")
>>> comparison = compare_structure_sets(auto_structures, manual_structures)
>>> print(comparison)

Source code in src/dosemetrics/metrics/geometric.py

def compare_structure_sets(
    structure_set1: StructureSet,
    structure_set2: StructureSet,
    structure_names: Optional[list] = None
) -> pd.DataFrame:
    """
    Compute geometric metrics between two structure sets.

    Args:
        structure_set1: First structure set (e.g., auto-segmentation)
        structure_set2: Second structure set (e.g., manual segmentation)
        structure_names: List of structure names to compare (optional)

    Returns:
        DataFrame with geometric metrics for each structure

    Examples:
        >>> auto_structures = load_structure_set("auto/")
        >>> manual_structures = load_structure_set("manual/")
        >>> comparison = compare_structure_sets(auto_structures, manual_structures)
        >>> print(comparison)
    """
    if structure_names is None:
        # Use common structures
        names1 = set(structure_set1.structure_names)
        names2 = set(structure_set2.structure_names)
        structure_names = list(names1.intersection(names2))

    results = []

    for name in structure_names:
        try:
            struct1 = structure_set1.get_structure(name)
            struct2 = structure_set2.get_structure(name)

            dice = compute_dice_coefficient(struct1, struct2)
            jaccard = compute_jaccard_index(struct1, struct2)
            vol_diff = compute_volume_difference(struct1, struct2)
            vol_ratio = compute_volume_ratio(struct1, struct2)
            sensitivity = compute_sensitivity(struct1, struct2)

            results.append({
                'Structure': name,
                'Dice': dice,
                'Jaccard': jaccard,
                'Volume_Difference_cc': vol_diff,
                'Volume_Ratio': vol_ratio,
                'Sensitivity': sensitivity,
            })
        except ValueError:
            # Structure not found in one of the sets
            continue

    return pd.DataFrame(results)

`homogeneity`¶

homogeneity ¶

Homogeneity indices for target dose uniformity.

This module provides metrics to assess the uniformity of dose distribution within target volumes. More homogeneous dose distributions are generally preferred for tumor control.

Classes¶

Functions:¶

compute_homogeneity_index ¶

compute_homogeneity_index(dose: Dose, target: Structure, d2_percentile: float = 2.0, d98_percentile: float = 98.0) -> float

Compute Homogeneity Index (HI).

HI = (D2 - D98) / D50

Where: - D2 = dose received by 2% of volume (near-maximum) - D98 = dose received by 98% of volume (near-minimum) - D50 = median dose

Measures dose uniformity within target. Lower values indicate more homogeneous dose distribution.

Typical acceptable range: 0.05 - 0.20

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`target`	`Structure`	Target structure (PTV, CTV, etc.)	required
`d2_percentile`	`float`	Upper percentile for near-max (typically 2%)	`2.0`
`d98_percentile`	`float`	Lower percentile for near-min (typically 98%)	`98.0`

Returns:

Type	Description
`float`	Homogeneity index (dimensionless)

References

ICRU Report 83 (2010)

Examples:

>>> hi = compute_homogeneity_index(dose, ptv)
>>> print(f"Homogeneity Index: {hi:.3f}")
>>> if hi < 0.15:
...     print("Excellent dose homogeneity")

Source code in src/dosemetrics/metrics/homogeneity.py

def compute_homogeneity_index(
    dose: Dose,
    target: Structure,
    d2_percentile: float = 2.0,
    d98_percentile: float = 98.0
) -> float:
    """
    Compute Homogeneity Index (HI).

    HI = (D2 - D98) / D50

    Where:
    - D2 = dose received by 2% of volume (near-maximum)
    - D98 = dose received by 98% of volume (near-minimum)
    - D50 = median dose

    Measures dose uniformity within target. Lower values indicate more
    homogeneous dose distribution.

    Typical acceptable range: 0.05 - 0.20

    Args:
        dose: Dose distribution object
        target: Target structure (PTV, CTV, etc.)
        d2_percentile: Upper percentile for near-max (typically 2%)
        d98_percentile: Lower percentile for near-min (typically 98%)

    Returns:
        Homogeneity index (dimensionless)

    References:
        ICRU Report 83 (2010)

    Examples:
        >>> hi = compute_homogeneity_index(dose, ptv)
        >>> print(f"Homogeneity Index: {hi:.3f}")
        >>> if hi < 0.15:
        ...     print("Excellent dose homogeneity")
    """
    dose_values = dose.get_dose_in_structure(target)

    if len(dose_values) == 0:
        return 0.0

    # Note: D2 means 2% of volume receives at least this dose
    # This corresponds to 98th percentile of dose array
    d2 = np.percentile(dose_values, 100 - d2_percentile)
    d98 = np.percentile(dose_values, 100 - d98_percentile)
    d50 = np.percentile(dose_values, 50)

    if d50 == 0:
        return float('inf')

    return float((d2 - d98) / d50)

compute_gradient_index ¶

compute_gradient_index(dose: Dose, target: Structure, prescription_dose: float, half_prescription_volume_method: bool = True) -> float

Compute Gradient Index (GI) for dose fall-off outside target.

The implementation uses the half-prescription-volume definition: GI = V_50% / V_100%.

Where: - V_100% = volume receiving >= prescription dose - V_50% = volume receiving >= 50% prescription dose

Lower values indicate steeper dose fall-off (better for sparing OARs).

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`target`	`Structure`	Target structure	required
`prescription_dose`	`float`	Prescription dose in Gy	required
`half_prescription_volume_method`	`bool`	Compatibility parameter. The current implementation always uses V_50%/V_100%.	`True`

Returns:

Type	Description
`float`	Gradient index (dimensionless, typically 2-8)

References

Paddick and Lippitz, J Neurosurg 2006

Examples:

>>> gi = compute_gradient_index(dose, ptv, prescription_dose=60.0)
>>> print(f"Gradient Index: {gi:.2f}")
>>> if gi < 3.0:
...     print("Excellent dose fall-off")

Source code in src/dosemetrics/metrics/homogeneity.py

def compute_gradient_index(
    dose: Dose,
    target: Structure,
    prescription_dose: float,
    half_prescription_volume_method: bool = True
) -> float:
    """
    Compute Gradient Index (GI) for dose fall-off outside target.

    The implementation uses the half-prescription-volume definition:
    GI = V_50% / V_100%.

    Where:
    - V_100% = volume receiving >= prescription dose
    - V_50% = volume receiving >= 50% prescription dose

    Lower values indicate steeper dose fall-off (better for sparing OARs).

    Args:
        dose: Dose distribution object
        target: Target structure
        prescription_dose: Prescription dose in Gy
        half_prescription_volume_method: Compatibility parameter. The current
            implementation always uses V_50%/V_100%.

    Returns:
        Gradient index (dimensionless, typically 2-8)

    References:
        Paddick and Lippitz, J Neurosurg 2006

    Examples:
        >>> gi = compute_gradient_index(dose, ptv, prescription_dose=60.0)
        >>> print(f"Gradient Index: {gi:.2f}")
        >>> if gi < 3.0:
        ...     print("Excellent dose fall-off")
    """
    v_100 = np.sum(dose.dose_array >= prescription_dose)
    v_50 = np.sum(dose.dose_array >= 0.5 * prescription_dose)

    if v_100 == 0:
        return float('inf')

    return float(v_50 / v_100)

compute_dose_homogeneity ¶

compute_dose_homogeneity(dose: Dose, target: Structure) -> float

Compute coefficient of variation (CV) of dose within target.

CV = std_dose / mean_dose

Alternative measure of dose homogeneity. Lower values indicate more uniform dose distribution.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`target`	`Structure`	Target structure	required

Returns:

Type	Description
`float`	Coefficient of variation (dimensionless)

Examples:

>>> cv = compute_dose_homogeneity(dose, ptv)
>>> print(f"Dose CV: {cv:.3f}")

Source code in src/dosemetrics/metrics/homogeneity.py

def compute_dose_homogeneity(
    dose: Dose,
    target: Structure
) -> float:
    """
    Compute coefficient of variation (CV) of dose within target.

    CV = std_dose / mean_dose

    Alternative measure of dose homogeneity. Lower values indicate
    more uniform dose distribution.

    Args:
        dose: Dose distribution object
        target: Target structure

    Returns:
        Coefficient of variation (dimensionless)

    Examples:
        >>> cv = compute_dose_homogeneity(dose, ptv)
        >>> print(f"Dose CV: {cv:.3f}")
    """
    dose_values = dose.get_dose_in_structure(target)

    if len(dose_values) == 0:
        return 0.0

    mean = np.mean(dose_values)
    if mean == 0:
        return float('inf')

    std = np.std(dose_values)
    return float(std / mean)

compute_uniformity_index ¶

compute_uniformity_index(dose: Dose, target: Structure) -> float

Compute uniformity index.

UI = 1 - (D_max - D_min) / D_ref

Values closer to 1.0 indicate better uniformity.

Parameters:

Name	Type	Description	Default
`dose`	`Dose`	Dose distribution object	required
`target`	`Structure`	Target structure	required

Returns:

Type	Description
`float`	Uniformity index (0-1)

Note

The current implementation uses median target dose as D_ref.

Examples:

>>> ui = compute_uniformity_index(dose, ptv)
>>> print(f"Uniformity Index: {ui:.3f}")

Source code in src/dosemetrics/metrics/homogeneity.py

def compute_uniformity_index(
    dose: Dose,
    target: Structure
) -> float:
    """
    Compute uniformity index.

    UI = 1 - (D_max - D_min) / D_ref

    Values closer to 1.0 indicate better uniformity.

    Args:
        dose: Dose distribution object
        target: Target structure

    Returns:
        Uniformity index (0-1)

    Note:
        The current implementation uses median target dose as D_ref.

    Examples:
        >>> ui = compute_uniformity_index(dose, ptv)
        >>> print(f"Uniformity Index: {ui:.3f}")
    """
    dose_values = dose.get_dose_in_structure(target)

    if len(dose_values) == 0:
        return 0.0

    d_max = np.max(dose_values)
    d_min = np.min(dose_values)
    d_ref = np.median(dose_values)  # Use median as reference

    if d_ref == 0:
        return 0.0

    return float(1.0 - (d_max - d_min) / d_ref)

Metrics API¶

Direct plan-comparison functions¶

compare_ptv_dose ¶

compare_paddick_conformity_index ¶

compare_paddick_gradient_index ¶

compare_homogeneity_index ¶

compare_body_rmse ¶

compare_gamma ¶

compare_oar_constraints ¶

compare_oar_dvh_auc ¶

compare_mean_oar_dvh_auc ¶

compare_dvh_score ¶

Comparison metadata¶

comparison ¶

Attributes¶

COMPARISON_METRICS module-attribute ¶

Classes¶

EvaluationTask ¶

MetricCategory ¶

MetricDefinition dataclass ¶

conformity¶

conformity ¶

Classes¶

Functions:¶

compute_conformity_index ¶

compute_conformity_number ¶

compute_paddick_conformity_index ¶

compute_coverage ¶

compute_spillage ¶

compute_rtog_conformity_index ¶

compute_prescription_mae ¶

dose_comparison¶

dose_comparison ¶

Classes¶

Functions:¶

compare_ssim ¶

Parameters¶

Returns¶

Notes¶

References¶

Examples¶

Raises¶

compare_mse ¶

Parameters¶

Returns¶

Raises¶

compare_mae ¶

Parameters¶

Returns¶

Notes¶

Raises¶

compare_psnr ¶

Parameters¶

Returns¶

Notes¶

Raises¶

compare_mutual_information ¶

Parameters¶

Returns¶

Notes¶

Raises¶

compare_normalized_cross_correlation ¶

Parameters¶

Returns¶

Notes¶

Raises¶

compare_dose_difference_map ¶

Parameters¶

Returns¶

Notes¶

Raises¶

compare_dose ¶

Parameters¶

Returns¶

Examples¶

Raises¶

compute_3d_dose_gradient ¶

Parameters¶

Returns¶

Notes¶

COMPARISON_METRICS `module-attribute` ¶

MetricDefinition `dataclass` ¶

`conformity`¶

`dose_comparison`¶

`dvh`¶

`gamma`¶

`geometric`¶

`homogeneity`¶