PyPlanAnalysis.io

DICOM auto-discovery and loading, CT/dose grid resampling, and RT Struct mask extraction (binary and fractional/supersampled).

PyPlanAnalysis.io

DICOM discovery and loading, CT/dose grid resampling, and RT Struct mask extraction (binary and fractional).

find_dicom_files(folder, rad_type)

Auto-discover RT Dose (physical dose, LET), RT Struct, RT Plan and CT files belonging to the same plan, by inspecting DICOM modality tags and cross-reference (Referenced UID) tags. The folder (and its sub-folders) may contain extra or unrelated files of any of these types; this function searches through all of them for the one set that is actually linked together.

Return type:

dict

Returns:

dict with keys – “dose”, “let”, “rtstruct”, “rtplan” : Path or None “CT” : Path to the directory holding the matched CT

series, or None

”Patient_ID” : str or None “linked” : True if RTPLAN -> RTSTRUCT -> CT was fully

confirmed via DICOM reference tags; False if a fallback/best-guess match had to be used; None if there wasn’t enough data to even attempt the check (e.g. no RTSTRUCT and no CT found at all).

”link_warnings”: list[str], one entry per fallback

that was needed, explaining what couldn’t be confirmed and what was used instead.

Parameters:

• folder (Path)

• rad_type (str)

load_ct_series(ct_folder)

Load a multi-slice CT DICOM series from a folder.

Slices are sorted by ImagePositionPatient z-coordinate. Spacing is taken from PixelSpacing of the first slice and the z-step between consecutive slice positions.

Return type:

tuple

Returns:

• sitk_image (SimpleITK.Image (x, y, z ordering internally))

• ct_geometry (dict with keys:) – “origin” : [x0, y0, z0] mm (corner of first voxel) “spacing” : [dx, dy, dz] mm “shape” : (nz, ny, nx) — numpy (z,y,x) convention “z_positions” : np.ndarray of slice z-coordinates length nz

Parameters:

ct_folder (str | Path)

resample_dose_to_new_grid(dose_sitk, dose_ds, new_spacing, interpolator=SimpleITK.sitkLinear)

Resample RTDOSE onto a new isotropic/anisotropic grid :type dose_sitk: sitk.Image Original dose image :param dose_sitk: :type dose_sitk: sitk.Image Original dose image :type dose_ds: pydicom Dataset RTDOSE dataset :param dose_ds: :type dose_ds: pydicom Dataset RTDOSE dataset :type new_spacing: tuple/list (sx, sy, sz) in mm :param new_spacing: :type new_spacing: tuple/list (sx, sy, sz) in mm

Returns:

• resampled_dose (sitk.Image)

• updated_info (dict)

• dose_ds (updated dataset)

resample_dose_on_ct(sitk_dose, sitk_ct)

Resample a dose (or LET) SimpleITK image onto the CT grid.

The CT image defines the output origin, spacing, direction, and size. This ensures that the resampled dose array is perfectly aligned with the CT grid on which contours will be rasterised.

Parameters:

• sitk_dose (SimpleITK.Image (dose or LET, in dose-grid coordinates))

• sitk_ct (SimpleITK.Image (full CT series))

Return type:

Image

Returns:

SimpleITK.Image same grid as sitk_ct

load_dose_grid(path)

Load a DICOM RT Dose file (dose or LET stored as dose grid).

Return type:

tuple

Returns:

• array (np.ndarray shape (z, y, x))

• ds (pydicom Dataset)

Parameters:

path (str | Path)

get_grid_geometry(ds)

Extract (origin, spacing) from an RT Dose dataset.

Return type:

tuple

Returns:

• origin ([x0, y0, z0] mm)

• spacing ([dx, dy, dz] mm)

get_all_structure_names(rtstruct_ds)

Return list of all structure names in an RT Struct dataset.

Return type:

list

get_structure_mask_on_grid(struct_name, rtstruct_ds, origin, spacing, shape, z_positions)

Rasterise RT Struct contours as BINARY MASKS for struct_name onto an arbitrary grid.

