Metrics API¶

This module contains core functions for dose analysis, DVH computation, quality metrics, and plan comparison.

DVH Module¶

dvh ¶

Dose-Volume Histogram (DVH) computation and analysis.

This module provides functions for computing DVHs and extracting DVH-based metrics such as volume at dose (VX) and dose at volume (DX).

Classes¶

Functions:¶

compute_dvh ¶

compute_dvh(dose: Dose, structure: Structure, max_dose: Optional[float] = None, step_size: float = 0.1) -> 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`

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)
>>>
>>> # 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,
) -> 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

    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)
        >>>
        >>> # 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])
        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
        ]
    )

    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) -> 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

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)
>>> print(f"Mean dose: {stats['mean_dose']:.2f} Gy")
>>> print(f"D95: {stats['D95']:.2f} Gy")

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

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

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

    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)
        >>> print(f"Mean dose: {stats['mean_dose']:.2f} Gy")
        >>> print(f"D95: {stats['D95']:.2f} Gy")
    """
    dose_values = dose.get_dose_in_structure(structure)

    if len(dose_values) == 0:
        return {
            "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,
        }

    return {
        "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
    }

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

compute_dvh_score ¶

compute_dvh_score(dose_reference: Dose, dose_evaluated: Dose, structure: Structure) -> float

Compute DVH Score: average absolute difference in D1, D95, D99 between two dose distributions.

DVH Score = (|D1_ref - D1_eval| + |D95_ref - D95_eval| + |D99_ref - D99_eval|) / 3

Where DX means X% of the volume receives at least this dose. This metric captures clinically relevant dose differences at near-maximum (D1), near-minimum (D99), and dose coverage (D95) levels.

Lower values indicate better agreement between the two distributions.

Parameters:

Name	Type	Description	Default
`dose_reference`	`Dose`	Reference dose distribution	required
`dose_evaluated`	`Dose`	Evaluated dose distribution to compare	required
`structure`	`Structure`	Structure to restrict comparison to	required

Returns:

Type	Description
`float`	Average absolute DVH difference in Gy

References

Adapted from GDP-HMM AAPM Challenge evaluation methodology.

Examples:

>>> score = compute_dvh_score(reference_dose, predicted_dose, ptv)
>>> print(f"DVH Score: {score:.3f} Gy")

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

def compute_dvh_score(
    dose_reference: Dose,
    dose_evaluated: Dose,
    structure: Structure,
) -> float:
    """
    Compute DVH Score: average absolute difference in D1, D95, D99 between two dose distributions.

    DVH Score = (|D1_ref - D1_eval| + |D95_ref - D95_eval| + |D99_ref - D99_eval|) / 3

    Where DX means X% of the volume receives at least this dose. This metric
    captures clinically relevant dose differences at near-maximum (D1),
    near-minimum (D99), and dose coverage (D95) levels.

    Lower values indicate better agreement between the two distributions.

    Args:
        dose_reference: Reference dose distribution
        dose_evaluated: Evaluated dose distribution to compare
        structure: Structure to restrict comparison to

    Returns:
        Average absolute DVH difference in Gy

    References:
        Adapted from GDP-HMM AAPM Challenge evaluation methodology.

    Examples:
        >>> score = compute_dvh_score(reference_dose, predicted_dose, ptv)
        >>> print(f"DVH Score: {score:.3f} Gy")
    """
    ref_values = dose_reference.get_dose_in_structure(structure)
    eval_values = dose_evaluated.get_dose_in_structure(structure)

    if len(ref_values) == 0 or len(eval_values) == 0:
        return float("nan")

    # DX: X% of volume receives at least this dose → (100-X)th percentile of dose array
    d1_ref = float(np.percentile(ref_values, 99))  # D1: near-max
    d95_ref = float(np.percentile(ref_values, 5))  # D95: coverage
    d99_ref = float(np.percentile(ref_values, 1))  # D99: near-min

    d1_eval = float(np.percentile(eval_values, 99))
    d95_eval = float(np.percentile(eval_values, 5))
    d99_eval = float(np.percentile(eval_values, 1))

    return (
        abs(d1_ref - d1_eval) + abs(d95_ref - d95_eval) + abs(d99_ref - d99_eval)
    ) / 3.0

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 value in Gy (dose × volume units).

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 value in Gy (dose × volume units).

    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))

Conformity Module¶

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)))

Homogeneity Module¶

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.

Two calculation methods: 1. Half-prescription volume: GI = V_50% / V_100% 2. Distance-based: Ratio of volumes at specific distances

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`	Use V_50%/V_100% method (default True)	`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.

    Two calculation methods:
    1. Half-prescription volume: GI = V_50% / V_100%
    2. Distance-based: Ratio of volumes at specific distances

    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: Use V_50%/V_100% method (default True)

    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_prescription

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

Requires prescription dose in target metadata or as parameter. Currently uses median dose as approximation.

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_prescription

    Values closer to 1.0 indicate better uniformity.

    Args:
        dose: Dose distribution object
        target: Target structure

    Returns:
        Uniformity index (0-1)

    Note:
        Requires prescription dose in target metadata or as parameter.
        Currently uses median dose as approximation.

    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)

Geometric Module¶

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)

Gamma Module¶

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:¶

compute_gamma_index ¶

compute_gamma_index(dose_reference: Dose, 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¶

dose_reference : Dose Reference (planned) dose distribution. dose_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 = compute_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 compute_gamma_index(
    dose_reference: Dose,
    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
    ----------
    dose_reference : Dose
        Reference (planned) dose distribution.
    dose_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 = compute_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 dose_reference.dose_array.shape != dose_evaluated.dose_array.shape:
        raise ValueError(
            f"Dose shapes must match: {dose_reference.dose_array.shape} vs "
            f"{dose_evaluated.dose_array.shape}"
        )

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

    # Get dose arrays and spatial information
    ref_dose = dose_reference.dose_array
    eval_dose = dose_evaluated.dose_array
    spacing = np.array(dose_reference.spacing)
    origin = np.array(dose_reference.origin)
    shape = dose_reference.shape

    # 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 = None  # Will use local normalization

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

    # Create coordinate grids for both distributions
    x = origin[0] + np.arange(shape[0]) * spacing[0]
    y = origin[1] + np.arange(shape[1]) * spacing[1]
    z = origin[2] + np.arange(shape[2]) * spacing[2]

    # Create interpolator for evaluated dose
    eval_interpolator = RegularGridInterpolator(
        (x, y, z), eval_dose, method="linear", bounds_error=False, fill_value=0.0
    )

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

    # 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]

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

                # Skip if below threshold
                if ref_value < dose_threshold:
                    continue

                # Use local normalization if requested
                if not global_normalization:
                    normalization_dose = ref_value if ref_value > 0 else 1.0

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

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

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

                if len(i_search) == 0:
                    continue

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

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

                # Compute gamma values using the Low et al. formula
                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, k] = np.min(gamma_values)

    return gamma_result

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 compute_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 compute_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

compute_2d_gamma ¶

compute_2d_gamma(dose_reference_slice: ndarray, dose_evaluated_slice: 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¶

dose_reference_slice : np.ndarray 2D reference dose slice. dose_evaluated_slice : 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 compute_2d_gamma(
    dose_reference_slice: np.ndarray,
    dose_evaluated_slice: 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
    ----------
    dose_reference_slice : np.ndarray
        2D reference dose slice.
    dose_evaluated_slice : 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 dose_reference_slice.ndim != 2:
        raise ValueError(
            f"Reference slice must be 2D, got shape {dose_reference_slice.shape}"
        )
    if dose_evaluated_slice.ndim != 2:
        raise ValueError(
            f"Evaluated slice must be 2D, got shape {dose_evaluated_slice.shape}"
        )
    if dose_reference_slice.shape != dose_evaluated_slice.shape:
        raise ValueError(
            f"Slice shapes must match: {dose_reference_slice.shape} vs "
            f"{dose_evaluated_slice.shape}"
        )

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

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

    # Create interpolator for evaluated dose
    x = np.arange(shape[0]) * spacing[0]
    y = np.arange(shape[1]) * spacing[1]
    eval_interpolator = RegularGridInterpolator(
        (x, y),
        dose_evaluated_slice,
        method="linear",
        bounds_error=False,
        fill_value=0.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 = dose_reference_slice[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 = dose_evaluated_slice[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

compute_gamma_index_gpu ¶

compute_gamma_index_gpu(dose_reference: Dose, 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¶

dose_reference : Dose Reference dose. 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 compute_gamma_index_gpu(
    dose_reference: Dose,
    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
    ----------
    dose_reference : Dose
        Reference dose.
    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 compute_gamma_index() for CPU-based calculation.",
        FutureWarning,
    )
    raise NotImplementedError(
        "GPU-accelerated gamma not implemented. Use compute_gamma_index()."
    )

Dose Comparison Module¶

dose_comparison ¶

Dose distribution comparison metrics beyond DVH.

This module provides image-based metrics for comparing 3D dose distributions, including SSIM, MSE, MAE, and other similarity measures.

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:¶

compute_ssim ¶

compute_ssim(dose1: Dose, dose2: 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¶

dose1 : Dose Reference dose distribution. dose2 : Dose Comparison 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 = compute_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 compute_ssim(
    dose1: Dose,
    dose2: 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
    ----------
    dose1 : Dose
        Reference dose distribution.
    dose2 : Dose
        Comparison 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 = compute_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 = dose1.dose_array
    arr2 = dose2.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)

compute_mse ¶

compute_mse(dose1: Dose, dose2: Dose, structure: Optional[Structure] = None) -> float

Compute Mean Squared Error between two dose distributions.

Parameters¶

dose1 : Dose Reference dose. dose2 : Dose Comparison 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 compute_mse(
    dose1: Dose, dose2: Dose, structure: Optional[Structure] = None
) -> float:
    """
    Compute Mean Squared Error between two dose distributions.

    Parameters
    ----------
    dose1 : Dose
        Reference dose.
    dose2 : Dose
        Comparison 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 = dose1.dose_array
    arr2 = dose2.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)

compute_mae ¶

compute_mae(dose1: Dose, dose2: Dose, structure: Optional[Structure] = None) -> float

Compute Mean Absolute Error between two dose distributions.

Parameters¶

dose1 : Dose Reference dose. dose2 : Dose Comparison 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 compute_mae(
    dose1: Dose, dose2: Dose, structure: Optional[Structure] = None
) -> float:
    """
    Compute Mean Absolute Error between two dose distributions.

    Parameters
    ----------
    dose1 : Dose
        Reference dose.
    dose2 : Dose
        Comparison 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 = dose1.dose_array
    arr2 = dose2.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)

compute_psnr ¶

compute_psnr(dose1: Dose, dose2: Dose, structure: Optional[Structure] = None, data_range: Optional[float] = None) -> float

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

Parameters¶

dose1 : Dose Reference dose. dose2 : Dose Comparison 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 compute_psnr(
    dose1: Dose,
    dose2: Dose,
    structure: Optional[Structure] = None,
    data_range: Optional[float] = None,
) -> float:
    """
    Compute Peak Signal-to-Noise Ratio between two dose distributions.

    Parameters
    ----------
    dose1 : Dose
        Reference dose.
    dose2 : Dose
        Comparison 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 = compute_mse(dose1, dose2, structure)

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

    # Determine data range
    if data_range is None:
        arr1 = dose1.dose_array
        arr2 = dose2.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)

compute_mutual_information ¶

compute_mutual_information(dose1: Dose, dose2: Dose, structure: Optional[Structure] = None, bins: int = 256) -> float

Compute Mutual Information between two dose distributions.

Parameters¶

dose1 : Dose First dose distribution. dose2 : Dose Second 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 compute_mutual_information(
    dose1: Dose, dose2: Dose, structure: Optional[Structure] = None, bins: int = 256
) -> float:
    """
    Compute Mutual Information between two dose distributions.

    Parameters
    ----------
    dose1 : Dose
        First dose distribution.
    dose2 : Dose
        Second 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 = dose1.dose_array.flatten()
    arr2 = dose2.dose_array.flatten()

    # Check shapes match
    if arr1.shape != arr2.shape:
        raise ValueError(f"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)

compute_normalized_cross_correlation ¶

compute_normalized_cross_correlation(dose1: Dose, dose2: Dose, structure: Optional[Structure] = None) -> float

Compute Normalized Cross-Correlation between two dose distributions.

Parameters¶

dose1 : Dose First dose distribution. dose2 : Dose Second 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 compute_normalized_cross_correlation(
    dose1: Dose, dose2: Dose, structure: Optional[Structure] = None
) -> float:
    """
    Compute Normalized Cross-Correlation between two dose distributions.

    Parameters
    ----------
    dose1 : Dose
        First dose distribution.
    dose2 : Dose
        Second 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 = dose1.dose_array.flatten()
    arr2 = dose2.dose_array.flatten()

    # Check shapes match
    if arr1.shape != arr2.shape:
        raise ValueError(f"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)

compute_dose_difference_map ¶

compute_dose_difference_map(dose1: Dose, dose2: Dose, absolute: bool = False) -> Dose

Compute voxel-wise dose difference map.

Parameters¶

dose1 : Dose Reference dose. dose2 : Dose Comparison 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 compute_dose_difference_map(
    dose1: Dose, dose2: Dose, absolute: bool = False
) -> Dose:
    """
    Compute voxel-wise dose difference map.

    Parameters
    ----------
    dose1 : Dose
        Reference dose.
    dose2 : Dose
        Comparison 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 dose1.dose_array.shape != dose2.dose_array.shape:
        raise ValueError(
            f"Dose shapes must match: {dose1.dose_array.shape} vs {dose2.dose_array.shape}"
        )

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

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

    return diff_dose

compute_dose_comparison_metrics ¶

compute_dose_comparison_metrics(dose1: Dose, dose2: Dose, structure: Optional[Structure] = None) -> Dict[str, float]

Compute comprehensive set of dose comparison metrics.

Parameters¶

dose1 : Dose Reference dose. dose2 : Dose Comparison 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 = compute_dose_comparison_metrics(dose1, dose2, 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 compute_dose_comparison_metrics(
    dose1: Dose, dose2: Dose, structure: Optional[Structure] = None
) -> Dict[str, float]:
    """
    Compute comprehensive set of dose comparison metrics.

    Parameters
    ----------
    dose1 : Dose
        Reference dose.
    dose2 : Dose
        Comparison 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 = compute_dose_comparison_metrics(dose1, dose2, 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"] = compute_mse(dose1, dose2, structure)
    except Exception as e:
        warnings.warn(f"MSE computation failed: {e}")
        metrics["mse"] = np.nan

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

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

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

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

    try:
        metrics["mi"] = compute_mutual_information(dose1, dose2, 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")

compute_normalized_mae ¶

compute_normalized_mae(dose_reference: Dose, 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
`dose_reference`	`Dose`	Reference dose distribution	required
`dose_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 compute_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 = compute_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 compute_normalized_mae(
    dose_reference: Dose,
    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:
        dose_reference: Reference dose distribution
        dose_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 compute_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 = compute_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 = dose_reference.get_dose_in_structure(structure)
        eval_arr = dose_evaluated.get_dose_in_structure(structure)
    else:
        ref_arr = dose_reference.dose_array.flatten()
        eval_arr = dose_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

Advanced DVH Module¶

advanced_dvh ¶

Advanced DVH metrics and comparison tools.

This module provides advanced DVH-based metrics for comparing dose distributions including Wasserstein distance, area between curves, and other statistical measures.

Future Implementation TODOs

Wasserstein distance (Earth Mover's Distance) between DVHs
Area between DVH curves (L1/L2 norms)
DVH bandwidth and confidence intervals
Chi-square and Kolmogorov-Smirnov tests for DVH comparison
DVH-based TCP/NTCP models

Classes¶

Functions:¶

compute_dvh_wasserstein_distance ¶

compute_dvh_wasserstein_distance(dose1: Dose, dose2: Dose, structure: Structure) -> float

Compute Wasserstein distance (Earth Mover's Distance) between two DVHs.

The Wasserstein distance quantifies the minimum "work" required to transform one DVH into another, providing a meaningful metric for DVH similarity.

Parameters¶

dose1 : Dose First dose distribution. dose2 : Dose Second dose distribution. structure : Structure Structure for which to compute DVHs.

Returns¶

distance : float Wasserstein distance between the two DVHs.

Notes¶

The Wasserstein distance is also known as: - Earth Mover's Distance (EMD) - Kantorovich-Rubinstein metric - Mallows distance

It satisfies the triangle inequality and is a true metric, unlike simple area-between-curves measures.

References¶

Rubner Y, Tomasi C, Guibas LJ. "The Earth Mover's Distance as a Metric for Image Retrieval." Int J Comput Vision. 2000;40(2):99-121.

Examples¶

from dosemetrics.metrics import advanced_dvh distance = advanced_dvh.compute_dvh_wasserstein_distance( ... planned_dose, delivered_dose, ptv ... ) print(f"DVH Wasserstein distance: {distance:.2f} Gy")

Raises¶

NotImplementedError This function is a stub for future implementation.

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

def compute_dvh_wasserstein_distance(
    dose1: Dose,
    dose2: Dose,
    structure: Structure
) -> float:
    """
    Compute Wasserstein distance (Earth Mover's Distance) between two DVHs.

    The Wasserstein distance quantifies the minimum "work" required to transform
    one DVH into another, providing a meaningful metric for DVH similarity.

    Parameters
    ----------
    dose1 : Dose
        First dose distribution.
    dose2 : Dose
        Second dose distribution.
    structure : Structure
        Structure for which to compute DVHs.

    Returns
    -------
    distance : float
        Wasserstein distance between the two DVHs.

    Notes
    -----
    The Wasserstein distance is also known as:
        - Earth Mover's Distance (EMD)
        - Kantorovich-Rubinstein metric
        - Mallows distance

    It satisfies the triangle inequality and is a true metric, unlike
    simple area-between-curves measures.

    References
    ----------
    - Rubner Y, Tomasi C, Guibas LJ. "The Earth Mover's Distance as a Metric
      for Image Retrieval." Int J Comput Vision. 2000;40(2):99-121.

    Examples
    --------
    >>> from dosemetrics.metrics import advanced_dvh
    >>> distance = advanced_dvh.compute_dvh_wasserstein_distance(
    ...     planned_dose, delivered_dose, ptv
    ... )
    >>> print(f"DVH Wasserstein distance: {distance:.2f} Gy")

    Raises
    ------
    NotImplementedError
        This function is a stub for future implementation.
    """
    # Get dose values within structure for both doses
    dose_values1 = dose1.get_dose_in_structure(structure)
    dose_values2 = dose2.get_dose_in_structure(structure)

    if len(dose_values1) == 0 or len(dose_values2) == 0:
        return 0.0

    # Compute Wasserstein distance directly on dose values
    distance = wasserstein_distance(dose_values1, dose_values2)

    return float(distance)

compute_area_between_dvh_curves ¶

compute_area_between_dvh_curves(dose1: Dose, dose2: Dose, structure: Structure, norm: str = 'L2') -> float

Compute area between two DVH curves.

Parameters¶

dose1 : Dose First dose distribution. dose2 : Dose Second dose distribution. structure : Structure Structure for which to compute DVHs. norm : {'L1', 'L2'}, optional Norm to use for area calculation: - 'L1': Sum of absolute differences - 'L2': Sum of squared differences (default)

Returns¶

area : float Area between the two DVH curves.

Notes¶

The L1 norm gives the Manhattan distance, while L2 gives Euclidean distance. For DVH comparison, L1 is often more interpretable.

Raises¶

ValueError If norm is not 'L1' or 'L2'.

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

def compute_area_between_dvh_curves(
    dose1: Dose,
    dose2: Dose,
    structure: Structure,
    norm: str = 'L2'
) -> float:
    """
    Compute area between two DVH curves.

    Parameters
    ----------
    dose1 : Dose
        First dose distribution.
    dose2 : Dose
        Second dose distribution.
    structure : Structure
        Structure for which to compute DVHs.
    norm : {'L1', 'L2'}, optional
        Norm to use for area calculation:
            - 'L1': Sum of absolute differences
            - 'L2': Sum of squared differences (default)

    Returns
    -------
    area : float
        Area between the two DVH curves.

    Notes
    -----
    The L1 norm gives the Manhattan distance, while L2 gives Euclidean distance.
    For DVH comparison, L1 is often more interpretable.

    Raises
    ------
    ValueError
        If norm is not 'L1' or 'L2'.
    """
    if norm not in ['L1', 'L2']:
        raise ValueError(f"norm must be 'L1' or 'L2', got '{norm}'")

    # Compute DVHs
    dose_bins1, volumes1 = compute_dvh(dose1, structure)
    dose_bins2, volumes2 = compute_dvh(dose2, structure)

    # Create common dose bins
    max_dose = max(np.max(dose_bins1), np.max(dose_bins2))
    step_size = min(
        dose_bins1[1] - dose_bins1[0] if len(dose_bins1) > 1 else 0.1,
        dose_bins2[1] - dose_bins2[0] if len(dose_bins2) > 1 else 0.1
    )
    common_bins = np.arange(0, max_dose + step_size, step_size)

    # Interpolate to common bins
    volumes1_interp = np.interp(common_bins, dose_bins1, volumes1)
    volumes2_interp = np.interp(common_bins, dose_bins2, volumes2)

    # Compute area based on norm
    if norm == 'L1':
        area = np.sum(np.abs(volumes1_interp - volumes2_interp)) * step_size
    else:  # L2
        area = np.sqrt(np.sum((volumes1_interp - volumes2_interp) ** 2)) * step_size

    return float(area)

compute_dvh_chi_square ¶

compute_dvh_chi_square(dose1: Dose, dose2: Dose, structure: Structure) -> Tuple[float, float]

Perform chi-square test comparing two DVHs.

Parameters¶

dose1 : Dose First (expected) dose distribution. dose2 : Dose Second (observed) dose distribution. structure : Structure Structure for DVH computation.

Returns¶

chi2_statistic : float Chi-square test statistic. p_value : float P-value for the test.

Notes¶

Tests the null hypothesis that the two DVHs come from the same distribution.

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

def compute_dvh_chi_square(
    dose1: Dose,
    dose2: Dose,
    structure: Structure
) -> Tuple[float, float]:
    """
    Perform chi-square test comparing two DVHs.

    Parameters
    ----------
    dose1 : Dose
        First (expected) dose distribution.
    dose2 : Dose
        Second (observed) dose distribution.
    structure : Structure
        Structure for DVH computation.

    Returns
    -------
    chi2_statistic : float
        Chi-square test statistic.
    p_value : float
        P-value for the test.

    Notes
    -----
    Tests the null hypothesis that the two DVHs come from the same distribution.
    """
    # Compute DVHs
    dose_bins1, volumes1 = compute_dvh(dose1, structure)
    dose_bins2, volumes2 = compute_dvh(dose2, structure)

    # Create common bins
    max_dose = max(np.max(dose_bins1), np.max(dose_bins2))
    step_size = min(
        dose_bins1[1] - dose_bins1[0] if len(dose_bins1) > 1 else 0.1,
        dose_bins2[1] - dose_bins2[0] if len(dose_bins2) > 1 else 0.1
    )
    common_bins = np.arange(0, max_dose + step_size, step_size)

    # Interpolate
    volumes1_interp = np.interp(common_bins, dose_bins1, volumes1)
    volumes2_interp = np.interp(common_bins, dose_bins2, volumes2)

    # Convert cumulative DVH to differential (histogram)
    diff_volumes1 = -np.diff(np.append(volumes1_interp, 0))
    diff_volumes2 = -np.diff(np.append(volumes2_interp, 0))

    # Ensure non-negative
    diff_volumes1 = np.maximum(diff_volumes1, 0)
    diff_volumes2 = np.maximum(diff_volumes2, 0)

    # Avoid zeros in expected values
    diff_volumes1 = diff_volumes1 + 1e-10

    # Compute chi-square
    chi2_stat, p_value = chisquare(diff_volumes2, diff_volumes1)

    return float(chi2_stat), float(p_value)

compute_dvh_ks_test ¶

compute_dvh_ks_test(dose1: Dose, dose2: Dose, structure: Structure) -> Tuple[float, float]

Perform Kolmogorov-Smirnov test comparing two DVHs.

Parameters¶

dose1 : Dose First dose distribution. dose2 : Dose Second dose distribution. structure : Structure Structure for DVH computation.

Returns¶

ks_statistic : float KS test statistic (maximum difference between CDFs). p_value : float P-value for the test.

Notes¶

The KS test is non-parametric and tests whether two samples come from the same distribution.

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

def compute_dvh_ks_test(
    dose1: Dose,
    dose2: Dose,
    structure: Structure
) -> Tuple[float, float]:
    """
    Perform Kolmogorov-Smirnov test comparing two DVHs.

    Parameters
    ----------
    dose1 : Dose
        First dose distribution.
    dose2 : Dose
        Second dose distribution.
    structure : Structure
        Structure for DVH computation.

    Returns
    -------
    ks_statistic : float
        KS test statistic (maximum difference between CDFs).
    p_value : float
        P-value for the test.

    Notes
    -----
    The KS test is non-parametric and tests whether two samples come from
    the same distribution.
    """
    from scipy.stats import ks_2samp

    # Get dose values in structure for both doses
    dose_values1 = dose1.get_dose_in_structure(structure)
    dose_values2 = dose2.get_dose_in_structure(structure)

    if len(dose_values1) == 0 or len(dose_values2) == 0:
        return np.nan, np.nan

    # Perform two-sample KS test
    ks_stat, p_value = ks_2samp(dose_values1, dose_values2)

    return float(ks_stat), float(p_value)

compute_dvh_confidence_interval ¶

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

Compute DVH confidence intervals from multiple dose distributions.

Useful for uncertainty quantification from multiple treatment plans or Monte Carlo dose simulations.

Parameters¶

doses : list of Dose Multiple dose distributions (e.g., from robust optimization). structure : Structure Structure for DVH computation. confidence : float, optional Confidence level (default: 0.95 for 95% CI).

Returns¶

dose_bins : np.ndarray Dose bin values. mean_dvh : np.ndarray Mean DVH curve. ci_lower : np.ndarray Lower confidence interval. ci_upper : np.ndarray Upper confidence interval.

Examples¶

dose_bins, mean, lower, upper = compute_dvh_confidence_interval( ... [dose1, dose2, dose3], ptv ... ) plt.fill_between(dose_bins, lower, upper, alpha=0.3) plt.plot(dose_bins, mean, 'k-', linewidth=2)

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

def compute_dvh_confidence_interval(
    doses: List[Dose],
    structure: Structure,
    confidence: float = 0.95
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """
    Compute DVH confidence intervals from multiple dose distributions.

    Useful for uncertainty quantification from multiple treatment plans
    or Monte Carlo dose simulations.

    Parameters
    ----------
    doses : list of Dose
        Multiple dose distributions (e.g., from robust optimization).
    structure : Structure
        Structure for DVH computation.
    confidence : float, optional
        Confidence level (default: 0.95 for 95% CI).

    Returns
    -------
    dose_bins : np.ndarray
        Dose bin values.
    mean_dvh : np.ndarray
        Mean DVH curve.
    ci_lower : np.ndarray
        Lower confidence interval.
    ci_upper : np.ndarray
        Upper confidence interval.

    Examples
    --------
    >>> dose_bins, mean, lower, upper = compute_dvh_confidence_interval(
    ...     [dose1, dose2, dose3], ptv
    ... )
    >>> plt.fill_between(dose_bins, lower, upper, alpha=0.3)
    >>> plt.plot(dose_bins, mean, 'k-', linewidth=2)
    """
    if len(doses) == 0:
        raise ValueError("At least one dose is required")

    # Compute DVH for each dose
    all_dvhs = []
    max_dose = 0
    min_step = float('inf')

    for dose in doses:
        dose_bins, volumes = compute_dvh(dose, structure)
        all_dvhs.append((dose_bins, volumes))
        max_dose = max(max_dose, np.max(dose_bins))
        if len(dose_bins) > 1:
            min_step = min(min_step, dose_bins[1] - dose_bins[0])

    if min_step == float('inf'):
        min_step = 0.1

    # Create common dose bins
    common_bins = np.arange(0, max_dose + min_step, min_step)

    # Interpolate all DVHs to common bins
    interpolated_dvhs = []
    for dose_bins, volumes in all_dvhs:
        volumes_interp = np.interp(common_bins, dose_bins, volumes)
        interpolated_dvhs.append(volumes_interp)

    # Stack into array (n_doses x n_bins)
    dvh_array = np.array(interpolated_dvhs)

    # Compute mean and confidence intervals
    mean_dvh = np.mean(dvh_array, axis=0)

    # Compute percentiles for confidence interval
    alpha = 1 - confidence
    lower_percentile = (alpha / 2) * 100
    upper_percentile = (1 - alpha / 2) * 100

    ci_lower = np.percentile(dvh_array, lower_percentile, axis=0)
    ci_upper = np.percentile(dvh_array, upper_percentile, axis=0)

    return common_bins, mean_dvh, ci_lower, ci_upper

compute_dvh_bandwidth ¶

compute_dvh_bandwidth(doses: List[Dose], structure: Structure) -> np.ndarray

Compute DVH bandwidth (maximum difference at each dose level).

Parameters¶

doses : list of Dose Multiple dose distributions. structure : Structure Structure for DVH computation.

Returns¶

bandwidth : np.ndarray Maximum difference between DVHs at each dose bin.

Notes¶

Useful for robust plan evaluation - smaller bandwidth indicates more robust plan.

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

def compute_dvh_bandwidth(
    doses: List[Dose],
    structure: Structure
) -> np.ndarray:
    """
    Compute DVH bandwidth (maximum difference at each dose level).

    Parameters
    ----------
    doses : list of Dose
        Multiple dose distributions.
    structure : Structure
        Structure for DVH computation.

    Returns
    -------
    bandwidth : np.ndarray
        Maximum difference between DVHs at each dose bin.

    Notes
    -----
    Useful for robust plan evaluation - smaller bandwidth indicates
    more robust plan.
    """
    if len(doses) == 0:
        raise ValueError("At least one dose is required")

    # Compute DVH for each dose
    all_dvhs = []
    max_dose = 0
    min_step = float('inf')

    for dose in doses:
        dose_bins, volumes = compute_dvh(dose, structure)
        all_dvhs.append((dose_bins, volumes))
        max_dose = max(max_dose, np.max(dose_bins))
        if len(dose_bins) > 1:
            min_step = min(min_step, dose_bins[1] - dose_bins[0])

    if min_step == float('inf'):
        min_step = 0.1

    # Create common dose bins
    common_bins = np.arange(0, max_dose + min_step, min_step)

    # Interpolate all DVHs to common bins
    interpolated_dvhs = []
    for dose_bins, volumes in all_dvhs:
        volumes_interp = np.interp(common_bins, dose_bins, volumes)
        interpolated_dvhs.append(volumes_interp)

    # Stack into array
    dvh_array = np.array(interpolated_dvhs)

    # Compute bandwidth (max - min at each dose)
    bandwidth = np.max(dvh_array, axis=0) - np.min(dvh_array, axis=0)

    return bandwidth

compute_dvh_similarity_index ¶

compute_dvh_similarity_index(dose1: Dose, dose2: Dose, structure: Structure, method: str = 'dice') -> float

Compute DVH similarity index using various methods.

Parameters¶

dose1 : Dose First dose distribution. dose2 : Dose Second dose distribution. structure : Structure Structure for DVH computation. method : {'dice', 'jaccard', 'correlation', 'cosine'}, optional Similarity metric to use (default: 'dice').

Returns¶

similarity : float Similarity score (0-1, higher is more similar).

Raises¶

ValueError If method is not recognized.

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

def compute_dvh_similarity_index(
    dose1: Dose,
    dose2: Dose,
    structure: Structure,
    method: str = 'dice'
) -> float:
    """
    Compute DVH similarity index using various methods.

    Parameters
    ----------
    dose1 : Dose
        First dose distribution.
    dose2 : Dose
        Second dose distribution.
    structure : Structure
        Structure for DVH computation.
    method : {'dice', 'jaccard', 'correlation', 'cosine'}, optional
        Similarity metric to use (default: 'dice').

    Returns
    -------
    similarity : float
        Similarity score (0-1, higher is more similar).

    Raises
    ------
    ValueError
        If method is not recognized.
    """
    if method not in ['dice', 'jaccard', 'correlation', 'cosine']:
        raise ValueError(f"Unknown method: {method}. Use 'dice', 'jaccard', 'correlation', or 'cosine'.")

    # Compute DVHs
    dose_bins1, volumes1 = compute_dvh(dose1, structure)
    dose_bins2, volumes2 = compute_dvh(dose2, structure)

    # Create common bins and interpolate
    max_dose = max(np.max(dose_bins1), np.max(dose_bins2))
    step_size = min(
        dose_bins1[1] - dose_bins1[0] if len(dose_bins1) > 1 else 0.1,
        dose_bins2[1] - dose_bins2[0] if len(dose_bins2) > 1 else 0.1
    )
    common_bins = np.arange(0, max_dose + step_size, step_size)

    volumes1_interp = np.interp(common_bins, dose_bins1, volumes1)
    volumes2_interp = np.interp(common_bins, dose_bins2, volumes2)

    if method == 'dice':
        # Treat DVH curves as binary masks at each dose level
        intersection = np.minimum(volumes1_interp, volumes2_interp)
        union = volumes1_interp + volumes2_interp
        if np.sum(union) == 0:
            return 0.0
        return float(2.0 * np.sum(intersection) / np.sum(union))

    elif method == 'jaccard':
        # Jaccard index (IoU)
        intersection = np.minimum(volumes1_interp, volumes2_interp)
        union = np.maximum(volumes1_interp, volumes2_interp)
        if np.sum(union) == 0:
            return 0.0
        return float(np.sum(intersection) / np.sum(union))

    elif method == 'correlation':
        # Pearson correlation
        if len(volumes1_interp) < 2:
            return 0.0
        corr = np.corrcoef(volumes1_interp, volumes2_interp)[0, 1]
        return float(corr) if not np.isnan(corr) else 0.0

    elif method == 'cosine':
        # Cosine similarity
        dot_product = np.dot(volumes1_interp, volumes2_interp)
        norm1 = np.linalg.norm(volumes1_interp)
        norm2 = np.linalg.norm(volumes2_interp)
        if norm1 == 0 or norm2 == 0:
            return 0.0
        return float(dot_product / (norm1 * norm2))

    return 0.0

Usage Examples¶

Computing a Basic DVH¶

from dosemetrics import Dose, Structure
from dosemetrics.metrics.dvh import compute_dvh

dose = Dose.from_nifti("dose.nii.gz")
ptv  = Structure.from_nifti("ptv.nii.gz", name="PTV")

dose_bins, volumes = compute_dvh(dose, ptv)

Conformity and Homogeneity Metrics¶

from dosemetrics.metrics.conformity import (
    compute_conformity_index,          # ICRU CI = V_target_rx / V_rx
    compute_rtog_conformity_index,     # RTOG CI = V_rx / V_target
    compute_paddick_conformity_index,  # Paddick / van't Riet CI = (V_target_rx)² / (V_target × V_rx)
    compute_conformity_number,         # van't Riet CN — same formula as Paddick CI
    compute_coverage,                  # V_target_rx / V_target
    compute_spillage,                  # (V_rx - V_target_rx) / V_rx
    compute_prescription_mae,          # mean |dose - prescription| in target
)
from dosemetrics.metrics.homogeneity import (
    compute_homogeneity_index,   # ICRU HI = (D2 - D98) / D50
    compute_gradient_index,      # Paddick-Lippitz GI = V_50% / V_100%
    compute_dose_homogeneity,    # coefficient of variation
    compute_uniformity_index,    # UI = 1 - (Dmax - Dmin) / Dref
)

prescription = 60.0  # Gy

ci_icru  = compute_conformity_index(dose, ptv, prescription)
ci_rtog  = compute_rtog_conformity_index(dose, ptv, prescription)
ci_pad   = compute_paddick_conformity_index(dose, ptv, prescription)
hi       = compute_homogeneity_index(dose, ptv)
gi       = compute_gradient_index(dose, ptv, prescription)
rx_mae   = compute_prescription_mae(dose, ptv, prescription)

print(f"ICRU CI:        {ci_icru:.3f}")
print(f"RTOG CI:        {ci_rtog:.3f}")
print(f"Paddick CI:     {ci_pad:.3f}")
print(f"HI (ICRU 83):   {hi:.3f}")
print(f"Gradient Index: {gi:.2f}")
print(f"Prescription MAE: {rx_mae:.2f} Gy")

DVH Comparison Metrics¶

from dosemetrics.metrics.dvh import compute_dvh_score, compute_dvh_auc
from dosemetrics.metrics.dose_comparison import (
    compute_normalized_mae,
    compute_variance_of_laplacian,
)

# DVH Score: average |D1|, |D95|, |D99| difference (Gy)
score = compute_dvh_score(dose_reference, dose_evaluated, ptv)

# DVH AUC: integral of DVH curve, normalised to [0, 1]
auc = compute_dvh_auc(dose, ptv, normalize=True)

# Normalized MAE with high-dose masking
n_mae = compute_normalized_mae(
    dose_reference,
    dose_evaluated,
    normalization_value=60.0,
    dose_threshold_gy=5.0,
)

# Dose sharpness (Variance of Laplacian)
vol = compute_variance_of_laplacian(dose)

print(f"DVH Score:      {score:.2f} Gy")
print(f"DVH AUC:        {auc:.3f}")
print(f"Normalized MAE: {n_mae:.4f}")
print(f"VoL (sharpness): {vol:.4f}")

Dose Statistics¶

from dosemetrics.metrics.dvh import compute_dose_statistics

stats = compute_dose_statistics(dose, ptv)
print(f"Mean dose: {stats['mean_dose']:.2f} Gy")
print(f"D95:       {stats['D95']:.2f} Gy")
print(f"D2 (near-max): {stats['D02']:.2f} Gy")