"""
PyPlanAnalysis.plan
==================
Main user-facing API: PatientPlan (loads DICOM data and runs the
analysis) and AnalysisResults (metrics, plots, NTCP, and exports).
See the project README for usage examples.
"""
import warnings
import numpy as np
import pandas as pd
from pathlib import Path
import re
import pydicom
from typing import Union
import matplotlib.pyplot as plt
from .io import (load_ct_series,load_dose_grid, get_grid_geometry, get_fractional_mask_on_grid,get_structure_mask_on_grid,
get_all_structure_names, find_dicom_files,resample_dose_on_ct,_np_to_sitk,resample_dose_to_new_grid,get_roi_center_of_mass)
from .rbe import RBEConfig, rbe_fixed, compute_rbe_dose
from .metrics import (RadiobiologyConfig, MetricConfig,
compute_dvh_metrics, compute_let_metrics,
compute_cumulative_histogram, compute_2d_histogram)
from .plots import (plot_cumulative_histogram, plot_dvh_comparison,
plot_dlvh_2d, plot_dlvh_3d, save_figures_to_pdf, save_figure)
from .NTCP import (NTCPConfig,NTCPModelBase)
import copy
import SimpleITK as sitk
# ============================================================
# PatientPlan
# ============================================================
[docs]
class PatientPlan:
"""
Represents a single patient's proton therapy plan.
Parameters
----------
patient_id : str
Identifier used in output filenames and CSV column.
dose_file : str or Path
DICOM RT Dose file containing physical dose [Gy].
let_file : str or Path
DICOM RT Dose file containing LETd [keV/µm] stored as dose grid.
rtstruct : str or Path
DICOM RT Struct file.
CT : str or Path
folder of DICOM CT files
Class method
------------
PatientPlan.from_folder(folder, patient_id=None)
Auto-discover DICOM files from a folder.
"""
def __init__(self,
patient_id : str | None = None,
plan_file : Union[str, Path, None] = None,
dose_file : Union[str, Path, None] = None,
let_file : Union[str, Path, None] = None,
rtstruct : Union[str, Path, None] = None,
CT_folder_path : Union[str, Path, None] = None,
n_fractions : int | None = None):
self.plan_file = Path(plan_file) if plan_file else None
self.dose_file = Path(dose_file) if dose_file else None
self.let_file = Path(let_file) if let_file else None
self.rtstruct = Path(rtstruct) if rtstruct else None
self.CT_folder = Path(CT_folder_path) if CT_folder_path else None
self.n_fractions = n_fractions
self._is_scaling_performed_from_RBEw_dose = False
if patient_id is not None:
self.patient_id = patient_id
elif dose_file is not None:
ds = pydicom.dcmread(str(dose_file), stop_before_pixels=True, force=True)
self.patient_id = ds.PatientID
elif rtstruct is not None:
ds = pydicom.dcmread(str(rtstruct), stop_before_pixels=True, force=True)
self.patient_id = ds.PatientID
else:
print("can´t resolve patient name, setting it to UNKNOWN")
self.patient_id = "UNKNOWN"
self._plan_ds = None
self._dose_arr = None
self._dose_ds = None
self._let_arr = None
self._let_ds = None
self._rtstruct_ds = None
self._CT_sitk = None
self._CT_geom = None
self._loaded = False
self._CT_loaded = False
# ----------------------------------------------------------
[docs]
@classmethod
def from_folder(cls,
folder: Union[str, Path],
n_fractions: int | None = None,
rad_type: str = "PROTON") -> "PatientPlan":
"""
Create a PatientPlan by auto-discovering DICOM files in folder.
Parameters
----------
folder : path containing .dcm files
patient_id : if None, uses the folder name
"""
folder = Path(folder)
files = find_dicom_files(folder, rad_type)
missing = [k for k, v in files.items() if v is None]
if missing:
if 'rtplan' in missing:
print(
f"⚠ Could not find {missing} in {folder}."
)
else:
print(
f"⚠ Could not find {missing} in {folder}. "
"Set file paths manually via PatientPlan(...)."
)
found = [k for k, v in files.items() if v is not None]
if found:
pat = files["Patient_ID"]
print(f"\nSearched for patient {pat}")
for f in found:
if f not in ["linked","link_warnings","Patient_ID"]:
print(f"✓ Ready to load {f}: {files[f]}")
return cls(
patient_id = files["Patient_ID"],
plan_file = files["rtplan"],
dose_file = files["dose"],
let_file = files["let"],
rtstruct = files["rtstruct"],
CT_folder_path= files["CT"],
n_fractions = n_fractions
)
# ----------------------------------------------------------
[docs]
def load(self):
"""Load and cache all DICOM data. Called automatically by analyse()."""
if self._loaded:
return
print("\n")
print(f"Loading DICOM files for {self.patient_id}...")
if self.plan_file is not None:
self._plan_ds = pydicom.dcmread(str(self.plan_file), force=True)
n_fractions_from_rtplan = self._plan_ds.FractionGroupSequence[0].NumberOfFractionsPlanned
print(f"Selected plan: {self._plan_ds.RTPlanName}")
try:
radtype= self._plan_ds.RadiationType
print(f"Radiation Type: {radtype}")
except:
print("Radiation Type: unknown")
if self.n_fractions is not None:
if n_fractions_from_rtplan != self.n_fractions:
print(f'⚠ Number of fractions specified does not match RTPLAN, overwriting with RTPLAN value: {n_fractions_from_rtplan} fractions')
self.n_fractions = n_fractions_from_rtplan
else:
self.n_fractions = n_fractions_from_rtplan
if self.plan_file is None and self.n_fractions is None:
raise AttributeError(
"No fraction number is given, please define it or select RTPLAN"
)
if self.dose_file is not None:
#if available dose is effective, scale down by 10%
self._dose_arr, self._dose_ds = load_dose_grid(self.dose_file)
dose_type = getattr(self._dose_ds, "DoseType", "").upper()
if dose_type == "EFFECTIVE":
self._is_scaling_performed_from_RBEw_dose = True
self._dose_arr = self._dose_arr/1.1
if self.let_file is not None:
self._let_arr, self._let_ds = load_dose_grid(self.let_file)
if self.rtstruct is not None:
self._rtstruct_ds = pydicom.dcmread(str(self.rtstruct), force=True)
if self.CT_folder is not None:
try:
self._CT_sitk, self._CT_geom = load_ct_series(str(self.CT_folder))
self._CT_loaded = True
except Exception as e:
warnings.warn(f"Could not load CT series: {e}")
self._CT_loaded = False
else:
self._CT_loaded = False
self._loaded = True
@property
def structure_names(self) -> list:
"""List of all structure names in the RT Struct."""
if not self._loaded:
self.load()
return get_all_structure_names(self._rtstruct_ds)
#Not used
@property
def voxel_volume_cc(self) -> float:
"""Volume of one dose voxel in cc."""
if not self._loaded:
self.load()
if self._CT_loaded:
spacing = self._CT_geom['spacing']
elif self._dose_ds is not None:
_, spacing = get_grid_geometry(self._dose_ds)
elif self._let_ds is not None:
_, spacing = get_grid_geometry(self._let_ds)
else:
raise ValueError(
"No grid geometry available — need at least one of "
"dose, LET, or CT to determine voxel volume."
)
return float(np.prod(spacing)) / 1000.0
# ----------------------------------------------------------
def _sanity_check_and_resample_dose(self,
structures : list = None,
metric_cfg : MetricConfig = None,
resample_on_CT : bool = False,
resample_on_custom_grid : bool = False):
self.load()
# ---------------------------------------------------------
# Validate loaded data
# ---------------------------------------------------------
if self._dose_ds is not None:
assert tuple(self._dose_arr.shape) == (
int(self._dose_ds.NumberOfFrames),
int(self._dose_ds.Rows),
int(self._dose_ds.Columns),
)
if self._dose_arr is None and self._let_arr is None:
warnings.warn("Neither dose nor LET could be loaded — nothing to analyse.")
return None
if self._rtstruct_ds is None:
warnings.warn("No RT Struct loaded — cannot extract any structures.")
return None
# ---------------------------------------------------------
# Local copies: never modify self's pristine data
# ---------------------------------------------------------
# Using local variables makes analyse() idempotent with respect to
# self: repeated calls, in any order, with any combination of flags,
# always start from the same pristine self._dose_arr/self._let_arr.
dose_arr = self._dose_arr.copy() if self._dose_arr is not None else None
dose_ds = copy.deepcopy(self._dose_ds) if self._dose_ds is not None else None
let_arr = self._let_arr.copy() if self._let_arr is not None else None
let_ds = copy.deepcopy(self._let_ds) if self._let_ds is not None else None
# --- resample dose/LET only if the array exists ---
sitk_dose = _np_to_sitk(dose_arr, dose_ds) if dose_arr is not None else None
sitk_let = _np_to_sitk(let_arr, let_ds) if let_arr is not None else None
grid_origin = grid_spacing = grid_shape = grid_z_positions = None
if resample_on_CT and not self._CT_loaded :
raise ValueError("No CT has been loaded for resampling")
if resample_on_custom_grid and resample_on_CT:
raise ValueError("Two resampling criteria have been selected. Please select only one.")
if self._CT_loaded and resample_on_CT:
print(" Resampling dose and LET onto CT grid...")
if sitk_dose is not None:
res_dose_sitk = resample_dose_on_ct(sitk_dose,self._CT_sitk)
dose_arr = sitk.GetArrayFromImage(res_dose_sitk) #sitk_dose
if sitk_let is not None:
res_let_sitk = resample_dose_on_ct(sitk_let, self._CT_sitk)
let_arr = sitk.GetArrayFromImage(res_let_sitk) #sitk_let
#Set grid for resampling and creating struct masks
grid_origin = self._CT_geom['origin']
grid_spacing = self._CT_geom['spacing'] #x,y,z
grid_shape = self._CT_geom['shape']
grid_z_positions = self._CT_geom['z_positions']
elif resample_on_custom_grid:
grid_new = metric_cfg.New_grid
print(f" Resampling dose and LET onto custom {grid_new} grid...")
if sitk_dose is not None:
res_dose_sitk, dose_arr, dose_geom, _ = resample_dose_to_new_grid(sitk_dose,
dose_ds,
grid_new)
if sitk_let is not None:
res_let_sitk, let_arr, _, let_ds = resample_dose_to_new_grid (sitk_let,
let_ds,
grid_new)
#Set grid for resampling and creating struct masks
grid_origin = dose_geom['origin']
grid_spacing = dose_geom['spacing'] #x,y,z
grid_shape = dose_geom['shape'] #z,y,x
grid_z_positions = dose_geom['z_positions']
else:
print("No Resampling for dose and LET")
# no CT, no resample requested — use whichever grid (dose preferred, else LET) is available
ref_ds = dose_ds if dose_ds is not None else let_ds
if ref_ds is not None:
grid_origin, grid_spacing = get_grid_geometry(ref_ds)
grid_shape = (dose_arr.shape if dose_arr is not None
else let_arr.shape)
z_offsets = [float(v) for v in ref_ds.GridFrameOffsetVector]
grid_z_positions = np.array([grid_origin[2] + o for o in z_offsets])
if grid_origin is None:
warnings.warn("Could not determine a voxel grid — aborting analysis.")
return None
return dose_arr, let_arr, grid_origin, grid_spacing, grid_shape, grid_z_positions
[docs]
def get_roi_mask(self,struct_name, rtstruct_ds, grid_origin, grid_spacing,
grid_shape, grid_z_positions, use_fractional=False,
supersample=4, max_supersample=12):
if use_fractional:
frac = get_fractional_mask_on_grid(
struct_name, rtstruct_ds,
grid_origin, grid_spacing,
grid_shape, grid_z_positions,
supersample=supersample,
max_supersample=max_supersample
)
mask = frac > 0
return frac[mask], mask
return None, get_structure_mask_on_grid(
struct_name, rtstruct_ds,
grid_origin, grid_spacing,
grid_shape, grid_z_positions
)
def _get_priority(self, struct_name: str, priority_map: dict) -> int:
"""Substring match on structure name, same as alpha_beta lookup."""
name = struct_name.lower()
for key, val in priority_map.items():
if key in name:
return val
return 999 # lowest priority when not listed
[docs]
def build_alpha_beta_map(self, rtstruct_ds, structures, grid_origin, grid_spacing,
grid_shape, grid_z_positions, radiobio_cfg,
use_fractional,supersample, max_supersample):
ab_map = np.full(grid_shape, radiobio_cfg.alpha_beta_default, dtype=np.float32)
priority_map = np.full(grid_shape, 1000, dtype=np.int32)
for struct_name in structures:
try:
_,mask = self.get_roi_mask(struct_name, rtstruct_ds, grid_origin, grid_spacing,
grid_shape, grid_z_positions, use_fractional,
supersample, max_supersample
)
except Exception:
continue
if not mask.any():
continue
priority = self._get_priority(
struct_name, radiobio_cfg.structure_priority
)
update = mask & (priority < priority_map)
ab_map[update] = radiobio_cfg.get_alpha_beta(struct_name)
priority_map[update] = priority
return ab_map, priority_map
[docs]
def compute_vRBE(self,
model : str,
structures : list = None,
metric_cfg : MetricConfig = None,
radiobio_cfg : RadiobiologyConfig = None,
resample_on_CT : bool = False,
resample_on_custom_grid : bool = False,
use_fractional: bool = False,
supersample : int = 4,
max_supersample: int = 12):
"""
Compute a voxel-wise variable-RBE dose volume for the whole patient grid.
It Works independently of analyse() — no DVH metrics are computed, only the 3D RBE-weighted dose map and the α/β map used to produce it.
It works independetly of RBEConfig
Parameters
----------
model : str — RBE model to apply: "linear", "mcnamara", "wedenberg", or "carabe"
structures : list of structure name strings used to build the α/β map, or None (= all).
Voxels not inside any listed structure receive the default α/β.
metric_cfg : MetricConfig — only New_grid is used, when resample_on_custom_grid=True
radiobio_cfg : RadiobiologyConfig — supplies per-structure α/β values and overlap
priority for building the voxel-wise α/β map
resample_on_CT : if True, resample dose and LET onto the CT grid before computing RBE.
Requires CT to have been loaded (ct_dir passed to PatientPlan)
resample_on_custom_grid : if True, resample dose and LET onto the isotropic grid defined
by metric_cfg.New_grid [mm]. Mutually exclusive with resample_on_CT.
use_fractional : if True, use fractional voxel membership for the α/β map boundary voxels.
More accurate for small structures. Default: False.
supersample : baseline in-plane supersampling factor (N×N sub-points per voxel).
Only used when use_fractional=True. Scaled up automatically for small
structures — see io._adaptive_supersample. Default: 4.
max_supersample : hard cap on the adaptive in-plane supersampling factor, to bound
compute cost for very small ROIs. Default: 12.
Returns
-------
vRBEResults
A container holding the 3D vRBE dose array, the α/β map, grid geometry,
and a to_nii() method to save both volumes as NIfTI files.
"""
metric_cfg = metric_cfg or MetricConfig()
radiobio_cfg = radiobio_cfg or RadiobiologyConfig()
if structures is None:
structures = self.structure_names
else:
available = self.structure_names
missing = [s for s in structures if s not in available]
if missing:
warnings.warn(f"Structures not in RT Struct: {missing}")
structures = [s for s in structures if s in available]
mode = f"fractional (supersample={supersample})" if use_fractional else "binary"
print(f"Analysing {len(structures)} structures [{mode} mask]...")
#--------------------
# Perform resampling and sanity checks
#-----------------------------
dose_arr, let_arr, grid_origin, grid_spacing, grid_shape, grid_z_positions = self._sanity_check_and_resample_dose( structures,
metric_cfg,
resample_on_CT,
resample_on_custom_grid)
ab_map,_ = self.build_alpha_beta_map( self._rtstruct_ds, structures, grid_origin, grid_spacing,
grid_shape, grid_z_positions, radiobio_cfg,
use_fractional,supersample,max_supersample)
vRBE = compute_rbe_dose(dose_arr,let_arr, self.n_fractions, ab_map, model)
return vRBEResults(patient_id = self.patient_id,
model = model.lower(),
dose_arr = dose_arr,
let_arr = let_arr,
vRBE = vRBE,
ab_map = ab_map,
geom = {"origin" : grid_origin,"spacing" : grid_spacing} )
[docs]
def analyse(self,
structures : list = None,
rbe_cfg : RBEConfig = None,
metric_cfg : MetricConfig = None,
radiobio_cfg : RadiobiologyConfig = None,
resample_on_CT : bool = False,
resample_on_custom_grid : bool = False,
use_fractional: bool = False,
supersample : int = 4,
max_supersample: int = 12) -> "AnalysisResults":
"""
Run the full analysis pipeline for this patient.
Parameters
----------
structures : list of structure name strings, or None (= all)
rbe_cfg : RBEConfig (default: all three models, RBE=1.1)
metric_cfg : MetricConfig (default values)
radiobio_cfg : RadiobiologyConfig (default per-tissue α/β and gEUD-a)
resample: if True, interpolate dose values to the new spacing defined in MetricConfig
use_fractional: if True, use fractional voxel membership mask instead
of binary mask. More accurate for small structures.
Each voxel is weighted by the fraction of its volume
inside the contour. Default: False.
supersample : baseline in-plane supersampling factor (N x N sub-points
per voxel). Only used when use_fractional=True.
Scaled up automatically for small structures — see
io._adaptive_supersample. Default: 4.
max_supersample: hard cap on the adaptive in-plane supersampling
factor, to bound compute cost for very small ROIs.
Returns
-------
AnalysisResults
"""
rbe_cfg = rbe_cfg or RBEConfig()
metric_cfg = metric_cfg or MetricConfig()
radiobio_cfg = radiobio_cfg or RadiobiologyConfig()
if structures is None:
structures = self.structure_names
else:
available = self.structure_names
missing = [s for s in structures if s not in available]
if missing:
warnings.warn(f"Structures not in RT Struct: {missing}")
structures = [s for s in structures if s in available]
mode = f"fractional (supersample={supersample})" if use_fractional else "binary"
print(f"Analysing {len(structures)} structures [{mode} mask]...")
#--------------------
# Perform resampling and sanity checks
#-----------------------------
dose_arr, let_arr, grid_origin, grid_spacing, grid_shape, grid_z_positions = self._sanity_check_and_resample_dose( structures,
metric_cfg,
resample_on_CT,
resample_on_custom_grid)
all_rows = []
dvh_data = {} # {struct: {dose_type: (edges, cum)}}
lvh_data = {} # {struct: (edges, cum)}
dlvh_data_diff = {} # {struct: (H, d_edges, l_edges)}
dlvh_data_cum = {} # {struct: (H, d_edges, l_edges)}
# Voxel volume MUST come from the grid actually used to build the
# mask (grid_spacing, set above per the branch actually taken), not
# from self.voxel_volume_cc.
vol_cc = float(np.prod(grid_spacing)) / 1000.0
for struct_name in structures:
print(f" → {struct_name}")
try:
weights,mask = self.get_roi_mask(struct_name,
self._rtstruct_ds,
grid_origin,
grid_spacing,
grid_shape,
grid_z_positions,
use_fractional,
supersample = supersample,
max_supersample = max_supersample)
except Exception as e:
warnings.warn(f"Skipping '{struct_name}': {e}")
continue
if mask.sum() == 0:
warnings.warn(f"Empty mask for '{struct_name}' – skipping.")
continue
ab = radiobio_cfg.get_alpha_beta(struct_name)
geud_a = radiobio_cfg.get_geud_a(struct_name)
dose_vox = dose_arr[mask] if dose_arr is not None else None
let_vox = let_arr[mask] if let_arr is not None else None
# effective volume: sum of fractional weights × voxel volume
if weights is not None:
eff_vol_cc = float(weights.sum() * vol_cc)
else:
eff_vol_cc = float(mask.sum() * vol_cc)
try:
com = get_roi_center_of_mass(struct_name, self._rtstruct_ds)
except Exception as e:
warnings.warn(f"Could not compute center of mass for '{struct_name}': {e}")
com = None
row = {
"patient_id" : self.patient_id,
"ROI_Name" : struct_name,
"volume_cc" : eff_vol_cc,
"alpha_beta" : ab,
"geud_a" : geud_a,
"COM_x" : com[0] if com is not None else np.nan,
"COM_y" : com[1] if com is not None else np.nan,
"COM_z" : com[2] if com is not None else np.nan,
}
dvh_data[struct_name] = {}
# ---- Physical dose ----
if dose_arr is not None:
row.update(compute_dvh_metrics(dose_vox, vol_cc, "Phys",
metric_cfg, ab, geud_a,
weights=weights,n_fractions = self.n_fractions))
edges, cum = compute_cumulative_histogram(dose_vox,
metric_cfg.dvh_bins,
weights=weights)
dvh_data[struct_name]["Phys"] = (edges, cum)
# ---- Fixed RBE ----
dose_fixed = rbe_fixed(dose_vox, rbe_cfg.fixed_rbe)
lbl_fixed = f"RBE{rbe_cfg.fixed_rbe}"
row.update(compute_dvh_metrics(dose_fixed, vol_cc, lbl_fixed ,
metric_cfg, ab, geud_a,
weights=weights,n_fractions = self.n_fractions))
edges_f, cum_f = compute_cumulative_histogram(dose_fixed,
metric_cfg.dvh_bins,
weights=weights)
dvh_data[struct_name][lbl_fixed] = (edges_f, cum_f)
if let_arr is not None:
# ---- Variable RBE (per model) ----
for model in rbe_cfg.models:
rbe_dose_vox = compute_rbe_dose(dose_vox, let_vox, self.n_fractions, ab, model)
row.update(compute_dvh_metrics(rbe_dose_vox, vol_cc, model,
metric_cfg, ab, geud_a,
weights=weights,n_fractions = self.n_fractions))
e_m, c_m = compute_cumulative_histogram(rbe_dose_vox,
metric_cfg.dvh_bins,
weights=weights)
dvh_data[struct_name][model] = (e_m, c_m)
# ---- LET-thresholded DVH metrics (high-LET subvolume on fixed RBE - dirty dose) ----
high_let_idx = let_vox >= metric_cfg.LET_thr
if high_let_idx.sum() > 0:
w_high = weights[high_let_idx] if weights is not None else None
row.update(compute_dvh_metrics(
dose_vox[high_let_idx], vol_cc,
"highLET_fixed_rbe", metric_cfg, ab, geud_a,
weights=w_high,n_fractions = self.n_fractions))
# ---- 2-D DLVH ----
H_diff, H_cum, d_edges, l_edges_2d = compute_2d_histogram(
dose_fixed, let_vox,
metric_cfg.dlvh_dose_bins,
metric_cfg.dlvh_let_bins,
weights=weights)
dlvh_data_diff[struct_name] = (H_diff, d_edges, l_edges_2d)
dlvh_data_cum[struct_name] = (H_cum, d_edges, l_edges_2d)
if let_arr is not None:
# ---- LET metrics & LVH ----
row.update(compute_let_metrics(let_vox, dose_vox, metric_cfg.lx,
weights=weights))
l_edges, l_cum = compute_cumulative_histogram(let_vox,
metric_cfg.lvh_bins,
weights=weights)
lvh_data[struct_name] = (l_edges, l_cum)
all_rows.append(row)
metrics_df = pd.DataFrame(all_rows)
return AnalysisResults(
patient_id = self.patient_id,
metrics_df = metrics_df,
dvh_data = dvh_data,
lvh_data = lvh_data,
dlvh_data_diff = dlvh_data_diff,
dlvh_data_cum = dlvh_data_cum,
)
# ============================================================
# AnalysisResults
# ============================================================
[docs]
class AnalysisResults:
"""
Container for all outputs from PatientPlan.analyse().
Attributes
----------
patient_id : str
metrics_df : pd.DataFrame – all DVH/LVH metrics, one row per structure
dvh_data : dict {struct: {dose_type: (edges, cum)}} e.g. # {"CTV": {"Physical": (edges, cum), "RBE×1.1": (edges, cum), ...}, "Brainstem": {"Physical": (edges, cum), ...},
lvh_data : dict {struct: (edges, cum)}
dlvh_data : dict {struct: (H, dose_edges, let_edges)}
"""
def __init__(self,
patient_id : str,
metrics_df : pd.DataFrame,
dvh_data : dict,
lvh_data : dict,
dlvh_data_diff : dict,
dlvh_data_cum : dict):
self.patient_id = patient_id
self.metrics_df = metrics_df
self.dvh_data = dvh_data
self.lvh_data = lvh_data
self.dlvh_data_diff = dlvh_data_diff
self.dlvh_data_cum = dlvh_data_cum
# ----------------------------------
# NTCP calculation
# ----------------------------------
[docs]
def CalcNTCP (self, NTCPConfig: NTCPConfig, path_models_xls: str = None):
# load list with model names and corresponding parameter names
list_RBE_calculated = list(self.dvh_data[list(self.dvh_data.keys())[0]].keys())
NTCP_results = {}
for NTCP_model in NTCPConfig.models:
roi_override = NTCPConfig.roi_overrides.get(NTCP_model) if NTCPConfig.roi_overrides else None
NTCP_base = NTCPModelBase(NTCP_model, path_models_xls,
roi_name = roi_override,
ctv_name = NTCPConfig.ctv_name)
NTCP_dict = {}
for RBE in list_RBE_calculated:
NTCP_dict[RBE] = NTCP_base.compute_NTCP(RBE, self.metrics_df)
NTCP_results[NTCP_model] = NTCP_dict
return NTCP_results
# ----------------------------------------------------------
# Export methods
# ----------------------------------------------------------
[docs]
def to_csv(self, path: Union[str, Path] = None) -> Path:
"""
Save metrics to CSV.
Parameters
----------
path : output path. Default: "<patient_id>/dvh_metrics.csv"
"""
path = Path(path or f"{self.patient_id}/dvh_metrics.csv")
path.parent.mkdir(parents=True, exist_ok=True)
self.metrics_df.to_csv(path, index=False)
print(f" Saved CSV: {path}")
return path
[docs]
def to_excel(self, path: Union[str, Path] = None) -> Path:
"""
Save metrics to Excel (.xlsx).
Parameters
----------
path : output path. Default: "<patient_id>/dvh_metrics.xlsx"
"""
path = Path(path or f"{self.patient_id}/dvh_metrics.xlsx")
path.parent.mkdir(parents=True, exist_ok=True)
self.metrics_df.to_excel(path, index=False, sheet_name=re.sub(r'[/*?\[\]<>|:]', '',self.patient_id))
print(f" Saved Excel: {path}")
return path
# ----------------------------------------------------------
# Plot methods
# ----------------------------------------------------------
[docs]
def plot_dvh(self,
save : Union[str, Path] = None,
dpi : int = 350) -> list:
"""
Generate DVH curve figures.
Creates one figure per dose type (Physical, RBE×1.1, each model)
showing all structures, PLUS one per-structure comparison figure.
Parameters
----------
save : path for output PDF. Default: "<patient_id>/dvh_curves.pdf"
dpi : figure DPI
Returns
-------
list of matplotlib Figures
"""
save = Path(save or f"{self.patient_id}/dvh_curves.pdf")
save.parent.mkdir(parents=True, exist_ok=True)
figs = []
# One figure per dose label (all structures overlaid)
if self.dvh_data:
all_labels = list(next(iter(self.dvh_data.values())).keys())
for dose_type in all_labels:
curves = {struct: data[dose_type]
for struct, data in self.dvh_data.items()
if dose_type in data}
fig = plot_cumulative_histogram(
curves,
title = f"{self.patient_id} – DVH ({dose_type})",
xlabel = "Dose [Gy / Gy(RBE)]",
dpi = dpi)
figs.append(fig)
# Per-structure comparison (all dose types on one plot)
for struct, label_curves in self.dvh_data.items():
fig = plot_dvh_comparison(label_curves, struct, dpi=dpi)
figs.append(fig)
save_figures_to_pdf(figs, save, dpi=dpi)
print(f" Saved DVH PDF: {save}")
return figs
[docs]
def plot_lvh(self,
save : Union[str, Path] = None,
dpi : int = 150) -> plt.Figure:
"""
Generate LVH curve figure (all structures on one plot).
Parameters
----------
save : path for output PDF. Default: "<patient_id>/lvh_curves.pdf"
Returns
-------
matplotlib Figure
"""
save = Path(save or f"{self.patient_id}/lvh_curves.pdf")
save.parent.mkdir(parents=True, exist_ok=True)
fig = plot_cumulative_histogram(
self.lvh_data,
title = f"{self.patient_id} – LVH",
xlabel = "LETd [keV/μm]",
ylabel = "Volume fraction",
dpi = dpi)
save_figures_to_pdf([fig], save, dpi=dpi)
print(f" Saved LVH PDF: {save}")
return fig
[docs]
def plot_dlvh(self,
save : Union[str, Path] = None,
dpi : int = 150) -> list:
"""
Generate 2-D DLVH figures – one PNG per structure.
Parameters
----------
save : output directory. Default: "<patient_id>/dlvh_2d/"
Returns
-------
list of (structure_name, Path) tuples
"""
save_dir = Path(save or f"{self.patient_id}/dlvh_2d")
save_dir.mkdir(parents=True, exist_ok=True)
saved = []
for struct, (H, d_edges, l_edges) in self.dlvh_data_diff.items():
fig = plot_dlvh_2d(H, d_edges, l_edges, struct, dpi=dpi)
name = struct.replace(" ", "_").replace("/", "-")
out = save_dir / f"{name}_Diff.png"
save_figure(fig, out, dpi=dpi)
saved.append((struct, out))
saved = []
for struct, (H, d_edges, l_edges) in self.dlvh_data_cum.items():
fig = plot_dlvh_2d(H, d_edges, l_edges, struct, dpi=dpi)
name = struct.replace(" ", "_").replace("/", "-")
out = save_dir / f"{name}_Cumulative.png"
save_figure(fig, out, dpi=dpi)
saved.append((struct, out))
for struct, (H, d_edges, l_edges) in self.dlvh_data_cum.items():
fig = plot_dlvh_3d(H, d_edges, l_edges, struct, dpi=dpi)
name = struct.replace(" ", "_").replace("/", "-")
out = save_dir / f"{name}_Cumulative3D.png"
save_figure(fig, out, dpi=dpi)
saved.append((struct, out))
#print(f" Saved {len(saved)} DLVH plots → {save_dir}/")
return saved
[docs]
def save_all(self,
output_dir : Union[str, Path] = None,
csv_name : str = "dvh_metrics.csv",
excel_name : str = "dvh_metrics.xlsx",
dpi : int = 150):
"""
Convenience method: save CSV, Excel, DVH PDF, LVH PDF, DLVH PNGs.
Parameters
----------
output_dir : base output directory. Default: patient_id/
csv_name : filename for the CSV output
excel_name : filename for the Excel output
dpi : figure DPI for all plots
"""
base = Path(output_dir or self.patient_id)
base.mkdir(parents=True, exist_ok=True)
self.to_csv( base / csv_name)
self.to_excel( base / excel_name)
self.plot_dvh( base / "dvh_curves.pdf", dpi=dpi)
self.plot_lvh( base / "lvh_curves.pdf", dpi=dpi)
self.plot_dlvh(base / "dlvh_2d", dpi=dpi)
# ============================================================
# vRBEResults
# ============================================================
[docs]
class vRBEResults:
"""
Container for all outputs from PatientPlan.compute_vRBE().
Attributes
----------
patient_id : str
model : str
dose_arr : np.array(), of 1.1-weighted dose
let_arr : np.array(), of Let in kev/um
vRBE : np.array (float) of variablRBE-weighted dose
ab_map : np.array (int)?
geom : dict (órigin,spacing) in (x,y,z). Reminder: SimpleITK uses x-y-z. Numpy uses z-y-x.
"""
def __init__(self,
patient_id : str,
model : str,
dose_arr : np.ndarray,
let_arr : np.ndarray,
vRBE : np.ndarray,
ab_map : np.ndarray,
geom : dict ):
self.patient_id = patient_id
self.model = model
self.dose_arr = dose_arr
self.let_arr = let_arr
self.vRBE = vRBE
self.ab_map = ab_map if not model.lower() == "linear" else None
self.geom = geom
# ----------------------------------------------------------
# Export methods
# ----------------------------------------------------------
[docs]
def to_nii(self,
path : Union[str, Path],
name_vRBE : str = "vRBE.nii.gz",
name_ab : str = "ab_map.nii.gz",
name_cRBE : str = "cRBE.nii.gz",
name_let : str = "LETd.nii.gz") -> None:
"""
Save vRBE dose and α/β map as NIfTI files.
Parameters
----------
path : output directory (created if it does not exist)
name_vrbe : filename for the vRBE dose volume
name_ab : filename for the α/β map volume
"""
out_dir = Path(path)
out_dir.mkdir(parents=True, exist_ok=True)
def _arr_to_sitk(arr, origin, spacing):
img = sitk.GetImageFromArray(arr.astype(np.float32))
img.SetOrigin([float(v) for v in origin])
img.SetSpacing([float(v) for v in spacing])
return img
origin = self.geom["origin"] # adjust if geom is a dataclass: self.geom.origin
spacing = self.geom["spacing"]
sitk.WriteImage(_arr_to_sitk(self.vRBE, origin, spacing),
str(out_dir / name_vRBE))
sitk.WriteImage(_arr_to_sitk(self.dose_arr, origin, spacing),
str(out_dir / name_cRBE))
sitk.WriteImage(_arr_to_sitk(self.let_arr, origin, spacing),
str(out_dir / name_let))
if self.ab_map is not None:
sitk.WriteImage(_arr_to_sitk(self.ab_map, origin, spacing),
str(out_dir / name_ab))
print(f" ✓ Saved {name_cRBE},{name_vRBE},{name_let} and {name_ab} → {out_dir}")