This is the core function used for both dose-grid and CT-grid masking. Contour z-values are matched to the nearest z in z_positions.

Parameters:

• struct_name (str)

• rtstruct_ds (pydicom Dataset)

• origin ([x0, y0, z0] mm — physical coordinate of voxel (0,0,0) corner, then converted to voxel center)

• spacing ([dx, dy, dz] mm)

• shape ((nz, ny, nx) — numpy array shape)

• z_positions (1-D array of z-coordinates for each slice (length nz))

Return type:

ndarray

Returns:

mask (np.ndarray bool, shape (nz, ny, nx))

get_roi_center_of_mass(struct_name, rtstruct_ds)

Compute the 3-D center of mass (mm, in the RT Struct’s patient coordinate system) of a structure directly from its contour polygons — no dose/CT grid required.

Each contour (one per slice) contributes its own 2-D polygon centroid, weighted by that contour’s area, so slices with more cross-sectional area count more towards the overall COM volume-weighted approximation without needing a full 3-D mask.

Parameters:

• struct_name (str) – Structure name as it appears in the RT Struct.

• rtstruct_ds (pydicom Dataset)

Return type:

ndarray

Returns:

np.ndarray, shape (3,), or None – [x, y, z] center of mass in mm. Returns None if no usable contour data is available for this structure.

contour_area_signed(xy)

Signed polygon area via the shoelace formula.

Parameters:

xy (np.ndarray, shape (N, 2)) – Polygon vertices [x, y] in mm.

Returns:

float – Signed area; positive for counter-clockwise vertex order, negative for clockwise. Used to detect contour holes.

smooth_contour(poly_xy, n_pts=300)

Resample a contour polygon onto n_pts evenly-spaced points using a periodic B-spline fit. Not used by default in the fractional-mask pipeline (see inline note in get_fractional_mask_on_grid); kept available for callers who want smoothed contours.

Parameters:

• poly_xy (np.ndarray, shape (N, 2)) – Polygon vertices [x, y] in mm.

• n_pts (int, default 300) – Number of points in the resampled output.

Returns:

np.ndarray, shape (n_pts, 2) – Smoothed polygon vertices. Falls back to the original polygon, unchanged, if the spline fit fails.

rasterize_supersampled(smooth_xy, x0, y0, dx, dy, ny, nx, N)

Rasterise one contour polygon onto a grid at N times finer resolution, then average back down — giving each output voxel a fractional [0, 1] membership value instead of a hard binary in/out.

Parameters:

• smooth_xy (np.ndarray, shape (N, 2)) – Polygon vertices [x, y] in mm.

• x0 (float) – Grid origin (mm) in x and y.

• y0 (float) – Grid origin (mm) in x and y.

• dx (float) – Grid voxel spacing (mm) in x and y.

• dy (float) – Grid voxel spacing (mm) in x and y.

• ny (int) – Output grid shape.

• nx (int) – Output grid shape.

• N (int) – Supersampling factor per side (N² sub-samples per voxel).

Returns:

np.ndarray, shape (ny, nx), float32 – Fractional coverage of each voxel by the polygon, in [0, 1].

compute_roi_volume_comparison(frac_mask, dx, dy, dz)

Print and return the total volume of a fractional mask, plus a per-slice breakdown — useful for sanity-checking get_fractional_mask_on_grid output against a TPS-reported volume.

Parameters:

• frac_mask (np.ndarray, shape (nz, ny, nx)) – Fractional voxel membership mask, values in [0, 1].

• dx (float) – Voxel spacing (mm) in x, y, z.

• dy (float) – Voxel spacing (mm) in x, y, z.

• dz (float) – Voxel spacing (mm) in x, y, z.

Returns:

float – Total volume in cc (sum of fractional weights x voxel volume).

prismatoid_volume(frac_mask, dx, dy, dz)

Volume of a fractional mask via the prismatoid (Simpson’s-rule-like) formula between consecutive slices, instead of a flat sum-of-slices approximation. Slightly more accurate for structures with rapidly changing cross-sectional area between slices.

Parameters:

• frac_mask (np.ndarray, shape (nz, ny, nx)) – Fractional voxel membership mask, values in [0, 1].

• dx (float) – Voxel spacing (mm) in x, y, z.

• dy (float) – Voxel spacing (mm) in x, y, z.

• dz (float) – Voxel spacing (mm) in x, y, z.

Returns:

float – Total volume in cc.

compute_roi_volume_from_contours(struct_name, rtstruct_ds)

Compute the ROI’s total volume (cc) directly from its contour polygon geometry, at the NATIVE contour z-spacing — independent of any dose/LET/analysis grid resolution.

TECHNICAL JUSTIFICATION: get_fractional_mask_on_grid (and hence compute_dvh_metrics/compute_let_metrics in metrics.py, which normalise Vx%/Dx%/Lx% by the sum of that mask’s fractional weights) reports volume on whatever analysis grid the mask was rasterised on (e.g. a resampled dose/LET grid, New_grid — commonly coarser than the CT the structure was actually contoured on). For structures that only span a handful of voxels across that grid (small serial OARs: optic nerves, chiasm, lenses, lacrimal glands), that grid-resolution dependence is the dominant source of volume error, which then propagates into every weighted DVH/LVH percentile metric for that structure (they all share the same total-weight normalisation).

This function instead reconstructs volume the way most TPS systems report “ROI volume” — directly from the contour polygons themselves: net in-plane area per true contour slice (_slice_net_area_mm2, correctly excluding holes), trapezoidal integration between consecutive TRUE contour z-planes at their native (possibly irregular) spacing, plus a half-native-spacing end-cap extension at each pole — the same physical convention get_fractional_mask_on_grid uses for its own z-extent padding, so the two stay consistent with each other even though this number is otherwise grid-independent.

Use this as the authoritative reported volume_cc for a structure. Note that it is NOT automatically a safe drop-in replacement for the weighted-percentile normalisation used inside compute_dvh_metrics/ compute_let_metrics (Vx%, Dx%, Lx%, …): those numerators are still computed from the analysis-grid mask, so swapping only the denominator to a different, higher-resolution source is only self-consistent once the grid-captured voxel population is itself a low-bias subsample of the true structure (see the exact z-integration in get_fractional_mask_on_grid, which is what makes that assumption reasonable in the first place).

Parameters:

• struct_name (str)

• rtstruct_ds (pydicom Dataset)

Return type:

float

Returns:

float — ROI volume in cc.

rasterize_slice_coverage(slice_polys, x0, y0, dx, dy, ny, nx, N)

Rasterise ALL polygon loops belonging to ONE true contour z-position (i.e. one DICOM ContourData slice) into a single in-plane fractional coverage map, islands added and holes subtracted per _classify_holes_by_nesting.

Factored out of get_fractional_mask_on_grid so that in-plane coverage can be computed once per ORIGINAL contour z (see that function’s docstring for why this must be decoupled from the output grid’s z-slices).

Parameters:

• slice_polys (list of np.ndarray, shape (M, 2)) – Polygon loops (islands + holes) for one contour slice, in mm.

• x0 (see rasterize_supersampled.)

• y0 (see rasterize_supersampled.)

• dx (see rasterize_supersampled.)

• dy (see rasterize_supersampled.)

• ny (see rasterize_supersampled.)

• nx (see rasterize_supersampled.)

• N (see rasterize_supersampled.)

Returns:

np.ndarray, shape (ny, nx), float32 – Fractional in-plane coverage, clipped to [0, 1].

get_fractional_mask_on_grid(struct_name, rtstruct_ds, origin, spacing, shape, z_positions, supersample=4, max_supersample=12)

Compute a fractional voxel membership mask on an arbitrary grid.

Each voxel receives a value in [0,1] — the fraction of its physical VOLUME (in-plane AND through-slice) that lies inside the RT Struct contour. In-plane coverage is estimated by supersampling (adaptive N² sub-points per voxel, see _adaptive_supersample). Through-slice (z) coverage is estimated by treating the true contour z-positions as control planes and interpolating the STRUCTURE’S BOUNDARY GEOMETRY — not its raw in-plane coverage/opacity — between the bracketing contour slices across each output voxel’s z-extent. This is done via shape-based interpolation: each contour slice’s fractional coverage raster is first converted to a 2-D signed distance field (SDF, see _coverage_to_sdf); the per-voxel z-slab is then integrated EXACTLY in SDF space (_slab_avg_sdf_exact — see that function’s docstring, and NOTE below), not by sub-sampling; and the resulting averaged SDF is converted back to a fractional coverage raster only once, at the end (_sdf_to_coverage).

NOTE on the z-integration being exact rather than sub-sampled: the SDF field is piecewise LINEAR in z by construction (linear between each pair of adjacent true contour planes, constant beyond the first/last). The integral average of a piecewise-linear function over any interval is available in closed form from its endpoint values at each linear segment (trapezoid rule, exact for linear functions) — so no Nz sub-sampling parameter is needed, and there is zero z-quadrature error regardless of how coarse the output grid’s dz is relative to the native contour spacing, or how many true contour planes fall inside one output voxel’s z-slab. (A generic fixed-Nz sub-sampling grid — the previous approach — is only exact when a slab happens to fall entirely within one linear segment; it carries residual error whenever a slab straddles an interior contour plane, which is exactly the case where the output grid is coarser than the contour spacing.) The exact integration also correctly handles voxels whose z-extent only PARTIALLY overlaps the structure’s true (padded) extent — e.g. the ROI’s first/last output slice — by scaling the covered region’s coverage down by the actual overlap fraction, rather than reporting the covered portion’s coverage as if it applied to the voxel’s full z-extent.

TECHNICAL JUSTIFICATION for the z-interpolation (this is the main change from the previous version): contours were previously snapped to the single nearest output z-index, so an output voxel’s z-extent was always either “fully in” or “fully out” of the structure — no different, in the z direction, from a hard binary mask. This reproduces the same staircase/quantization error along the cranial-caudal axis that fractional in-plane weighting was written to avoid in-plane, and it gets worse whenever the analysis grid’s z spacing does not match the original contour spacing (e.g. resampling dose/LET/mask onto a coarser custom grid).

A direct linear blend of the raw coverage rasters (alpha cross-dissolving) was tried first, but that interpolates OPACITY, not GEOMETRY: wherever the cross-section shifts, rotates, or changes size between two contour planes (tapering structures, branching anatomy, slightly mis-registered slices), an opacity blend produces a faint, eroded double-exposure of both shapes rather than a smoothly moving boundary, understating volume in exactly the transition regions where accuracy matters most. Interpolating signed distance fields instead — the classical “shape-based interpolation” approach (Raya & Udupa, 1990) — averages DISTANCE TO THE BOUNDARY rather than opacity, so the reconstructed boundary at any intermediate depth is the correct (piecewise-linear-in-distance) locus of points equidistant between the two true contours; this is the discretised equivalent of lofting a triangulated surface between them, without an opacity blend’s tendency to thin the structure out.

Outside the structure’s own z-extent, the boundary geometry (i.e. the SDF, hence coverage) is held constant at the nearest edge slice’s value (rather than tapering to zero), because RTSTRUCT contours conventionally represent the slice-thickness worth of structure centred on each contour, so the ROI is assumed to occupy a half-slice-thickness beyond the first/last contour plane, not to end exactly on it.

Parameters:

• struct_name (str)

• rtstruct_ds (pydicom Dataset)

• origin ([x0, y0, z0] mm)

• spacing ([dx, dy, dz] mm)

• shape ((nz, ny, nx))

• z_positions (1-D array length nz)

• supersample (baseline in-plane N (default 4). Scaled up) – per-structure by _adaptive_supersample.

• max_supersample (hard cap for the adaptive in-plane N.)

Return type:

ndarray

Returns:

frac_mask (np.ndarray float32, shape (nz, ny, nx), values in [0, 1])

save_nifti(arr, ref_img, path)