From 916a05e175e12f491879bc30ad677a6083452dd1 Mon Sep 17 00:00:00 2001 From: Denis Samatov <98321287+denis-samatov@users.noreply.github.com> Date: Tue, 13 Oct 2026 14:03:05 +0700 Subject: [PATCH 1/5] Add validated CT/SEG tensor and patient-split tutorial --- README.md | 1 + .../advanced_topics/ct_seg_to_tensors.ipynb | 880 ++++++++++++++++++ test/src/test_ct_seg_to_tensors.py | 184 ++++ 3 files changed, 1065 insertions(+) create mode 100644 notebooks/advanced_topics/ct_seg_to_tensors.ipynb create mode 100644 test/src/test_ct_seg_to_tensors.py diff --git a/README.md b/README.md index b9b6fed..fa88d90 100644 --- a/README.md +++ b/README.md @@ -50,6 +50,7 @@ Zenodo permanent citeable record: https://zenodo.org/records/22691908. Notebooks in this folder focus on topics that will require understanding of the basics, and aim to address more narrow use cases of IDC usage. * [Searching DICOM private tags](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/dicom_private_tags_intro.ipynb): all of DICOM attributes for the imaging data in IDC are searchable using BigQuery. DICOM private tags often contain critical information, such as diffusion b-values, but are a bit more tricky to access from BigQuery. In this tutorial you will learn how to accomplish this. +* [Validated CT/DICOM SEG tensors and patient-disjoint splits](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/ct_seg_to_tensors.ipynb): select a small C4KC-KiTS cohort from the pinned IDC catalog, validate source references and spatial correspondence, export HU images and separate mask channels for PyTorch, and create a reproducible split by patient. * [Working with IDC clinical data](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/clinical_data_intro.ipynb): this tutorial will help you become familiar with how clinical and other tabular data accompanying images in IDC is organized, and how you can use it to build cohorts. * [Direct access to tags/frames from GCS/AWS buckets](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/gcs_aws_direct_access.ipynb): learn how to access individual frames or tags of large DICOM files from the bucket without having to download the entire file (this notebook accompanies documentation article here: https://learn.canceridc.dev/data/downloading-data/direct-loading) * [Using DICOMweb to access IDC data](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/idc_dicomweb_access.ipynb): both metadata and pixel data can be access using DICOMweb, which is particularly important while working with digital pathology, as it enables granular access to the individual pyramid tiles (frames) diff --git a/notebooks/advanced_topics/ct_seg_to_tensors.ipynb b/notebooks/advanced_topics/ct_seg_to_tensors.ipynb new file mode 100644 index 0000000..7e347c2 --- /dev/null +++ b/notebooks/advanced_topics/ct_seg_to_tensors.ipynb @@ -0,0 +1,880 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "87781114", + "metadata": {}, + "source": [ + "\"Open\n", + "\n", + "# Validated CT/DICOM SEG tensors and patient-disjoint splits\n", + "\n", + "Updated: October 2026\n", + "\n", + "## Goal\n", + "\n", + "Turn a small public **C4KC-KiTS** cohort into traceable PyTorch inputs: an HU image\n", + "`[1, D, H, W]`, separate binary mask channels `[K, D, H, W]`, an affine in DICOM\n", + "patient coordinates, and a patient-disjoint train/test manifest.\n", + "\n", + "The [IDC Hugging Face dataset](https://huggingface.co/datasets/ImagingDataCommons/idc-index-data)\n", + "is a **metadata catalog**, not pixel data. Its `train` split is a Hub convention,\n", + "not a medical ML training partition. This notebook queries a pinned catalog and\n", + "uses `idc-index` to download only selected DICOM series. No cloud account or GPU\n", + "is required. Start with [Getting Started Part 2](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/getting_started/part2_searching_basics.ipynb)\n", + "if you are new to IDC.\n", + "\n", + "### Scope and assumptions\n", + "\n", + "- Four patients, selected deterministically by the smallest combined CT/SEG download.\n", + " This deliberately biased convenience sample is for learning, not model evaluation.\n", + "- Legacy single-frame CT on a regular grid, linear rescaling to HU, and one BINARY\n", + " SEG instance per series. Enhanced CT, non-regular spacing, gantry tilt, and\n", + " resampled masks are rejected. Nothing is silently resampled.\n", + "- CT is selected by the SEG's **referenced series**, then verified against source\n", + " SOP Instance UIDs and per-frame geometry. A patient/study match is insufficient.\n", + "- Separate channels retain overlaps. Coded category/type, modifiers and anatomic\n", + " region define a stable channel order across cases; free-text labels and local\n", + " segment numbers alone do not define classes.\n", + "- Missing CT files are errors. Omitted mask frames are also errors unless the\n", + " caller explicitly enables an empty-frame policy for a verified dataset.\n", + "- The instructional split is within the public KiTS19 training collection; it is\n", + " not the official challenge test set. Cross-collection identity reconciliation\n", + " and statistically representative sampling are outside this example." + ] + }, + { + "cell_type": "markdown", + "id": "af3feb50", + "metadata": {}, + "source": [ + "## Setup\n", + "\n", + "Install the versions used for this example in a fresh Python 3.11/3.12 environment.\n", + "The PyTorch wheel is CPU-only. Local Jupyter users also need `ipykernel` and\n", + "`nbclient` to execute the notebook programmatically." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cc43facd", + "metadata": {}, + "outputs": [], + "source": [ + "%pip install -q \"idc-index==0.13.0\" \"pydicom==3.0.2\" \"highdicom==0.28.2\" \"numpy==2.3.5\" \"pandas==2.2.3\" \"matplotlib==3.10.8\"\n", + "%pip install -q \"torch==2.14.1\" --index-url https://download.pytorch.org/whl/cpu" + ] + }, + { + "cell_type": "markdown", + "id": "b66ef835", + "metadata": {}, + "source": [ + "### Parameters\n", + "\n", + "The default uses Hugging Face revision `93983d810bea0454ae1a400ab62c48dd01f8d0ca`\n", + "of the IDC v25 catalog. Set `CATALOG_SOURCE = \"idc-index\"` to use the corresponding\n", + "package indices instead. Both routes are checked for matching release metadata.\n", + "\n", + "Metadata downloads are separate from DICOM download size estimates. All paths\n", + "below are relative to the working directory. Rerunning reuses existing DICOM\n", + "files and replaces generated outputs for the same sample." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "af68526c", + "metadata": {}, + "outputs": [], + "source": [ + "from importlib.metadata import version\n", + "from pathlib import Path\n", + "from time import perf_counter\n", + "import json\n", + "import hashlib\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", + "import pyarrow.parquet as pq\n", + "import requests\n", + "import torch\n", + "from idc_index import IDCClient\n", + "from IPython.display import display\n", + "\n", + "CATALOG_SOURCE = \"huggingface\" # or \"idc-index\"\n", + "HF_REVISION = \"93983d810bea0454ae1a400ab62c48dd01f8d0ca\"\n", + "COLLECTION = \"c4kc_kits\"\n", + "N_PATIENTS = 4\n", + "SPLIT_SEED = 2026\n", + "TEST_FRACTION = 0.25\n", + "MAX_DICOM_MB = 200.0\n", + "ALLOW_OMITTED_EMPTY_FRAMES = False\n", + "WORK_DIR = Path(\"idc_ct_seg_data\")\n", + "DICOM_DIR = WORK_DIR / \"dicom\"\n", + "OUTPUT_DIR = WORK_DIR / \"outputs\"\n", + "OUTPUT_DIR.mkdir(parents=True, exist_ok=True)\n", + "run_started = perf_counter()" + ] + }, + { + "cell_type": "markdown", + "id": "ff60e33c", + "metadata": {}, + "source": [ + "## Steps\n", + "\n", + "### 1. Define the validation boundary\n", + "\n", + "The following functions validate DICOM objects before turning them into arrays.\n", + "Geometric tolerances are 0.01 mm for positions/slice spacing and `1e-4` for\n", + "direction cosines/in-plane spacing. The affine maps array indices `[D, H, W]`\n", + "to **LPS millimeters**, not NIfTI RAS coordinates.\n", + "\n", + "C4KC-KiTS SEG objects do not set `SpatialLocationsPreserved`. The code tolerates\n", + "that absent flag only after independently checking every frame's geometry.\n", + "An explicit `NO` is always rejected. `SliceThickness` is not used as the slice\n", + "step: positions determine it." + ] + }, + { + "cell_type": "markdown", + "id": "5c7be955", + "metadata": {}, + "source": [ + "### Validation helpers" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "da62205f", + "metadata": { + "tags": [ + "definitions" + ] + }, + "outputs": [], + "source": [ + "from __future__ import annotations\n", + "\n", + "import hashlib\n", + "import json\n", + "import math\n", + "import random\n", + "from pathlib import Path\n", + "\n", + "import highdicom as hd\n", + "import numpy as np\n", + "import pydicom\n", + "from pydicom.uid import CTImageStorage\n", + "\n", + "\n", + "class ValidationError(ValueError):\n", + " \"\"\"A CT/SEG pair cannot be interpreted safely within this tutorial's scope.\"\"\"\n", + "\n", + "\n", + "def require(condition: bool, message: str) -> None:\n", + " if not condition:\n", + " raise ValidationError(message)\n", + "\n", + "\n", + "def numeric_vector(value, length: int, name: str) -> np.ndarray:\n", + " result = np.asarray(value, dtype=np.float64)\n", + " require(result.shape == (length,) and np.isfinite(result).all(), f\"Invalid {name}\")\n", + " return result" + ] + }, + { + "cell_type": "markdown", + "id": "b6631e97", + "metadata": {}, + "source": [ + "### Validate and order the CT volume" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c9dd0578", + "metadata": { + "tags": [ + "definitions" + ] + }, + "outputs": [], + "source": [ + "def sort_and_validate_ct(datasets, expected_count: int):\n", + " \"\"\"Return geometrically sorted single-frame CT instances and an LPS affine.\"\"\"\n", + " require(len(datasets) == expected_count, \"CT instance count differs from the catalog\")\n", + " require(len(datasets) >= 2, \"A volume needs at least two CT slices\")\n", + " first = datasets[0]\n", + " uids = [str(d.SOPInstanceUID) for d in datasets]\n", + " require(len(set(uids)) == len(uids), \"Duplicate CT SOP Instance UID\")\n", + " orientation = numeric_vector(first.ImageOrientationPatient, 6, \"CT orientation\")\n", + " row_direction, column_direction = orientation[:3], orientation[3:]\n", + " require(np.isclose(np.linalg.norm(row_direction), 1, atol=1e-4), \"Non-unit CT row direction\")\n", + " require(np.isclose(np.linalg.norm(column_direction), 1, atol=1e-4), \"Non-unit CT column direction\")\n", + " require(abs(np.dot(row_direction, column_direction)) < 1e-4, \"Non-orthogonal CT directions\")\n", + " normal = np.cross(row_direction, column_direction)\n", + " normal /= np.linalg.norm(normal)\n", + " spacing = numeric_vector(first.PixelSpacing, 2, \"CT pixel spacing\")\n", + " require((spacing > 0).all(), \"CT pixel spacing must be positive\")\n", + " positions = []\n", + " for ds in datasets:\n", + " require(str(ds.SOPClassUID) == str(CTImageStorage), \"Only legacy single-frame CT is supported\")\n", + " require(int(ds.get(\"NumberOfFrames\", 1)) == 1, \"Enhanced/multiframe CT is unsupported\")\n", + " require(ds.Modality == \"CT\", \"Non-CT instance in the CT series\")\n", + " for key in (\"PatientID\", \"StudyInstanceUID\", \"SeriesInstanceUID\", \"FrameOfReferenceUID\"):\n", + " require(str(ds.get(key, \"\")) == str(first.get(key, \"\")) != \"\", f\"Inconsistent CT {key}\")\n", + " require((ds.Rows, ds.Columns) == (first.Rows, first.Columns), \"Inconsistent CT dimensions\")\n", + " require(np.allclose(numeric_vector(ds.ImageOrientationPatient, 6, \"CT orientation\"), orientation, rtol=0, atol=1e-4), \"CT orientation varies\")\n", + " require(np.allclose(numeric_vector(ds.PixelSpacing, 2, \"CT pixel spacing\"), spacing, rtol=0, atol=1e-4), \"CT pixel spacing varies\")\n", + " positions.append(numeric_vector(ds.ImagePositionPatient, 3, \"CT position\"))\n", + " positions = np.stack(positions)\n", + " order = np.argsort(positions @ normal, kind=\"stable\")\n", + " positions = positions[order]\n", + " steps = np.diff(positions @ normal)\n", + " require((steps > 0.01).all(), \"Duplicate or nearly duplicate CT slice positions\")\n", + " slice_spacing = float(np.median(steps))\n", + " require(np.allclose(steps, slice_spacing, rtol=0, atol=0.01), \"Irregular CT slice spacing or missing slices\")\n", + " expected_positions = positions[0] + np.arange(len(datasets))[:, None] * slice_spacing * normal\n", + " require(np.allclose(positions, expected_positions, rtol=0, atol=0.01), \"Gantry tilt or non-regular CT grid is unsupported\")\n", + " affine = np.eye(4)\n", + " # Array indices [slice, row, column] -> DICOM patient coordinates [L, P, S], mm.\n", + " affine[:3, 0] = normal * slice_spacing\n", + " affine[:3, 1] = column_direction * spacing[0]\n", + " affine[:3, 2] = row_direction * spacing[1]\n", + " affine[:3, 3] = positions[0]\n", + " return [datasets[int(i)] for i in order], affine" + ] + }, + { + "cell_type": "markdown", + "id": "715888b1", + "metadata": {}, + "source": [ + "### Check every SEG frame against its referenced CT" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "01ae8c29", + "metadata": { + "tags": [ + "definitions" + ] + }, + "outputs": [], + "source": [ + "def functional_group(frame, shared, sequence_name):\n", + " sequence = frame.get(sequence_name, shared.get(sequence_name, None))\n", + " require(sequence is not None and len(sequence) == 1, f\"Missing/ambiguous SEG {sequence_name}\")\n", + " return sequence[0]\n", + "\n", + "\n", + "def validate_seg_geometry(seg, ct_datasets, allow_omitted_empty_frames: bool = False):\n", + " \"\"\"Verify source references and every SEG frame on the CT grid before decoding.\"\"\"\n", + " require(seg.Modality == \"SEG\" and str(seg.SegmentationType) == \"BINARY\", \"Only BINARY DICOM SEG is supported\")\n", + " ct = ct_datasets[0]\n", + " for key in (\"PatientID\", \"StudyInstanceUID\", \"FrameOfReferenceUID\"):\n", + " require(str(seg.get(key, \"\")) == str(ct.get(key, \"\")) != \"\", f\"CT/SEG {key} mismatch\")\n", + " require((seg.Rows, seg.Columns) == (ct.Rows, ct.Columns), \"CT/SEG dimensions differ\")\n", + " series_refs = seg.get(\"ReferencedSeriesSequence\", [])\n", + " require(len(series_refs) == 1, \"SEG must reference exactly one source series\")\n", + " require(str(series_refs[0].SeriesInstanceUID) == str(ct.SeriesInstanceUID), \"SEG references a different CT series\")\n", + " referenced_instances = series_refs[0].get(\"ReferencedInstanceSequence\", [])\n", + " referenced_uids = [str(item.ReferencedSOPInstanceUID) for item in referenced_instances]\n", + " require(len(referenced_uids) == len(set(referenced_uids)), \"Duplicate SEG source reference\")\n", + " ct_by_uid = {str(ds.SOPInstanceUID): ds for ds in ct_datasets}\n", + " require(set(referenced_uids) == set(ct_by_uid), \"SEG source references and complete CT series do not match\")\n", + " require(all(str(item.ReferencedSOPClassUID) == str(CTImageStorage) for item in referenced_instances), \"SEG references unsupported source classes\")\n", + " shared_sequence = seg.get(\"SharedFunctionalGroupsSequence\", [])\n", + " require(len(shared_sequence) == 1, \"SEG needs one Shared Functional Groups item\")\n", + " shared = shared_sequence[0]\n", + " frames = seg.get(\"PerFrameFunctionalGroupsSequence\", [])\n", + " require(len(frames) == int(seg.NumberOfFrames), \"SEG functional-group frame count mismatch\")\n", + " segment_numbers = [int(item.SegmentNumber) for item in seg.SegmentSequence]\n", + " require(segment_numbers and len(set(segment_numbers)) == len(segment_numbers), \"Invalid SEG segment numbering\")\n", + " seen = set()\n", + " preserved_flags = []\n", + " for frame in frames:\n", + " derivations = frame.get(\"DerivationImageSequence\", shared.get(\"DerivationImageSequence\", []))\n", + " sources = [source for derivation in derivations for source in derivation.get(\"SourceImageSequence\", [])]\n", + " require(len(sources) == 1, \"Each SEG frame must reference one single-frame CT instance\")\n", + " source = sources[0]\n", + " uid = str(source.ReferencedSOPInstanceUID)\n", + " require(uid in ct_by_uid, \"SEG frame references a CT instance that is missing\")\n", + " require(str(source.ReferencedSOPClassUID) == str(CTImageStorage), \"SEG frame source is not single-frame CT\")\n", + " require(\"ReferencedFrameNumber\" not in source, \"Multiframe source references are unsupported\")\n", + " preserved = str(source.get(\"SpatialLocationsPreserved\", \"\"))\n", + " require(preserved in (\"\", \"YES\"), \"SEG explicitly does not preserve spatial locations\")\n", + " preserved_flags.append(preserved)\n", + " source_ct = ct_by_uid[uid]\n", + " orientation = functional_group(frame, shared, \"PlaneOrientationSequence\").ImageOrientationPatient\n", + " position = functional_group(frame, shared, \"PlanePositionSequence\").ImagePositionPatient\n", + " pixel_spacing = functional_group(frame, shared, \"PixelMeasuresSequence\").PixelSpacing\n", + " require(np.allclose(numeric_vector(orientation, 6, \"SEG orientation\"), source_ct.ImageOrientationPatient, rtol=0, atol=1e-4), \"SEG orientation differs from its source CT\")\n", + " require(np.allclose(numeric_vector(position, 3, \"SEG position\"), source_ct.ImagePositionPatient, rtol=0, atol=0.01), \"SEG position differs from its source CT\")\n", + " require(np.allclose(numeric_vector(pixel_spacing, 2, \"SEG pixel spacing\"), source_ct.PixelSpacing, rtol=0, atol=1e-4), \"SEG pixel spacing differs from its source CT\")\n", + " number = int(functional_group(frame, shared, \"SegmentIdentificationSequence\").ReferencedSegmentNumber)\n", + " require(number in segment_numbers, \"SEG frame names an undefined segment\")\n", + " require((uid, number) not in seen, \"Duplicate SEG frame for one source instance and segment\")\n", + " seen.add((uid, number))\n", + " missing = {(uid, number) for uid in ct_by_uid for number in segment_numbers} - seen\n", + " require(not missing or allow_omitted_empty_frames, \"SEG omits source/segment frames; explicit empty-frame policy is required\")\n", + " return {\n", + " \"segment_numbers\": segment_numbers,\n", + " \"omitted_source_segment_frames\": len(missing),\n", + " \"spatial_locations_flag_absent\": any(flag == \"\" for flag in preserved_flags),\n", + " }" + ] + }, + { + "cell_type": "markdown", + "id": "615e8f2a", + "metadata": {}, + "source": [ + "### Preserve coded segment semantics" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2751c843", + "metadata": { + "tags": [ + "definitions" + ] + }, + "outputs": [], + "source": [ + "def code_metadata(code):\n", + " return {\n", + " \"value\": str(code.CodeValue),\n", + " \"scheme\": str(code.CodingSchemeDesignator),\n", + " \"meaning\": str(code.CodeMeaning),\n", + " }\n", + "\n", + "\n", + "def segment_metadata(seg):\n", + " result = []\n", + " for item in seg.SegmentSequence:\n", + " category = code_metadata(item.SegmentedPropertyCategoryCodeSequence[0])\n", + " property_type = code_metadata(item.SegmentedPropertyTypeCodeSequence[0])\n", + " modifiers = [code_metadata(c) for c in item.SegmentedPropertyTypeCodeSequence[0].get(\"SegmentedPropertyTypeModifierCodeSequence\", [])]\n", + " regions = []\n", + " for region in item.get(\"AnatomicRegionSequence\", []):\n", + " regions.append({\"code\": code_metadata(region), \"modifiers\": [code_metadata(c) for c in region.get(\"AnatomicRegionModifierSequence\", [])]})\n", + " # Include modifiers/anatomic region in the semantic key, not just the free-text label.\n", + " def identity(code):\n", + " return (code[\"scheme\"], code[\"value\"])\n", + "\n", + " semantic_key = json.dumps({\n", + " \"category\": identity(category), \"type\": identity(property_type),\n", + " \"type_modifiers\": sorted(identity(code) for code in modifiers),\n", + " \"regions\": sorted((identity(region[\"code\"]), sorted(identity(code) for code in region[\"modifiers\"])) for region in regions),\n", + " }, sort_keys=True)\n", + " result.append({\n", + " \"segment_number\": int(item.SegmentNumber),\n", + " \"label\": str(item.SegmentLabel),\n", + " \"description\": str(item.get(\"SegmentDescription\", \"\")),\n", + " \"algorithm_type\": str(item.SegmentAlgorithmType),\n", + " \"algorithm_name\": str(item.get(\"SegmentAlgorithmName\", \"\")),\n", + " \"category\": category, \"property_type\": property_type,\n", + " \"type_modifiers\": modifiers, \"anatomic_regions\": regions,\n", + " \"semantic_key\": semantic_key,\n", + " })\n", + " require(len({item[\"semantic_key\"] for item in result}) == len(result), \"Ambiguous duplicate semantic segment definitions\")\n", + " return sorted(result, key=lambda item: item[\"semantic_key\"])" + ] + }, + { + "cell_type": "markdown", + "id": "10236b4e", + "metadata": {}, + "source": [ + "### Convert a validated pair" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cf41526b", + "metadata": { + "tags": [ + "definitions" + ] + }, + "outputs": [], + "source": [ + "def validate_and_convert(ct_datasets, seg, expected_count: int, allow_omitted_empty_frames: bool = False):\n", + " \"\"\"Validate one real DICOM pair; return HU image, separate masks, and provenance.\"\"\"\n", + " ct_datasets, affine = sort_and_validate_ct(ct_datasets, expected_count)\n", + " checks = validate_seg_geometry(seg, ct_datasets, allow_omitted_empty_frames)\n", + " segments = segment_metadata(seg)\n", + " source_uids = [str(ds.SOPInstanceUID) for ds in ct_datasets]\n", + " # C4KC-KiTS does not set SpatialLocationsPreserved. Geometry was checked independently\n", + " # above; only that missing flag is ignored. Explicit NO has already been rejected.\n", + " masks_dhwk = seg.get_pixels_by_source_instance(\n", + " source_sop_instance_uids=source_uids,\n", + " segment_numbers=[item[\"segment_number\"] for item in segments],\n", + " combine_segments=False,\n", + " ignore_spatial_locations=checks[\"spatial_locations_flag_absent\"],\n", + " assert_missing_frames_are_empty=allow_omitted_empty_frames,\n", + " dtype=np.uint8,\n", + " )\n", + " require(masks_dhwk.shape == (len(ct_datasets), ct_datasets[0].Rows, ct_datasets[0].Columns, len(segments)), \"Decoded SEG shape differs from CT\")\n", + " require(np.isin(masks_dhwk, (0, 1)).all(), \"Non-binary decoded mask values\")\n", + " image_hu = np.empty((1, len(ct_datasets), ct_datasets[0].Rows, ct_datasets[0].Columns), dtype=np.float32)\n", + " for i, ds in enumerate(ct_datasets):\n", + " require(\"ModalityLUTSequence\" not in ds, \"This example supports linear CT rescaling only\")\n", + " require(str(ds.get(\"RescaleType\", \"HU\")).upper() == \"HU\", \"CT rescale units are not HU\")\n", + " require(\"RescaleSlope\" in ds and \"RescaleIntercept\" in ds, \"CT rescale parameters are missing\")\n", + " slope, intercept = float(ds.RescaleSlope), float(ds.RescaleIntercept)\n", + " require(math.isfinite(slope) and math.isfinite(intercept) and slope != 0, \"Invalid CT rescale parameters\")\n", + " image_hu[0, i] = ds.pixel_array.astype(np.float32) * slope + intercept\n", + " require(np.isfinite(image_hu).all(), \"Non-finite CT intensities after rescaling\")\n", + " masks = np.ascontiguousarray(np.moveaxis(masks_dhwk, -1, 0))\n", + " first = ct_datasets[0]\n", + " metadata = {\n", + " \"patient_id\": str(first.PatientID),\n", + " \"study_uid\": str(first.StudyInstanceUID),\n", + " \"ct_series_uid\": str(first.SeriesInstanceUID),\n", + " \"seg_series_uid\": str(seg.SeriesInstanceUID),\n", + " \"seg_sop_instance_uid\": str(seg.SOPInstanceUID),\n", + " \"frame_of_reference_uid\": str(first.FrameOfReferenceUID),\n", + " \"source_sop_instance_uids_in_slice_order\": source_uids,\n", + " \"array_axes\": [\"slice\", \"row\", \"column\"],\n", + " \"world_coordinate_system\": \"DICOM LPS, millimeters\",\n", + " \"affine_dhw_to_lps_mm\": affine.tolist(),\n", + " \"image_units\": \"HU\", \"segments\": segments, \"checks\": checks,\n", + " \"image_shape\": list(image_hu.shape), \"mask_shape\": list(masks.shape),\n", + " \"nonzero_voxels_by_channel\": [int(np.count_nonzero(mask)) for mask in masks],\n", + " }\n", + " return image_hu, masks, metadata" + ] + }, + { + "cell_type": "markdown", + "id": "1152e91f", + "metadata": {}, + "source": [ + "### Patient grouping and source checksums" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "286f15e7", + "metadata": { + "tags": [ + "definitions" + ] + }, + "outputs": [], + "source": [ + "def patient_disjoint_split(records, seed: int = 2026, test_fraction: float = 0.25):\n", + " \"\"\"Assign all cases for each collection/patient group to the same partition.\"\"\"\n", + " require(0 < test_fraction < 1, \"test_fraction must be between 0 and 1\")\n", + " for record in records:\n", + " require(bool(record.get(\"collection_id\")) and bool(record.get(\"patient_id\")), \"Missing patient group identity\")\n", + " groups = sorted({(record[\"collection_id\"], record[\"patient_id\"]) for record in records})\n", + " require(len(groups) >= 2, \"At least two patient groups are needed\")\n", + " random.Random(seed).shuffle(groups)\n", + " n_test = min(len(groups) - 1, max(1, math.ceil(len(groups) * test_fraction)))\n", + " test_groups = set(groups[:n_test])\n", + " result = [dict(record, split=\"test\" if (record[\"collection_id\"], record[\"patient_id\"]) in test_groups else \"train\") for record in records]\n", + " train = {(r[\"collection_id\"], r[\"patient_id\"]) for r in result if r[\"split\"] == \"train\"}\n", + " test = {(r[\"collection_id\"], r[\"patient_id\"]) for r in result if r[\"split\"] == \"test\"}\n", + " require(train and test and train.isdisjoint(test), \"Patient split validation failed\")\n", + " return result\n", + "\n", + "\n", + "def file_sha256(path: Path) -> str:\n", + " digest = hashlib.sha256()\n", + " with path.open(\"rb\") as stream:\n", + " for chunk in iter(lambda: stream.read(1024 * 1024), b\"\"):\n", + " digest.update(chunk)\n", + " return digest.hexdigest()" + ] + }, + { + "cell_type": "markdown", + "id": "1a6c69cc", + "metadata": {}, + "source": [ + "### 2. Load a pinned metadata catalog\n", + "\n", + "The main catalog has one row per series. `seg_index` supplies the referenced\n", + "CT series. We compare selected cloud object identifiers with the download\n", + "client's package index, so a catalog/client mismatch fails before downloading.\n", + "The Hugging Face option downloads two Parquet metadata files and records their\n", + "SHA256 checksums; this can be larger than the small DICOM example itself." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e7480611", + "metadata": {}, + "outputs": [], + "source": [ + "client = IDCClient()\n", + "require(client.get_idc_version() == \"v25\", \"This example is pinned to IDC v25\")\n", + "catalog_paths = {}\n", + "if CATALOG_SOURCE == \"huggingface\":\n", + " cache_dir = WORK_DIR / \"catalog\" / HF_REVISION\n", + " cache_dir.mkdir(parents=True, exist_ok=True)\n", + " for name in (\"idc_index\", \"seg_index\"):\n", + " path = cache_dir / f\"{name}.parquet\"\n", + " if not path.exists():\n", + " url = f\"https://huggingface.co/datasets/ImagingDataCommons/idc-index-data/resolve/{HF_REVISION}/{name}.parquet\"\n", + " temporary = path.with_suffix(\".partial\")\n", + " with requests.get(url, stream=True, timeout=(15, 90)) as response:\n", + " response.raise_for_status()\n", + " with temporary.open(\"wb\") as stream:\n", + " for chunk in response.iter_content(1024 * 1024):\n", + " stream.write(chunk)\n", + " temporary.replace(path)\n", + " catalog_paths[name] = path\n", + " index_table = pd.read_parquet(catalog_paths[\"idc_index\"])\n", + " seg_table = pd.read_parquet(catalog_paths[\"seg_index\"])\n", + "elif CATALOG_SOURCE == \"idc-index\":\n", + " client.fetch_index(\"seg_index\")\n", + " require(getattr(client, \"seg_index\", None) is not None, \"SEG index download failed\")\n", + " index_table, seg_table = client.index, client.seg_index\n", + " catalog_paths = {\"idc_index\": Path(client.indices_overview[\"index\"][\"file_path\"]),\n", + " \"seg_index\": Path(client.indices_overview[\"seg_index\"][\"file_path\"])}\n", + "else:\n", + " raise ValidationError(\"Unknown CATALOG_SOURCE\")\n", + "\n", + "catalog_provenance = {}\n", + "for name, path in catalog_paths.items():\n", + " generated_version = pq.read_metadata(path).metadata[b\"idc_index_data_version\"].decode()\n", + " require(generated_version.split(\"-\")[0] == version(\"idc-index-data\"), \"Catalog release differs from the download client\")\n", + " catalog_provenance[name] = {\"sha256\": file_sha256(path), \"generated_version\": generated_version,\n", + " \"size_bytes\": path.stat().st_size}\n", + "print(f\"Catalog source: {CATALOG_SOURCE}; IDC release: {client.get_idc_version()}\")" + ] + }, + { + "cell_type": "markdown", + "id": "2a5f05b7", + "metadata": {}, + "source": [ + "### 3. Select four referenced pairs and inspect the download estimate\n", + "\n", + "Restrict both modalities to this collection and the documented CC BY 3.0 license.\n", + "Require one SEG instance per series, retain one pair per patient, and show the\n", + "estimate before any DICOM download. Raise an error if the configured size cap\n", + "or patient count cannot be met. Check the [collection's usage and citation requirements](https://www.cancerimagingarchive.net/collection/c4kc-kits/)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2e65c739", + "metadata": {}, + "outputs": [], + "source": [ + "cohort = index_table.loc[index_table.collection_id.eq(COLLECTION)].copy()\n", + "seg_rows = cohort.loc[cohort.Modality.eq(\"SEG\") & cohort.instanceCount.eq(1)].merge(\n", + " seg_table[[\"SeriesInstanceUID\", \"segmented_SeriesInstanceUID\", \"SegmentationType\"]],\n", + " on=\"SeriesInstanceUID\", validate=\"one_to_one\")\n", + "seg_rows = seg_rows.loc[seg_rows.SegmentationType.eq(\"BINARY\")]\n", + "ct_rows = cohort.loc[cohort.Modality.eq(\"CT\")]\n", + "pairs = seg_rows.merge(ct_rows, left_on=\"segmented_SeriesInstanceUID\", right_on=\"SeriesInstanceUID\",\n", + " suffixes=(\"_seg\", \"_ct\"), validate=\"many_to_one\")\n", + "pairs = pairs.loc[pairs.PatientID_seg.eq(pairs.PatientID_ct)\n", + " & pairs.StudyInstanceUID_seg.eq(pairs.StudyInstanceUID_ct)\n", + " & pairs.license_short_name_seg.eq(\"CC BY 3.0\")\n", + " & pairs.license_short_name_ct.eq(\"CC BY 3.0\")].copy()\n", + "pairs[\"estimated_mb\"] = pairs.series_size_MB_seg + pairs.series_size_MB_ct\n", + "pairs = pairs.sort_values([\"estimated_mb\", \"PatientID_ct\", \"SeriesInstanceUID_seg\"], kind=\"stable\")\n", + "pairs = pairs.drop_duplicates(\"PatientID_ct\").head(N_PATIENTS)\n", + "require(N_PATIENTS >= 2 and len(pairs) == N_PATIENTS, \"Not enough eligible patient pairs\")\n", + "require(np.isfinite(pairs.estimated_mb).all() and (pairs.estimated_mb > 0).all(), \"Invalid size estimates\")\n", + "estimated_mb = float(pairs.estimated_mb.sum())\n", + "require(estimated_mb <= MAX_DICOM_MB, \"Cohort exceeds MAX_DICOM_MB; reduce N_PATIENTS\")\n", + "selected_uids = pairs.SeriesInstanceUID_ct.tolist() + pairs.SeriesInstanceUID_seg.tolist()\n", + "selected_catalog = cohort.loc[cohort.SeriesInstanceUID.isin(selected_uids)]\n", + "client_rows = client.index.set_index(\"SeriesInstanceUID\")\n", + "for row in selected_catalog.to_dict(\"records\"):\n", + " require(row[\"SeriesInstanceUID\"] in client_rows.index, \"Selected series missing from download client\")\n", + " download_row = client_rows.loc[row[\"SeriesInstanceUID\"]]\n", + " require(str(download_row.crdc_series_uuid) == str(row[\"crdc_series_uuid\"])\n", + " and int(download_row.instanceCount) == int(row[\"instanceCount\"]), \"Catalog/download client mismatch\")\n", + "selected_catalog.to_json(OUTPUT_DIR / \"cohort_catalog.json\", orient=\"records\", indent=2)\n", + "print(f\"DICOM estimate: {estimated_mb:.2f} MB for {len(pairs)} patients; cap {MAX_DICOM_MB:.0f} MB\")\n", + "display(pairs[[\"PatientID_ct\", \"instanceCount_ct\", \"estimated_mb\", \"license_short_name_ct\"]]\n", + " .rename(columns={\"PatientID_ct\": \"Patient\", \"instanceCount_ct\": \"CT slices\", \"estimated_mb\": \"CT + SEG (MB)\", \"license_short_name_ct\": \"License\"}))" + ] + }, + { + "cell_type": "markdown", + "id": "2bb962c5", + "metadata": {}, + "source": [ + "### 4. Download, validate and export one patient at a time\n", + "\n", + "Every file's SHA256 is recorded. The catalog's CT instance count must match the\n", + "downloaded series and the SEG's global source references. All source/segment\n", + "frames must be present by default; zero-filled omitted frames require an explicit\n", + "dataset assumption, controlled by `ALLOW_OMITTED_EMPTY_FRAMES`.\n", + "\n", + "The preview uses a fixed HU display window and the slice with the largest total\n", + "segmented area. This changes only the display, not stored HU values. No patient\n", + "identifiers are drawn into the image. Each exported mask channel has its original\n", + "segment metadata; DICOM algorithm type is preserved without reinterpretation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "528ac35b", + "metadata": {}, + "outputs": [], + "source": [ + "records = []\n", + "canonical_channel_keys = None\n", + "for row in pairs.to_dict(\"records\"):\n", + " patient = row[\"PatientID_ct\"]\n", + " ct_uid, seg_uid = row[\"SeriesInstanceUID_ct\"], row[\"SeriesInstanceUID_seg\"]\n", + " client.download_from_selection(downloadDir=str(DICOM_DIR), seriesInstanceUID=[ct_uid, seg_uid],\n", + " show_progress_bar=False, use_s5cmd_sync=True)\n", + " study_dir = DICOM_DIR / COLLECTION / patient / row[\"StudyInstanceUID_ct\"]\n", + " ct_paths = sorted((study_dir / f\"CT_{ct_uid}\").glob(\"*.dcm\"))\n", + " seg_paths = sorted((study_dir / f\"SEG_{seg_uid}\").glob(\"*.dcm\"))\n", + " require(len(seg_paths) == 1, \"Expected exactly one downloaded SEG instance\")\n", + " ct_datasets = [pydicom.dcmread(path) for path in ct_paths]\n", + " seg = hd.seg.segread(seg_paths[0])\n", + " require(all(str(ds.SeriesInstanceUID) == ct_uid and str(ds.PatientID) == patient for ds in ct_datasets), \"Downloaded CT differs from selection\")\n", + " require(str(seg.SeriesInstanceUID) == seg_uid, \"Downloaded SEG differs from selection\")\n", + " image_hu, masks, metadata = validate_and_convert(ct_datasets, seg, int(row[\"instanceCount_ct\"]),\n", + " ALLOW_OMITTED_EMPTY_FRAMES)\n", + " channel_keys = [item[\"semantic_key\"] for item in metadata[\"segments\"]]\n", + " if canonical_channel_keys is None:\n", + " canonical_channel_keys = channel_keys\n", + " require(channel_keys == canonical_channel_keys, \"Cases have different semantic channel definitions\")\n", + " sample_id = hashlib.sha256(ct_uid.encode()).hexdigest()[:16]\n", + " tensor_file = f\"{sample_id}.pt\"\n", + " torch.save({\"image_hu\": torch.from_numpy(image_hu), \"masks\": torch.from_numpy(masks)}, OUTPUT_DIR / tensor_file)\n", + " metadata[\"collection_id\"] = COLLECTION\n", + " metadata[\"sources\"] = [{\"path\": str(path.relative_to(DICOM_DIR)), \"sha256\": file_sha256(path),\n", + " \"size_bytes\": path.stat().st_size} for path in ct_paths + seg_paths]\n", + " metadata[\"license_ct\"] = row[\"license_short_name_ct\"]\n", + " metadata[\"license_seg\"] = row[\"license_short_name_seg\"]\n", + " metadata[\"source_dois\"] = sorted({row[\"source_DOI_ct\"], row[\"source_DOI_seg\"]})\n", + " metadata_file = f\"{sample_id}.json\"\n", + " (OUTPUT_DIR / metadata_file).write_text(json.dumps(metadata, indent=2))\n", + " records.append({\"collection_id\": COLLECTION, \"patient_id\": patient,\n", + " \"ct_series_uid\": ct_uid, \"seg_series_uid\": seg_uid,\n", + " \"tensor_file\": tensor_file, \"metadata_file\": metadata_file,\n", + " \"ct_slices\": image_hu.shape[1], \"image_shape\": list(image_hu.shape),\n", + " \"mask_shape\": list(masks.shape), \"nonzero_voxels\": metadata[\"nonzero_voxels_by_channel\"],\n", + " \"dicom_bytes\": sum(item[\"size_bytes\"] for item in metadata[\"sources\"])})\n", + " slice_index = int(np.argmax(masks.sum(axis=(0, 2, 3))))\n", + " fig, axes = plt.subplots(1, len(masks) + 1, figsize=(5 * (len(masks) + 1), 4))\n", + " axes[0].imshow(image_hu[0, slice_index], cmap=\"gray\", vmin=-200, vmax=300)\n", + " axes[0].set_title(f\"CT (HU); slice {slice_index}\")\n", + " for axis, mask, segment in zip(axes[1:], masks, metadata[\"segments\"]):\n", + " axis.imshow(image_hu[0, slice_index], cmap=\"gray\", vmin=-200, vmax=300)\n", + " axis.imshow(np.ma.masked_equal(mask[slice_index], 0), cmap=\"autumn\", alpha=0.45, vmin=0, vmax=1)\n", + " axis.set_title(f\"{segment['label']} — {segment['algorithm_type']}\")\n", + " for axis in axes:\n", + " axis.axis(\"off\")\n", + " fig.suptitle(f\"Validated CT/SEG pair · {sample_id}\")\n", + " fig.tight_layout()\n", + " fig.savefig(OUTPUT_DIR / f\"{sample_id}_overlay.png\", dpi=120)\n", + " plt.show()\n", + " plt.close(fig)\n", + " print(f\"Validated {patient}: {image_hu.shape} image; {masks.shape} masks; omitted frames {metadata['checks']['omitted_source_segment_frames']}\")\n", + " del ct_datasets, seg, image_hu, masks" + ] + }, + { + "cell_type": "markdown", + "id": "19a5fc07", + "metadata": {}, + "source": [ + "### 5. Write and verify the patient-disjoint split\n", + "\n", + "Grouping uses `(collection_id, PatientID)`, since PatientID is scoped to a\n", + "collection. The fixed seed is applied to a sorted patient list, so input row\n", + "order does not affect assignments. All cases for a patient share one partition.\n", + "This grouping cannot identify aliases or duplicate people across collections." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "280de258", + "metadata": {}, + "outputs": [], + "source": [ + "split_records = patient_disjoint_split(records, SPLIT_SEED, TEST_FRACTION)\n", + "manifest = {\n", + " \"idc_release\": client.get_idc_version(), \"catalog_source\": CATALOG_SOURCE,\n", + " \"hf_revision\": HF_REVISION if CATALOG_SOURCE == \"huggingface\" else None,\n", + " \"catalog_provenance\": catalog_provenance,\n", + " \"packages\": {name: version(name) for name in (\"idc-index\", \"idc-index-data\", \"highdicom\", \"pydicom\", \"numpy\", \"pandas\", \"torch\", \"matplotlib\")},\n", + " \"split_seed\": SPLIT_SEED, \"test_fraction\": TEST_FRACTION,\n", + " \"sampling\": \"smallest referenced CT/SEG pair per patient; tutorial convenience sample\",\n", + " \"allow_omitted_empty_frames\": ALLOW_OMITTED_EMPTY_FRAMES,\n", + " \"channel_semantic_keys\": canonical_channel_keys,\n", + " \"estimated_dicom_mb\": estimated_mb,\n", + " \"actual_dicom_bytes\": sum(record[\"dicom_bytes\"] for record in records),\n", + " \"runtime_seconds\": round(perf_counter() - run_started, 2),\n", + " \"records\": split_records,\n", + "}\n", + "(OUTPUT_DIR / \"split_manifest.json\").write_text(json.dumps(manifest, indent=2))\n", + "display(pd.DataFrame(split_records)[[\"patient_id\", \"split\", \"ct_slices\", \"image_shape\", \"mask_shape\", \"nonzero_voxels\"]])\n", + "print(f\"Actual selected DICOM: {manifest['actual_dicom_bytes'] / 1e6:.2f} MB; runtime {manifest['runtime_seconds']:.2f} seconds\")" + ] + }, + { + "cell_type": "markdown", + "id": "42271c87", + "metadata": {}, + "source": [ + "### 6. Consume saved tensors with PyTorch\n", + "\n", + "This Dataset opens one case at a time. Batch size 1 accommodates different\n", + "volume depths. The optional fixed window `[-200, 300] HU -> [0, 1]` is an\n", + "illustrative input transform, not a recommended training configuration or a\n", + "transformation fitted on the test set. Stored tensors always retain HU.\n", + "Only load `.pt` files generated by this notebook, using `weights_only=True`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "49a882e8", + "metadata": {}, + "outputs": [], + "source": [ + "class CTSegDataset(torch.utils.data.Dataset):\n", + " def __init__(self, output_dir, partition, normalize=False):\n", + " require(partition in (\"train\", \"test\"), \"Unknown dataset partition\")\n", + " self.output_dir = Path(output_dir)\n", + " saved_manifest = json.loads((self.output_dir / \"split_manifest.json\").read_text())\n", + " self.records = [record for record in saved_manifest[\"records\"] if record[\"split\"] == partition]\n", + " self.normalize = normalize\n", + "\n", + " def __len__(self):\n", + " return len(self.records)\n", + "\n", + " def __getitem__(self, index):\n", + " record = self.records[index]\n", + " tensors = torch.load(self.output_dir / record[\"tensor_file\"], map_location=\"cpu\", weights_only=True)\n", + " image = tensors[\"image_hu\"]\n", + " if self.normalize:\n", + " image = (image.clamp(-200, 300) + 200) / 500\n", + " return {\"image\": image, \"masks\": tensors[\"masks\"], \"patient_id\": record[\"patient_id\"]}\n", + "\n", + "train_dataset = CTSegDataset(OUTPUT_DIR, \"train\", normalize=True)\n", + "test_dataset = CTSegDataset(OUTPUT_DIR, \"test\", normalize=True)\n", + "loader = torch.utils.data.DataLoader(train_dataset, batch_size=1, shuffle=False, num_workers=0)\n", + "batch = next(iter(loader))\n", + "print(\"Train/test patients:\", len(train_dataset), len(test_dataset))\n", + "print(\"Batch shapes:\", tuple(batch[\"image\"].shape), tuple(batch[\"masks\"].shape))" + ] + }, + { + "cell_type": "markdown", + "id": "e23bf602", + "metadata": {}, + "source": [ + "## Checks\n", + "\n", + "Reload saved files and verify shapes, binary values, coded channel consistency,\n", + "source count, and patient separation. These checks demonstrate data preparation;\n", + "they do not establish annotation quality or predict segmentation accuracy." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4fb10c3c", + "metadata": {}, + "outputs": [], + "source": [ + "saved = json.loads((OUTPUT_DIR / \"split_manifest.json\").read_text())\n", + "train_groups = {(r[\"collection_id\"], r[\"patient_id\"]) for r in saved[\"records\"] if r[\"split\"] == \"train\"}\n", + "test_groups = {(r[\"collection_id\"], r[\"patient_id\"]) for r in saved[\"records\"] if r[\"split\"] == \"test\"}\n", + "require(train_groups and test_groups and train_groups.isdisjoint(test_groups), \"Saved split leaks patients\")\n", + "for record in saved[\"records\"]:\n", + " tensors = torch.load(OUTPUT_DIR / record[\"tensor_file\"], map_location=\"cpu\", weights_only=True)\n", + " metadata = json.loads((OUTPUT_DIR / record[\"metadata_file\"]).read_text())\n", + " require(list(tensors[\"image_hu\"].shape) == record[\"image_shape\"], \"Saved image shape mismatch\")\n", + " require(list(tensors[\"masks\"].shape) == record[\"mask_shape\"], \"Saved mask shape mismatch\")\n", + " require(tensors[\"image_hu\"].dtype == torch.float32 and torch.isfinite(tensors[\"image_hu\"]).all().item(), \"Invalid HU tensor\")\n", + " require(tensors[\"masks\"].dtype == torch.uint8 and torch.all((tensors[\"masks\"] == 0) | (tensors[\"masks\"] == 1)).item(), \"Invalid binary masks\")\n", + " require([segment[\"semantic_key\"] for segment in metadata[\"segments\"]] == saved[\"channel_semantic_keys\"], \"Saved channels differ between cases\")\n", + " require(len(metadata[\"source_sop_instance_uids_in_slice_order\"]) == record[\"ct_slices\"], \"Saved source count mismatch\")\n", + " del tensors\n", + "print(\"PASS: validated real CT/SEG pairs, consistent tensor channels, and disjoint patient groups\")" + ] + }, + { + "cell_type": "markdown", + "id": "3d0ac3a8", + "metadata": {}, + "source": [ + "## Next Steps\n", + "\n", + "Use a larger, representative cohort before training or assessing a model.\n", + "Keep patient grouping when adding studies, slices or patches. Fit any learned\n", + "preprocessing on training patients only. Validate resampling and class semantics\n", + "explicitly before extending this example to other collections or modalities.\n", + "\n", + "Outputs are under `idc_ct_seg_data/outputs`: `.pt` tensors, per-case JSON geometry\n", + "and source checksums, overlay PNGs, the selected catalog rows, and the split\n", + "manifest. The four-patient split is an educational artifact, not a benchmark.\n", + "\n", + "## Support\n", + "\n", + "For questions, use the [IDC user forum](https://discourse.canceridc.dev) or open\n", + "an issue in [IDC-Tutorials](https://github.com/ImagingDataCommons/IDC-Tutorials/issues).\n", + "More examples are available in the [IDC-Tutorials repository](https://github.com/ImagingDataCommons/IDC-Tutorials).\n", + "\n", + "## Acknowledgments\n", + "\n", + "Imaging Data Commons has been funded in whole or in part with Federal funds from\n", + "the National Cancer Institute, National Institutes of Health, under Task Order\n", + "No. HHSN26110071 under Contract No. HHSN261201500003I.\n", + "\n", + "If you use IDC in research, cite Fedorov et al., *National Cancer Institute Imaging\n", + "Data Commons: Toward Transparency, Reproducibility, and Scalability in Imaging\n", + "Artificial Intelligence*, RadioGraphics (2023), [doi:10.1148/rg.230180](https://doi.org/10.1148/rg.230180).\n", + "\n", + "The C4KC-KiTS data are provided under CC BY 3.0. Cite Heller et al., *Data from\n", + "C4KC-KiTS*, The Cancer Imaging Archive (2019), [doi:10.7937/TCIA.2019.IX49E8NX](https://doi.org/10.7937/TCIA.2019.IX49E8NX),\n", + "and follow the [collection's full citation instructions](https://www.cancerimagingarchive.net/collection/c4kc-kits/),\n", + "including Heller et al., *The KiTS19 Challenge Data: 300 Kidney Tumor Cases with\n", + "Clinical Context, CT Semantic Segmentations, and Surgical Outcomes*,\n", + "[arXiv:1904.00445](https://arxiv.org/abs/1904.00445).\n", + "\n", + "Libraries used: [idc-index](https://github.com/ImagingDataCommons/idc-index),\n", + "[highdicom](https://highdicom.readthedocs.io/en/latest/seg.html),\n", + "[pydicom](https://pydicom.github.io/), and [PyTorch](https://pytorch.org/)." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/test/src/test_ct_seg_to_tensors.py b/test/src/test_ct_seg_to_tensors.py new file mode 100644 index 0000000..293dbe6 --- /dev/null +++ b/test/src/test_ct_seg_to_tensors.py @@ -0,0 +1,184 @@ +"""Behavior checks for the self-contained CT/SEG tutorial definitions. + +Run after installing the notebook's dependencies: + python -m pytest test/src/test_ct_seg_to_tensors.py +""" +import copy +from pathlib import Path + +import highdicom as hd +import nbformat +import numpy as np +import pydicom +import pytest +from pydicom.dataset import Dataset, FileMetaDataset +from pydicom.uid import CTImageStorage, ExplicitVRLittleEndian, generate_uid + + +NOTEBOOK = Path(__file__).resolve().parents[2] / "notebooks/advanced_topics/ct_seg_to_tensors.ipynb" +namespace = {} +for cell in nbformat.read(NOTEBOOK, as_version=4).cells: + if "definitions" in cell.metadata.get("tags", []): + exec(compile(cell.source, str(NOTEBOOK), "exec"), namespace) +convert = namespace["validate_and_convert"] +split_patients = namespace["patient_disjoint_split"] +ValidationError = namespace["ValidationError"] + + +@pytest.fixture +def pair(): + study, series, frame_of_reference = generate_uid(), generate_uid(), generate_uid() + images = [] + for z in (0, 5, 10): + image = Dataset() + image.file_meta = FileMetaDataset() + image.file_meta.TransferSyntaxUID = ExplicitVRLittleEndian + image.SOPClassUID = CTImageStorage + image.SOPInstanceUID = generate_uid() + image.PatientID = "tutorial-patient" + image.PatientName = "Synthetic^Patient" + image.PatientBirthDate = "" + image.PatientSex = "" + image.StudyInstanceUID = study + image.SeriesInstanceUID = series + image.FrameOfReferenceUID = frame_of_reference + image.StudyID = "1" + image.SeriesNumber = 1 + image.InstanceNumber = z + 1 + image.StudyDate = "20260101" + image.StudyTime = "120000" + image.AccessionNumber = "" + image.Modality = "CT" + image.ImageOrientationPatient = [1, 0, 0, 0, 1, 0] + image.ImagePositionPatient = [0, 0, z] + image.PixelSpacing = [2, 3] + image.SliceThickness = 5 + image.Rows, image.Columns = 4, 4 + image.SamplesPerPixel = 1 + image.PhotometricInterpretation = "MONOCHROME2" + image.BitsAllocated = 16 + image.BitsStored = 16 + image.HighBit = 15 + image.PixelRepresentation = 0 + image.RescaleSlope = 2 + image.RescaleIntercept = -1000 + image.RescaleType = "HU" + image.PixelData = np.full((4, 4), z + 1, dtype=np.uint16).tobytes() + images.append(image) + masks = np.zeros((3, 4, 4, 2), dtype=np.uint8) + masks[0, 1, 2, 0] = 1 + masks[2, 2, 1, 0] = 1 + masks[2, 2, 1, 1] = 1 # Intentional overlap must survive conversion. + descriptions = [ + hd.seg.SegmentDescription( + segment_number=number, segment_label=label, + segmented_property_category=pydicom.sr.codedict.codes.SCT.AnatomicalStructure, + segmented_property_type=code, + algorithm_type=hd.seg.SegmentAlgorithmTypeValues.MANUAL, + ) + for number, label, code in [ + (1, "Kidney", pydicom.sr.codedict.codes.SCT.Kidney), + (2, "Liver", pydicom.sr.codedict.codes.SCT.Liver), + ] + ] + seg = hd.seg.Segmentation( + source_images=images, pixel_array=masks, + segmentation_type=hd.seg.SegmentationTypeValues.BINARY, + segment_descriptions=descriptions, omit_empty_frames=False, + series_instance_uid=generate_uid(), series_number=2, + sop_instance_uid=generate_uid(), instance_number=1, + manufacturer="Tutorial", manufacturer_model_name="Synthetic", + software_versions="1", device_serial_number="1", + ) + return images, copy.deepcopy(seg) + + +def test_ct_order_hu_affine_and_overlapping_mask_channels(pair): + images, seg = pair + image, masks, metadata = convert([images[2], images[0], images[1]], seg, 3) + assert image.shape == (1, 3, 4, 4) + assert image.dtype == np.float32 + assert image[0, :, 0, 0].tolist() == [-998, -988, -978] + np.testing.assert_allclose(metadata["affine_dhw_to_lps_mm"], + [[0, 0, 3, 0], [0, 2, 0, 0], [5, 0, 0, 0], [0, 0, 0, 1]]) + assert masks.dtype == np.uint8 + assert masks[:, 2, 2, 1].tolist() == [1, 1] + assert sorted(np.count_nonzero(mask) for mask in masks) == [1, 2] + + +@pytest.mark.parametrize("change, message", [ + ("patient", "PatientID mismatch"), + ("series", "different CT series"), + ("frame_of_reference", "FrameOfReferenceUID mismatch"), + ("position", "SEG position differs"), + ("spacing", "SEG pixel spacing differs"), + ("orientation", "SEG orientation differs"), + ("spatial_no", "does not preserve"), +]) +def test_incompatible_pair_is_rejected(pair, change, message): + images, seg = pair + if change == "patient": + seg.PatientID = "wrong-patient" + elif change == "series": + seg.ReferencedSeriesSequence[0].SeriesInstanceUID = generate_uid() + elif change == "frame_of_reference": + seg.FrameOfReferenceUID = generate_uid() + elif change == "position": + seg.PerFrameFunctionalGroupsSequence[0].PlanePositionSequence[0].ImagePositionPatient = [1, 0, 0] + elif change == "spacing": + seg.SharedFunctionalGroupsSequence[0].PixelMeasuresSequence[0].PixelSpacing = [1, 1] + elif change == "orientation": + seg.SharedFunctionalGroupsSequence[0].PlaneOrientationSequence[0].ImageOrientationPatient = [0, 1, 0, 1, 0, 0] + else: + seg.PerFrameFunctionalGroupsSequence[0].DerivationImageSequence[0].SourceImageSequence[0].SpatialLocationsPreserved = "NO" + with pytest.raises(ValidationError, match=message): + convert(images, seg, 3) + + +@pytest.mark.parametrize("change, message", [ + ("missing", "instance count"), + ("duplicate_uid", "Duplicate CT SOP"), + ("irregular", "Irregular CT slice spacing"), + ("gantry_tilt", "Gantry tilt"), +]) +def test_incomplete_or_unsupported_ct_is_rejected(pair, change, message): + images, seg = pair + if change == "missing": + images = images[:-1] + elif change == "duplicate_uid": + images[1].SOPInstanceUID = images[0].SOPInstanceUID + elif change == "irregular": + images[2].ImagePositionPatient = [0, 0, 11] + else: + images[2].ImagePositionPatient = [1, 0, 10] + with pytest.raises(ValidationError, match=message): + convert(images, seg, 3) + + +def test_missing_spatial_flag_is_allowed_only_after_geometry_validation(pair): + images, seg = pair + for frame in seg.PerFrameFunctionalGroupsSequence: + source = frame.DerivationImageSequence[0].SourceImageSequence[0] + if "SpatialLocationsPreserved" in source: + del source.SpatialLocationsPreserved + _, masks, metadata = convert(images, seg, 3) + assert metadata["checks"]["spatial_locations_flag_absent"] is True + assert masks[:, 2, 2, 1].tolist() == [1, 1] + + +def test_patient_split_is_grouped_and_independent_of_record_order(): + records = [{"collection_id": "c4kc_kits", "patient_id": patient, "study": study} + for patient, study in [("A", "1"), ("A", "2"), ("B", "3"), ("C", "4"), ("D", "5")]] + first = split_patients(records, seed=2026) + second = split_patients(list(reversed(records)), seed=2026) + assignment = {record["patient_id"]: record["split"] for record in first} + assert assignment == {record["patient_id"]: record["split"] for record in second} + assert len({record["split"] for record in first if record["patient_id"] == "A"}) == 1 + assert sum(partition == "test" for partition in assignment.values()) == 1 + assert {record["patient_id"] for record in first if record["split"] == "train"}.isdisjoint( + {record["patient_id"] for record in first if record["split"] == "test"}) + + +def test_one_patient_cannot_form_train_and_test(): + with pytest.raises(ValidationError, match="At least two patient groups"): + split_patients([{"collection_id": "c4kc_kits", "patient_id": "A"}]) From c57646a0d433d26296b5a635dc272adedd23e484 Mon Sep 17 00:00:00 2001 From: Denis Samatov <98321287+denis-samatov@users.noreply.github.com> Date: Thu, 15 Oct 2026 12:22:41 +0700 Subject: [PATCH 2/5] Use IDC geometry preflight and highdicom volumes for CT/SEG tensors --- README.md | 2 +- .../advanced_topics/ct_seg_to_tensors.ipynb | 233 +++++++++++------- test/src/test_ct_seg_to_tensors.py | 125 +++++++++- 3 files changed, 260 insertions(+), 100 deletions(-) diff --git a/README.md b/README.md index fa88d90..b79ab68 100644 --- a/README.md +++ b/README.md @@ -50,7 +50,7 @@ Zenodo permanent citeable record: https://zenodo.org/records/22691908. Notebooks in this folder focus on topics that will require understanding of the basics, and aim to address more narrow use cases of IDC usage. * [Searching DICOM private tags](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/dicom_private_tags_intro.ipynb): all of DICOM attributes for the imaging data in IDC are searchable using BigQuery. DICOM private tags often contain critical information, such as diffusion b-values, but are a bit more tricky to access from BigQuery. In this tutorial you will learn how to accomplish this. -* [Validated CT/DICOM SEG tensors and patient-disjoint splits](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/ct_seg_to_tensors.ipynb): select a small C4KC-KiTS cohort from the pinned IDC catalog, validate source references and spatial correspondence, export HU images and separate mask channels for PyTorch, and create a reproducible split by patient. +* [Validated CT/DICOM SEG tensors and patient-disjoint splits](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/ct_seg_to_tensors.ipynb): select referenced regular-grid C4KC-KiTS pairs with seg_index and volume_geometry_index, load HU volumes and separate mask channels with highdicom, validate source correspondence, and create a reproducible patient-disjoint split for PyTorch. * [Working with IDC clinical data](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/clinical_data_intro.ipynb): this tutorial will help you become familiar with how clinical and other tabular data accompanying images in IDC is organized, and how you can use it to build cohorts. * [Direct access to tags/frames from GCS/AWS buckets](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/gcs_aws_direct_access.ipynb): learn how to access individual frames or tags of large DICOM files from the bucket without having to download the entire file (this notebook accompanies documentation article here: https://learn.canceridc.dev/data/downloading-data/direct-loading) * [Using DICOMweb to access IDC data](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/idc_dicomweb_access.ipynb): both metadata and pixel data can be access using DICOMweb, which is particularly important while working with digital pathology, as it enables granular access to the individual pyramid tiles (frames) diff --git a/notebooks/advanced_topics/ct_seg_to_tensors.ipynb b/notebooks/advanced_topics/ct_seg_to_tensors.ipynb index 7e347c2..bfd3af3 100644 --- a/notebooks/advanced_topics/ct_seg_to_tensors.ipynb +++ b/notebooks/advanced_topics/ct_seg_to_tensors.ipynb @@ -31,8 +31,11 @@ "- Legacy single-frame CT on a regular grid, linear rescaling to HU, and one BINARY\n", " SEG instance per series. Enhanced CT, non-regular spacing, gantry tilt, and\n", " resampled masks are rejected. Nothing is silently resampled.\n", - "- CT is selected by the SEG's **referenced series**, then verified against source\n", - " SOP Instance UIDs and per-frame geometry. A patient/study match is insufficient.\n", + "- `seg_index.segmented_SeriesInstanceUID` selects the CT, and\n", + " `volume_geometry_index.regularly_spaced_3d_volume` filters it before download.\n", + " `highdicom` loads HU volumes and aligns separate SEG channels without resampling.\n", + " Source SOP Instance UIDs and SEG-to-source frame correspondence remain local\n", + " integrity checks; a patient/study match is insufficient.\n", "- Separate channels retain overlaps. Coded category/type, modifiers and anatomic\n", " region define a stable channel order across cases; free-text labels and local\n", " segment numbers alone do not define classes.\n", @@ -62,8 +65,19 @@ "metadata": {}, "outputs": [], "source": [ - "%pip install -q \"idc-index==0.13.0\" \"pydicom==3.0.2\" \"highdicom==0.28.2\" \"numpy==2.3.5\" \"pandas==2.2.3\" \"matplotlib==3.10.8\"\n", - "%pip install -q \"torch==2.14.1\" --index-url https://download.pytorch.org/whl/cpu" + "import subprocess\n", + "import sys\n", + "\n", + "# check_call stops execution if either installation fails.\n", + "subprocess.check_call([\n", + " sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"idc-index==0.13.0\",\n", + " \"pydicom==3.0.2\", \"highdicom==0.28.2\", \"numpy==2.3.5\",\n", + " \"pandas==2.2.3\", \"matplotlib==3.10.8\",\n", + "])\n", + "subprocess.check_call([\n", + " sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"torch==2.14.1\",\n", + " \"--index-url\", \"https://download.pytorch.org/whl/cpu\",\n", + "])" ] }, { @@ -125,17 +139,34 @@ "source": [ "## Steps\n", "\n", - "### 1. Define the validation boundary\n", - "\n", - "The following functions validate DICOM objects before turning them into arrays.\n", - "Geometric tolerances are 0.01 mm for positions/slice spacing and `1e-4` for\n", - "direction cosines/in-plane spacing. The affine maps array indices `[D, H, W]`\n", - "to **LPS millimeters**, not NIfTI RAS coordinates.\n", - "\n", - "C4KC-KiTS SEG objects do not set `SpatialLocationsPreserved`. The code tolerates\n", - "that absent flag only after independently checking every frame's geometry.\n", - "An explicit `NO` is always rejected. `SliceThickness` is not used as the slice\n", - "step: positions determine it." + "### 1. Use the catalog and highdicom at their intended boundaries\n", + "\n", + "The [IDC index reference](https://idc-index.readthedocs.io/en/latest/indices_reference.html)\n", + "and [tensor-loading example](https://huggingface.co/datasets/ImagingDataCommons/idc-index-data#loading-images-as-tensors)\n", + "provide the building blocks used below:\n", + "\n", + "- `seg_index` supplies CT/SEG correspondence, without reconstructing it manually.\n", + "- `volume_geometry_index` excludes CT series that are not regular 3D grids\n", + " before any DICOM download.\n", + "- `hd.get_volume_from_series()` orders CT slices, checks the actual downloaded\n", + " geometry, applies the modality transform to HU, and retains the LPS affine.\n", + "- `seg.get_volume(combine_segments=False).match_geometry(ct_volume)` aligns\n", + " separate binary channels by padding/flipping/permuting; it does not interpolate.\n", + "\n", + "The remaining local guards check download completeness, declared HU units,\n", + "source UIDs, and each SEG frame's correspondence to its referenced CT. The CT\n", + "geometry index does not establish those pair-specific properties. C4KC-KiTS\n", + "omits `SpatialLocationsPreserved`; absence is accepted only after these checks,\n", + "while explicit `NO` is rejected. SEG position tolerance is 0.01 mm and orientation\n", + "and pixel-spacing tolerance is `1e-4`. Highdicom CT spacing tolerance is 0.01 mm;\n", + "the final grid match uses `1e-4` voxels to tolerate rounding only.\n", + "\n", + "The exported affine maps the library's array indices `[D, H, W]` to **DICOM LPS\n", + "millimeters**. Slice direction comes from highdicom and may differ from file or\n", + "InstanceNumber order. No display window or real-world-value transform is applied\n", + "to the saved HU image. HU units must be explicit or guaranteed by the\n", + "[original CT rule in DICOM PS3.3 C.8.2.1](https://dicom.nema.org/medical/Dicom/current/output/chtml/part03/sect_C.8.2.html);\n", + "missing units on derived/localizer/multi-energy images are not assumed to be HU." ] }, { @@ -165,6 +196,7 @@ "import random\n", "from pathlib import Path\n", "\n", + "import pandas as pd\n", "import highdicom as hd\n", "import numpy as np\n", "import pydicom\n", @@ -191,7 +223,7 @@ "id": "b6631e97", "metadata": {}, "source": [ - "### Validate and order the CT volume" + "### Load the CT volume with highdicom" ] }, { @@ -205,49 +237,75 @@ }, "outputs": [], "source": [ - "def sort_and_validate_ct(datasets, expected_count: int):\n", - " \"\"\"Return geometrically sorted single-frame CT instances and an LPS affine.\"\"\"\n", + "def load_ct_volume(datasets, expected_count: int):\n", + " \"\"\"Check identity/completeness and delegate volume geometry and HU to highdicom.\"\"\"\n", " require(len(datasets) == expected_count, \"CT instance count differs from the catalog\")\n", " require(len(datasets) >= 2, \"A volume needs at least two CT slices\")\n", " first = datasets[0]\n", - " uids = [str(d.SOPInstanceUID) for d in datasets]\n", + " uids = [str(ds.SOPInstanceUID) for ds in datasets]\n", " require(len(set(uids)) == len(uids), \"Duplicate CT SOP Instance UID\")\n", - " orientation = numeric_vector(first.ImageOrientationPatient, 6, \"CT orientation\")\n", - " row_direction, column_direction = orientation[:3], orientation[3:]\n", - " require(np.isclose(np.linalg.norm(row_direction), 1, atol=1e-4), \"Non-unit CT row direction\")\n", - " require(np.isclose(np.linalg.norm(column_direction), 1, atol=1e-4), \"Non-unit CT column direction\")\n", - " require(abs(np.dot(row_direction, column_direction)) < 1e-4, \"Non-orthogonal CT directions\")\n", - " normal = np.cross(row_direction, column_direction)\n", - " normal /= np.linalg.norm(normal)\n", - " spacing = numeric_vector(first.PixelSpacing, 2, \"CT pixel spacing\")\n", - " require((spacing > 0).all(), \"CT pixel spacing must be positive\")\n", - " positions = []\n", " for ds in datasets:\n", " require(str(ds.SOPClassUID) == str(CTImageStorage), \"Only legacy single-frame CT is supported\")\n", " require(int(ds.get(\"NumberOfFrames\", 1)) == 1, \"Enhanced/multiframe CT is unsupported\")\n", " require(ds.Modality == \"CT\", \"Non-CT instance in the CT series\")\n", " for key in (\"PatientID\", \"StudyInstanceUID\", \"SeriesInstanceUID\", \"FrameOfReferenceUID\"):\n", " require(str(ds.get(key, \"\")) == str(first.get(key, \"\")) != \"\", f\"Inconsistent CT {key}\")\n", - " require((ds.Rows, ds.Columns) == (first.Rows, first.Columns), \"Inconsistent CT dimensions\")\n", - " require(np.allclose(numeric_vector(ds.ImageOrientationPatient, 6, \"CT orientation\"), orientation, rtol=0, atol=1e-4), \"CT orientation varies\")\n", - " require(np.allclose(numeric_vector(ds.PixelSpacing, 2, \"CT pixel spacing\"), spacing, rtol=0, atol=1e-4), \"CT pixel spacing varies\")\n", - " positions.append(numeric_vector(ds.ImagePositionPatient, 3, \"CT position\"))\n", - " positions = np.stack(positions)\n", - " order = np.argsort(positions @ normal, kind=\"stable\")\n", - " positions = positions[order]\n", - " steps = np.diff(positions @ normal)\n", - " require((steps > 0.01).all(), \"Duplicate or nearly duplicate CT slice positions\")\n", - " slice_spacing = float(np.median(steps))\n", - " require(np.allclose(steps, slice_spacing, rtol=0, atol=0.01), \"Irregular CT slice spacing or missing slices\")\n", - " expected_positions = positions[0] + np.arange(len(datasets))[:, None] * slice_spacing * normal\n", - " require(np.allclose(positions, expected_positions, rtol=0, atol=0.01), \"Gantry tilt or non-regular CT grid is unsupported\")\n", - " affine = np.eye(4)\n", - " # Array indices [slice, row, column] -> DICOM patient coordinates [L, P, S], mm.\n", - " affine[:3, 0] = normal * slice_spacing\n", - " affine[:3, 1] = column_direction * spacing[0]\n", - " affine[:3, 2] = row_direction * spacing[1]\n", - " affine[:3, 3] = positions[0]\n", - " return [datasets[int(i)] for i in order], affine" + " require(\"ModalityLUTSequence\" not in ds, \"This example supports linear CT rescaling only\")\n", + " # DICOM PS3.3 C.8.2.1: original non-localizer, non-multi-energy CT is HU,\n", + " # and RescaleType is optional in that case. Do not infer units for derived CT.\n", + " image_type = list(ds.get(\"ImageType\", []))\n", + " implicit_hu = (len(image_type) >= 3 and image_type[0] == \"ORIGINAL\"\n", + " and image_type[2] != \"LOCALIZER\"\n", + " and str(ds.get(\"MultienergyCTAcquisition\", \"NO\")) == \"NO\")\n", + " units = str(ds.get(\"RescaleType\", \"HU\" if implicit_hu else \"\")).upper()\n", + " require(units == \"HU\", \"CT rescale units cannot be verified as HU\")\n", + " require(\"RescaleSlope\" in ds and \"RescaleIntercept\" in ds, \"CT rescale parameters are missing\")\n", + " slope, intercept = float(ds.RescaleSlope), float(ds.RescaleIntercept)\n", + " require(math.isfinite(slope) and math.isfinite(intercept) and slope != 0, \"Invalid CT rescale parameters\")\n", + " try:\n", + " volume = hd.get_volume_from_series(\n", + " datasets, dtype=np.float32, apply_real_world_transform=False,\n", + " apply_modality_transform=True, apply_voi_transform=False,\n", + " apply_presentation_lut=False, atol=0.01,\n", + " )\n", + " except (ValueError, RuntimeError) as error:\n", + " raise ValidationError(f\"highdicom CT volume validation failed: {error}\") from error\n", + " require(np.isfinite(volume.array).all(), \"Non-finite CT intensities after rescaling\")\n", + " # Use the library's affine to record provenance in the exported slice order.\n", + " positions = np.asarray([ds.ImagePositionPatient for ds in datasets], dtype=np.float64)\n", + " slice_indices = volume.map_reference_to_indices(positions)[:, 0]\n", + " ordered = [datasets[int(i)] for i in np.argsort(slice_indices, kind=\"stable\")]\n", + " return volume, ordered\n", + "\n", + "\n", + "def select_catalog_pairs(index_table, seg_table, geometry_table, collection, n_patients, max_dicom_mb):\n", + " \"\"\"Select referenced regular CT/SEG pairs before downloading any DICOM.\"\"\"\n", + " require(isinstance(n_patients, int) and n_patients >= 2, \"At least two patient pairs are required\")\n", + " require(math.isfinite(max_dicom_mb) and max_dicom_mb > 0, \"Invalid DICOM size cap\")\n", + " cohort = index_table.loc[index_table.collection_id.eq(collection)].copy()\n", + " seg_rows = cohort.loc[cohort.Modality.eq(\"SEG\") & cohort.instanceCount.eq(1)].merge(\n", + " seg_table[[\"SeriesInstanceUID\", \"segmented_SeriesInstanceUID\", \"SegmentationType\"]],\n", + " on=\"SeriesInstanceUID\", validate=\"one_to_one\")\n", + " seg_rows = seg_rows.loc[seg_rows.SegmentationType.eq(\"BINARY\")]\n", + " eligible_geometry = geometry_table.loc[\n", + " geometry_table.regularly_spaced_3d_volume.eq(True).fillna(False),\n", + " [\"SeriesInstanceUID\", \"regularly_spaced_3d_volume\"],\n", + " ]\n", + " ct_rows = cohort.loc[cohort.Modality.eq(\"CT\")].merge(\n", + " eligible_geometry, on=\"SeriesInstanceUID\", validate=\"one_to_one\")\n", + " pairs = seg_rows.merge(ct_rows, left_on=\"segmented_SeriesInstanceUID\", right_on=\"SeriesInstanceUID\",\n", + " suffixes=(\"_seg\", \"_ct\"), validate=\"many_to_one\")\n", + " pairs = pairs.loc[pairs.PatientID_seg.eq(pairs.PatientID_ct)\n", + " & pairs.StudyInstanceUID_seg.eq(pairs.StudyInstanceUID_ct)\n", + " & pairs.license_short_name_seg.eq(\"CC BY 3.0\")\n", + " & pairs.license_short_name_ct.eq(\"CC BY 3.0\")].copy()\n", + " pairs[\"estimated_mb\"] = pairs.series_size_MB_seg + pairs.series_size_MB_ct\n", + " pairs = pairs.sort_values([\"estimated_mb\", \"PatientID_ct\", \"SeriesInstanceUID_seg\"], kind=\"stable\")\n", + " pairs = pairs.drop_duplicates(\"PatientID_ct\").head(n_patients)\n", + " require(len(pairs) == n_patients, \"Not enough eligible patient pairs with regular CT geometry\")\n", + " require(np.isfinite(pairs.estimated_mb).all() and (pairs.estimated_mb > 0).all(), \"Invalid size estimates\")\n", + " require(float(pairs.estimated_mb.sum()) <= max_dicom_mb, \"Cohort exceeds MAX_DICOM_MB; reduce N_PATIENTS\")\n", + " return pairs" ] }, { @@ -412,32 +470,28 @@ "source": [ "def validate_and_convert(ct_datasets, seg, expected_count: int, allow_omitted_empty_frames: bool = False):\n", " \"\"\"Validate one real DICOM pair; return HU image, separate masks, and provenance.\"\"\"\n", - " ct_datasets, affine = sort_and_validate_ct(ct_datasets, expected_count)\n", + " ct_volume, ct_datasets = load_ct_volume(ct_datasets, expected_count)\n", " checks = validate_seg_geometry(seg, ct_datasets, allow_omitted_empty_frames)\n", " segments = segment_metadata(seg)\n", " source_uids = [str(ds.SOPInstanceUID) for ds in ct_datasets]\n", - " # C4KC-KiTS does not set SpatialLocationsPreserved. Geometry was checked independently\n", - " # above; only that missing flag is ignored. Explicit NO has already been rejected.\n", - " masks_dhwk = seg.get_pixels_by_source_instance(\n", - " source_sop_instance_uids=source_uids,\n", - " segment_numbers=[item[\"segment_number\"] for item in segments],\n", - " combine_segments=False,\n", - " ignore_spatial_locations=checks[\"spatial_locations_flag_absent\"],\n", - " assert_missing_frames_are_empty=allow_omitted_empty_frames,\n", - " dtype=np.uint8,\n", - " )\n", - " require(masks_dhwk.shape == (len(ct_datasets), ct_datasets[0].Rows, ct_datasets[0].Columns, len(segments)), \"Decoded SEG shape differs from CT\")\n", + " try:\n", + " mask_volume = seg.get_volume(\n", + " segment_numbers=[item[\"segment_number\"] for item in segments],\n", + " combine_segments=False, dtype=np.uint8,\n", + " allow_missing_positions=allow_omitted_empty_frames, atol=0.01,\n", + " ).match_geometry(ct_volume, tol=1e-4)\n", + " except (ValueError, RuntimeError) as error:\n", + " raise ValidationError(f\"highdicom SEG grid alignment failed: {error}\") from error\n", + " masks_dhwk = mask_volume.array\n", + " require(masks_dhwk.shape == (*ct_volume.array.shape, len(segments)), \"Decoded SEG shape differs from CT\")\n", " require(np.isin(masks_dhwk, (0, 1)).all(), \"Non-binary decoded mask values\")\n", - " image_hu = np.empty((1, len(ct_datasets), ct_datasets[0].Rows, ct_datasets[0].Columns), dtype=np.float32)\n", - " for i, ds in enumerate(ct_datasets):\n", - " require(\"ModalityLUTSequence\" not in ds, \"This example supports linear CT rescaling only\")\n", - " require(str(ds.get(\"RescaleType\", \"HU\")).upper() == \"HU\", \"CT rescale units are not HU\")\n", - " require(\"RescaleSlope\" in ds and \"RescaleIntercept\" in ds, \"CT rescale parameters are missing\")\n", - " slope, intercept = float(ds.RescaleSlope), float(ds.RescaleIntercept)\n", - " require(math.isfinite(slope) and math.isfinite(intercept) and slope != 0, \"Invalid CT rescale parameters\")\n", - " image_hu[0, i] = ds.pixel_array.astype(np.float32) * slope + intercept\n", - " require(np.isfinite(image_hu).all(), \"Non-finite CT intensities after rescaling\")\n", + " image_hu = np.ascontiguousarray(ct_volume.array[np.newaxis], dtype=np.float32)\n", " masks = np.ascontiguousarray(np.moveaxis(masks_dhwk, -1, 0))\n", + " affine = ct_volume.affine\n", + " checks[\"ct_volume_loader\"] = \"highdicom.get_volume_from_series\"\n", + " checks[\"hu_units_from_original_ct_rule\"] = any(\"RescaleType\" not in ds for ds in ct_datasets)\n", + " checks[\"seg_volume_loader\"] = \"highdicom.seg.Segmentation.get_volume + Volume.match_geometry\"\n", + " checks[\"grid_match_tolerance_voxels\"] = 1e-4\n", " first = ct_datasets[0]\n", " metadata = {\n", " \"patient_id\": str(first.PatientID),\n", @@ -509,9 +563,9 @@ "### 2. Load a pinned metadata catalog\n", "\n", "The main catalog has one row per series. `seg_index` supplies the referenced\n", - "CT series. We compare selected cloud object identifiers with the download\n", + "CT series, and `volume_geometry_index` supplies precomputed regular-grid checks. We compare selected cloud object identifiers with the download\n", "client's package index, so a catalog/client mismatch fails before downloading.\n", - "The Hugging Face option downloads two Parquet metadata files and records their\n", + "The Hugging Face option downloads three Parquet metadata files and records their\n", "SHA256 checksums; this can be larger than the small DICOM example itself." ] }, @@ -528,7 +582,7 @@ "if CATALOG_SOURCE == \"huggingface\":\n", " cache_dir = WORK_DIR / \"catalog\" / HF_REVISION\n", " cache_dir.mkdir(parents=True, exist_ok=True)\n", - " for name in (\"idc_index\", \"seg_index\"):\n", + " for name in (\"idc_index\", \"seg_index\", \"volume_geometry_index\"):\n", " path = cache_dir / f\"{name}.parquet\"\n", " if not path.exists():\n", " url = f\"https://huggingface.co/datasets/ImagingDataCommons/idc-index-data/resolve/{HF_REVISION}/{name}.parquet\"\n", @@ -542,12 +596,16 @@ " catalog_paths[name] = path\n", " index_table = pd.read_parquet(catalog_paths[\"idc_index\"])\n", " seg_table = pd.read_parquet(catalog_paths[\"seg_index\"])\n", + " geometry_table = pd.read_parquet(catalog_paths[\"volume_geometry_index\"])\n", "elif CATALOG_SOURCE == \"idc-index\":\n", " client.fetch_index(\"seg_index\")\n", + " client.fetch_index(\"volume_geometry_index\")\n", " require(getattr(client, \"seg_index\", None) is not None, \"SEG index download failed\")\n", - " index_table, seg_table = client.index, client.seg_index\n", + " require(getattr(client, \"volume_geometry_index\", None) is not None, \"Geometry index download failed\")\n", + " index_table, seg_table, geometry_table = client.index, client.seg_index, client.volume_geometry_index\n", " catalog_paths = {\"idc_index\": Path(client.indices_overview[\"index\"][\"file_path\"]),\n", - " \"seg_index\": Path(client.indices_overview[\"seg_index\"][\"file_path\"])}\n", + " \"seg_index\": Path(client.indices_overview[\"seg_index\"][\"file_path\"]),\n", + " \"volume_geometry_index\": Path(client.indices_overview[\"volume_geometry_index\"][\"file_path\"])}\n", "else:\n", " raise ValidationError(\"Unknown CATALOG_SOURCE\")\n", "\n", @@ -568,7 +626,8 @@ "### 3. Select four referenced pairs and inspect the download estimate\n", "\n", "Restrict both modalities to this collection and the documented CC BY 3.0 license.\n", - "Require one SEG instance per series, retain one pair per patient, and show the\n", + "Require one SEG instance per series and a CT marked `regularly_spaced_3d_volume`\n", + "in the geometry index, retain one pair per patient, and show the\n", "estimate before any DICOM download. Raise an error if the configured size cap\n", "or patient count cannot be met. Check the [collection's usage and citation requirements](https://www.cancerimagingarchive.net/collection/c4kc-kits/)." ] @@ -581,22 +640,7 @@ "outputs": [], "source": [ "cohort = index_table.loc[index_table.collection_id.eq(COLLECTION)].copy()\n", - "seg_rows = cohort.loc[cohort.Modality.eq(\"SEG\") & cohort.instanceCount.eq(1)].merge(\n", - " seg_table[[\"SeriesInstanceUID\", \"segmented_SeriesInstanceUID\", \"SegmentationType\"]],\n", - " on=\"SeriesInstanceUID\", validate=\"one_to_one\")\n", - "seg_rows = seg_rows.loc[seg_rows.SegmentationType.eq(\"BINARY\")]\n", - "ct_rows = cohort.loc[cohort.Modality.eq(\"CT\")]\n", - "pairs = seg_rows.merge(ct_rows, left_on=\"segmented_SeriesInstanceUID\", right_on=\"SeriesInstanceUID\",\n", - " suffixes=(\"_seg\", \"_ct\"), validate=\"many_to_one\")\n", - "pairs = pairs.loc[pairs.PatientID_seg.eq(pairs.PatientID_ct)\n", - " & pairs.StudyInstanceUID_seg.eq(pairs.StudyInstanceUID_ct)\n", - " & pairs.license_short_name_seg.eq(\"CC BY 3.0\")\n", - " & pairs.license_short_name_ct.eq(\"CC BY 3.0\")].copy()\n", - "pairs[\"estimated_mb\"] = pairs.series_size_MB_seg + pairs.series_size_MB_ct\n", - "pairs = pairs.sort_values([\"estimated_mb\", \"PatientID_ct\", \"SeriesInstanceUID_seg\"], kind=\"stable\")\n", - "pairs = pairs.drop_duplicates(\"PatientID_ct\").head(N_PATIENTS)\n", - "require(N_PATIENTS >= 2 and len(pairs) == N_PATIENTS, \"Not enough eligible patient pairs\")\n", - "require(np.isfinite(pairs.estimated_mb).all() and (pairs.estimated_mb > 0).all(), \"Invalid size estimates\")\n", + "pairs = select_catalog_pairs(index_table, seg_table, geometry_table, COLLECTION, N_PATIENTS, MAX_DICOM_MB)\n", "estimated_mb = float(pairs.estimated_mb.sum())\n", "require(estimated_mb <= MAX_DICOM_MB, \"Cohort exceeds MAX_DICOM_MB; reduce N_PATIENTS\")\n", "selected_uids = pairs.SeriesInstanceUID_ct.tolist() + pairs.SeriesInstanceUID_seg.tolist()\n", @@ -663,6 +707,7 @@ " tensor_file = f\"{sample_id}.pt\"\n", " torch.save({\"image_hu\": torch.from_numpy(image_hu), \"masks\": torch.from_numpy(masks)}, OUTPUT_DIR / tensor_file)\n", " metadata[\"collection_id\"] = COLLECTION\n", + " metadata[\"catalog_regularly_spaced_3d_volume\"] = bool(row[\"regularly_spaced_3d_volume\"])\n", " metadata[\"sources\"] = [{\"path\": str(path.relative_to(DICOM_DIR)), \"sha256\": file_sha256(path),\n", " \"size_bytes\": path.stat().st_size} for path in ct_paths + seg_paths]\n", " metadata[\"license_ct\"] = row[\"license_short_name_ct\"]\n", @@ -722,7 +767,8 @@ " \"catalog_provenance\": catalog_provenance,\n", " \"packages\": {name: version(name) for name in (\"idc-index\", \"idc-index-data\", \"highdicom\", \"pydicom\", \"numpy\", \"pandas\", \"torch\", \"matplotlib\")},\n", " \"split_seed\": SPLIT_SEED, \"test_fraction\": TEST_FRACTION,\n", - " \"sampling\": \"smallest referenced CT/SEG pair per patient; tutorial convenience sample\",\n", + " \"sampling\": \"smallest referenced regular-grid CT/SEG pair per patient; tutorial convenience sample\",\n", + " \"geometry_preflight\": \"volume_geometry_index.regularly_spaced_3d_volume\",\n", " \"allow_omitted_empty_frames\": ALLOW_OMITTED_EMPTY_FRAMES,\n", " \"channel_semantic_keys\": canonical_channel_keys,\n", " \"estimated_dicom_mb\": estimated_mb,\n", @@ -815,6 +861,7 @@ " require(tensors[\"masks\"].dtype == torch.uint8 and torch.all((tensors[\"masks\"] == 0) | (tensors[\"masks\"] == 1)).item(), \"Invalid binary masks\")\n", " require([segment[\"semantic_key\"] for segment in metadata[\"segments\"]] == saved[\"channel_semantic_keys\"], \"Saved channels differ between cases\")\n", " require(len(metadata[\"source_sop_instance_uids_in_slice_order\"]) == record[\"ct_slices\"], \"Saved source count mismatch\")\n", + " require(metadata[\"catalog_regularly_spaced_3d_volume\"] is True, \"Missing catalog geometry preflight\")\n", " del tensors\n", "print(\"PASS: validated real CT/SEG pairs, consistent tensor channels, and disjoint patient groups\")" ] diff --git a/test/src/test_ct_seg_to_tensors.py b/test/src/test_ct_seg_to_tensors.py index 293dbe6..e4ff02c 100644 --- a/test/src/test_ct_seg_to_tensors.py +++ b/test/src/test_ct_seg_to_tensors.py @@ -9,6 +9,7 @@ import highdicom as hd import nbformat import numpy as np +import pandas as pd import pydicom import pytest from pydicom.dataset import Dataset, FileMetaDataset @@ -23,6 +24,7 @@ convert = namespace["validate_and_convert"] split_patients = namespace["patient_disjoint_split"] ValidationError = namespace["ValidationError"] +select_pairs = namespace["select_catalog_pairs"] @pytest.fixture @@ -98,11 +100,14 @@ def test_ct_order_hu_affine_and_overlapping_mask_channels(pair): image, masks, metadata = convert([images[2], images[0], images[1]], seg, 3) assert image.shape == (1, 3, 4, 4) assert image.dtype == np.float32 - assert image[0, :, 0, 0].tolist() == [-998, -988, -978] + # highdicom's default volume order runs from superior to inferior. + assert image[0, :, 0, 0].tolist() == [-978, -988, -998] np.testing.assert_allclose(metadata["affine_dhw_to_lps_mm"], - [[0, 0, 3, 0], [0, 2, 0, 0], [5, 0, 0, 0], [0, 0, 0, 1]]) + [[0, 0, 3, 0], [0, 2, 0, 0], [-5, 0, 0, 10], [0, 0, 0, 1]]) + assert metadata["source_sop_instance_uids_in_slice_order"] == [ + str(ds.SOPInstanceUID) for ds in reversed(images)] assert masks.dtype == np.uint8 - assert masks[:, 2, 2, 1].tolist() == [1, 1] + assert masks[:, 0, 2, 1].tolist() == [1, 1] assert sorted(np.count_nonzero(mask) for mask in masks) == [1, 2] @@ -138,8 +143,8 @@ def test_incompatible_pair_is_rejected(pair, change, message): @pytest.mark.parametrize("change, message", [ ("missing", "instance count"), ("duplicate_uid", "Duplicate CT SOP"), - ("irregular", "Irregular CT slice spacing"), - ("gantry_tilt", "Gantry tilt"), + ("irregular", "highdicom CT volume validation failed"), + ("gantry_tilt", "highdicom CT volume validation failed"), ]) def test_incomplete_or_unsupported_ct_is_rejected(pair, change, message): images, seg = pair @@ -163,7 +168,115 @@ def test_missing_spatial_flag_is_allowed_only_after_geometry_validation(pair): del source.SpatialLocationsPreserved _, masks, metadata = convert(images, seg, 3) assert metadata["checks"]["spatial_locations_flag_absent"] is True - assert masks[:, 2, 2, 1].tolist() == [1, 1] + assert masks[:, 0, 2, 1].tolist() == [1, 1] + + +def test_hu_values_do_not_receive_display_or_presentation_transforms(pair): + images, seg = pair + for image in images: + image.PhotometricInterpretation = "MONOCHROME1" + image.PresentationLUTShape = "INVERSE" + image.WindowCenter = -500 + image.WindowWidth = 100 + image, _, _ = convert(images, seg, 3) + assert image[0, :, 0, 0].tolist() == [-978, -988, -998] + + +def test_unspecified_rescale_units_are_not_silently_labelled_hu(pair): + images, seg = pair + del images[0].RescaleType + with pytest.raises(ValidationError, match="cannot be verified as HU"): + convert(images, seg, 3) + + +def test_original_ct_can_omit_rescale_type_under_dicom_hu_rule(pair): + images, seg = pair + for image in images: + del image.RescaleType + image.ImageType = ["ORIGINAL", "PRIMARY", "AXIAL"] + image, _, metadata = convert(images, seg, 3) + assert image[0, :, 0, 0].tolist() == [-978, -988, -998] + assert metadata["checks"]["hu_units_from_original_ct_rule"] is True + + +@pytest.mark.parametrize("image_type,multi_energy", [ + (["DERIVED", "PRIMARY", "AXIAL"], "NO"), + (["ORIGINAL", "PRIMARY", "LOCALIZER"], "NO"), + (["ORIGINAL", "PRIMARY", "AXIAL"], "YES"), +]) +def test_missing_rescale_units_are_not_inferred_outside_original_ct_rule(pair, image_type, multi_energy): + images, seg = pair + for image in images: + del image.RescaleType + image.ImageType = image_type + image.MultienergyCTAcquisition = multi_energy + with pytest.raises(ValidationError, match="cannot be verified as HU"): + convert(images, seg, 3) + + +def test_omitted_empty_seg_frames_need_explicit_policy_and_align_to_full_ct(pair): + images, seg = pair + pixels = seg.get_pixels_by_source_instance( + [str(image.SOPInstanceUID) for image in images], combine_segments=False) + sparse_seg = hd.seg.Segmentation( + source_images=images, pixel_array=pixels, + segmentation_type=hd.seg.SegmentationTypeValues.BINARY, + segment_descriptions=list(seg.SegmentSequence), omit_empty_frames=True, + series_instance_uid=generate_uid(), series_number=2, + sop_instance_uid=generate_uid(), instance_number=1, + manufacturer="Tutorial", manufacturer_model_name="Synthetic", + software_versions="1", device_serial_number="1", + ) + with pytest.raises(ValidationError, match="explicit empty-frame policy"): + convert(images, sparse_seg, 3) + image, masks, metadata = convert(images, sparse_seg, 3, allow_omitted_empty_frames=True) + assert image.shape == (1, 3, 4, 4) + assert masks.shape == (2, 3, 4, 4) + assert not masks[:, 1].any() + assert masks[:, 0, 2, 1].tolist() == [1, 1] + assert metadata["checks"]["omitted_source_segment_frames"] > 0 + + +@pytest.fixture +def catalog(): + rows, segmentations, geometries = [], [], [] + for patient, size in [("A", 1.0), ("B", 2.0), ("C", 3.0)]: + for modality in ["CT", "SEG"]: + rows.append({ + "collection_id": "c4kc_kits", "PatientID": patient, + "StudyInstanceUID": f"study-{patient}", + "SeriesInstanceUID": f"{modality}-{patient}", "Modality": modality, + "instanceCount": 3 if modality == "CT" else 1, + "series_size_MB": size, "license_short_name": "CC BY 3.0", + }) + segmentations.append({"SeriesInstanceUID": f"SEG-{patient}", + "segmented_SeriesInstanceUID": f"CT-{patient}", + "SegmentationType": "BINARY"}) + geometries.append({"SeriesInstanceUID": f"CT-{patient}", + "regularly_spaced_3d_volume": True}) + return pd.DataFrame(rows), pd.DataFrame(segmentations), pd.DataFrame(geometries) + + +@pytest.mark.parametrize("geometry_value", [False, None, pd.NA]) +def test_catalog_geometry_filter_excludes_failed_or_unknown_ct(catalog, geometry_value): + index, segmentations, geometry = catalog + geometry["regularly_spaced_3d_volume"] = geometry.regularly_spaced_3d_volume.astype(object) + geometry.loc[0, "regularly_spaced_3d_volume"] = geometry_value + pairs = select_pairs(index, segmentations, geometry, "c4kc_kits", 2, 100) + assert pairs.PatientID_ct.tolist() == ["B", "C"] + + +def test_catalog_missing_geometry_and_wrong_reference_do_not_fall_back_to_patient(catalog): + index, segmentations, geometry = catalog + geometry = geometry.iloc[1:] + segmentations.loc[1, "segmented_SeriesInstanceUID"] = "CT-A" + with pytest.raises(ValidationError, match="Not enough eligible"): + select_pairs(index, segmentations, geometry, "c4kc_kits", 2, 100) + + +def test_catalog_size_cap_fails_before_download(catalog): + with pytest.raises(ValidationError, match="exceeds MAX_DICOM_MB"): + select_pairs(*catalog, "c4kc_kits", 2, 1) def test_patient_split_is_grouped_and_independent_of_record_order(): From 3c9d62b93f64036fcb5b9df91349aecd00a34501 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Wed, 21 Oct 2026 16:48:14 -0400 Subject: [PATCH 3/5] Simplify tutorial notebook --- .../advanced_topics/ct_seg_to_tensors.ipynb | 966 +++++------------- 1 file changed, 240 insertions(+), 726 deletions(-) diff --git a/notebooks/advanced_topics/ct_seg_to_tensors.ipynb b/notebooks/advanced_topics/ct_seg_to_tensors.ipynb index bfd3af3..fe4880a 100644 --- a/notebooks/advanced_topics/ct_seg_to_tensors.ipynb +++ b/notebooks/advanced_topics/ct_seg_to_tensors.ipynb @@ -2,924 +2,438 @@ "cells": [ { "cell_type": "markdown", - "id": "87781114", + "id": "951549ba", "metadata": {}, "source": [ "\"Open\n", "\n", - "# Validated CT/DICOM SEG tensors and patient-disjoint splits\n", + "# Loading CT images and DICOM segmentations into PyTorch tensors\n", "\n", "Updated: October 2026\n", "\n", - "## Goal\n", - "\n", - "Turn a small public **C4KC-KiTS** cohort into traceable PyTorch inputs: an HU image\n", - "`[1, D, H, W]`, separate binary mask channels `[K, D, H, W]`, an affine in DICOM\n", - "patient coordinates, and a patient-disjoint train/test manifest.\n", - "\n", - "The [IDC Hugging Face dataset](https://huggingface.co/datasets/ImagingDataCommons/idc-index-data)\n", - "is a **metadata catalog**, not pixel data. Its `train` split is a Hub convention,\n", - "not a medical ML training partition. This notebook queries a pinned catalog and\n", - "uses `idc-index` to download only selected DICOM series. No cloud account or GPU\n", - "is required. Start with [Getting Started Part 2](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/getting_started/part2_searching_basics.ipynb)\n", - "if you are new to IDC.\n", - "\n", - "### Scope and assumptions\n", - "\n", - "- Four patients, selected deterministically by the smallest combined CT/SEG download.\n", - " This deliberately biased convenience sample is for learning, not model evaluation.\n", - "- Legacy single-frame CT on a regular grid, linear rescaling to HU, and one BINARY\n", - " SEG instance per series. Enhanced CT, non-regular spacing, gantry tilt, and\n", - " resampled masks are rejected. Nothing is silently resampled.\n", - "- `seg_index.segmented_SeriesInstanceUID` selects the CT, and\n", - " `volume_geometry_index.regularly_spaced_3d_volume` filters it before download.\n", - " `highdicom` loads HU volumes and aligns separate SEG channels without resampling.\n", - " Source SOP Instance UIDs and SEG-to-source frame correspondence remain local\n", - " integrity checks; a patient/study match is insufficient.\n", - "- Separate channels retain overlaps. Coded category/type, modifiers and anatomic\n", - " region define a stable channel order across cases; free-text labels and local\n", - " segment numbers alone do not define classes.\n", - "- Missing CT files are errors. Omitted mask frames are also errors unless the\n", - " caller explicitly enables an empty-frame policy for a verified dataset.\n", - "- The instructional split is within the public KiTS19 training collection; it is\n", - " not the official challenge test set. Cross-collection identity reconciliation\n", - " and statistically representative sampling are outside this example." + "Many collections in [NCI Imaging Data Commons (IDC)](https://imaging.datacommons.cancer.gov/) contain CT\n", + "images together with segmentations stored as [DICOM Segmentation (SEG)](https://highdicom.readthedocs.io/en/latest/seg.html)\n", + "objects. These are a great resource for training deep learning segmentation models, but there are a few\n", + "steps between \"a folder of DICOM files\" and \"a pair of tensors ready for a model\".\n", + "\n", + "In this tutorial we will:\n", + "\n", + "1. Use `idc-index` to find CT series that have a corresponding SEG, and download a few of them.\n", + "2. Use `highdicom` to load the CT series as a 3D volume in Hounsfield units (HU).\n", + "3. Use `highdicom` to load the SEG as a 3D mask array that lines up voxel-for-voxel with the CT.\n", + "4. Convert both to PyTorch tensors and wrap them in a simple `torch.utils.data.Dataset`.\n", + "\n", + "We use the [C4KC-KiTS](https://www.cancerimagingarchive.net/collection/c4kc-kits/) collection, which contains\n", + "abdominal CTs with manual segmentations of the kidneys and kidney tumors.\n", + "\n", + "This tutorial assumes basic familiarity with DICOM, `highdicom` and PyTorch. If you are new to IDC, start with\n", + "the [Getting Started](https://github.com/ImagingDataCommons/IDC-Tutorials/tree/master/notebooks/getting_started)\n", + "notebooks. No cloud account or GPU is needed." ] }, { "cell_type": "markdown", - "id": "af3feb50", + "id": "239b5820", "metadata": {}, "source": [ "## Setup\n", "\n", - "Install the versions used for this example in a fresh Python 3.11/3.12 environment.\n", - "The PyTorch wheel is CPU-only. Local Jupyter users also need `ipykernel` and\n", - "`nbclient` to execute the notebook programmatically." + "Install the packages we need. `idc-index` is used to search and download IDC data, and `highdicom`\n", + "(which builds on `pydicom`) is used to read the images and segmentations. If this upgrades a package that was\n", + "already installed, restart the kernel before continuing." ] }, { "cell_type": "code", "execution_count": null, - "id": "cc43facd", + "id": "091bbfb2", "metadata": {}, "outputs": [], "source": [ - "import subprocess\n", - "import sys\n", - "\n", - "# check_call stops execution if either installation fails.\n", - "subprocess.check_call([\n", - " sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"idc-index==0.13.0\",\n", - " \"pydicom==3.0.2\", \"highdicom==0.28.2\", \"numpy==2.3.5\",\n", - " \"pandas==2.2.3\", \"matplotlib==3.10.8\",\n", - "])\n", - "subprocess.check_call([\n", - " sys.executable, \"-m\", \"pip\", \"install\", \"-q\", \"torch==2.14.1\",\n", - " \"--index-url\", \"https://download.pytorch.org/whl/cpu\",\n", - "])" - ] - }, - { - "cell_type": "markdown", - "id": "b66ef835", - "metadata": {}, - "source": [ - "### Parameters\n", - "\n", - "The default uses Hugging Face revision `93983d810bea0454ae1a400ab62c48dd01f8d0ca`\n", - "of the IDC v25 catalog. Set `CATALOG_SOURCE = \"idc-index\"` to use the corresponding\n", - "package indices instead. Both routes are checked for matching release metadata.\n", - "\n", - "Metadata downloads are separate from DICOM download size estimates. All paths\n", - "below are relative to the working directory. Rerunning reuses existing DICOM\n", - "files and replaces generated outputs for the same sample." + "%pip install -q \"idc-index>=0.13\" \"highdicom>=0.28\" torch matplotlib" ] }, { "cell_type": "code", "execution_count": null, - "id": "af68526c", - "metadata": {}, - "outputs": [], - "source": [ - "from importlib.metadata import version\n", - "from pathlib import Path\n", - "from time import perf_counter\n", - "import json\n", - "import hashlib\n", - "\n", - "import matplotlib.pyplot as plt\n", - "import pandas as pd\n", - "import pyarrow.parquet as pq\n", - "import requests\n", - "import torch\n", - "from idc_index import IDCClient\n", - "from IPython.display import display\n", - "\n", - "CATALOG_SOURCE = \"huggingface\" # or \"idc-index\"\n", - "HF_REVISION = \"93983d810bea0454ae1a400ab62c48dd01f8d0ca\"\n", - "COLLECTION = \"c4kc_kits\"\n", - "N_PATIENTS = 4\n", - "SPLIT_SEED = 2026\n", - "TEST_FRACTION = 0.25\n", - "MAX_DICOM_MB = 200.0\n", - "ALLOW_OMITTED_EMPTY_FRAMES = False\n", - "WORK_DIR = Path(\"idc_ct_seg_data\")\n", - "DICOM_DIR = WORK_DIR / \"dicom\"\n", - "OUTPUT_DIR = WORK_DIR / \"outputs\"\n", - "OUTPUT_DIR.mkdir(parents=True, exist_ok=True)\n", - "run_started = perf_counter()" - ] - }, - { - "cell_type": "markdown", - "id": "ff60e33c", - "metadata": {}, - "source": [ - "## Steps\n", - "\n", - "### 1. Use the catalog and highdicom at their intended boundaries\n", - "\n", - "The [IDC index reference](https://idc-index.readthedocs.io/en/latest/indices_reference.html)\n", - "and [tensor-loading example](https://huggingface.co/datasets/ImagingDataCommons/idc-index-data#loading-images-as-tensors)\n", - "provide the building blocks used below:\n", - "\n", - "- `seg_index` supplies CT/SEG correspondence, without reconstructing it manually.\n", - "- `volume_geometry_index` excludes CT series that are not regular 3D grids\n", - " before any DICOM download.\n", - "- `hd.get_volume_from_series()` orders CT slices, checks the actual downloaded\n", - " geometry, applies the modality transform to HU, and retains the LPS affine.\n", - "- `seg.get_volume(combine_segments=False).match_geometry(ct_volume)` aligns\n", - " separate binary channels by padding/flipping/permuting; it does not interpolate.\n", - "\n", - "The remaining local guards check download completeness, declared HU units,\n", - "source UIDs, and each SEG frame's correspondence to its referenced CT. The CT\n", - "geometry index does not establish those pair-specific properties. C4KC-KiTS\n", - "omits `SpatialLocationsPreserved`; absence is accepted only after these checks,\n", - "while explicit `NO` is rejected. SEG position tolerance is 0.01 mm and orientation\n", - "and pixel-spacing tolerance is `1e-4`. Highdicom CT spacing tolerance is 0.01 mm;\n", - "the final grid match uses `1e-4` voxels to tolerate rounding only.\n", - "\n", - "The exported affine maps the library's array indices `[D, H, W]` to **DICOM LPS\n", - "millimeters**. Slice direction comes from highdicom and may differ from file or\n", - "InstanceNumber order. No display window or real-world-value transform is applied\n", - "to the saved HU image. HU units must be explicit or guaranteed by the\n", - "[original CT rule in DICOM PS3.3 C.8.2.1](https://dicom.nema.org/medical/Dicom/current/output/chtml/part03/sect_C.8.2.html);\n", - "missing units on derived/localizer/multi-energy images are not assumed to be HU." - ] - }, - { - "cell_type": "markdown", - "id": "5c7be955", + "id": "51d509ad", "metadata": {}, - "source": [ - "### Validation helpers" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "da62205f", - "metadata": { - "tags": [ - "definitions" - ] - }, "outputs": [], "source": [ - "from __future__ import annotations\n", - "\n", - "import hashlib\n", - "import json\n", - "import math\n", - "import random\n", "from pathlib import Path\n", "\n", - "import pandas as pd\n", "import highdicom as hd\n", + "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pydicom\n", - "from pydicom.uid import CTImageStorage\n", - "\n", - "\n", - "class ValidationError(ValueError):\n", - " \"\"\"A CT/SEG pair cannot be interpreted safely within this tutorial's scope.\"\"\"\n", - "\n", - "\n", - "def require(condition: bool, message: str) -> None:\n", - " if not condition:\n", - " raise ValidationError(message)\n", - "\n", - "\n", - "def numeric_vector(value, length: int, name: str) -> np.ndarray:\n", - " result = np.asarray(value, dtype=np.float64)\n", - " require(result.shape == (length,) and np.isfinite(result).all(), f\"Invalid {name}\")\n", - " return result" - ] - }, - { - "cell_type": "markdown", - "id": "b6631e97", - "metadata": {}, - "source": [ - "### Load the CT volume with highdicom" + "import torch\n", + "from idc_index import IDCClient" ] }, { "cell_type": "code", "execution_count": null, - "id": "c9dd0578", - "metadata": { - "tags": [ - "definitions" - ] - }, - "outputs": [], - "source": [ - "def load_ct_volume(datasets, expected_count: int):\n", - " \"\"\"Check identity/completeness and delegate volume geometry and HU to highdicom.\"\"\"\n", - " require(len(datasets) == expected_count, \"CT instance count differs from the catalog\")\n", - " require(len(datasets) >= 2, \"A volume needs at least two CT slices\")\n", - " first = datasets[0]\n", - " uids = [str(ds.SOPInstanceUID) for ds in datasets]\n", - " require(len(set(uids)) == len(uids), \"Duplicate CT SOP Instance UID\")\n", - " for ds in datasets:\n", - " require(str(ds.SOPClassUID) == str(CTImageStorage), \"Only legacy single-frame CT is supported\")\n", - " require(int(ds.get(\"NumberOfFrames\", 1)) == 1, \"Enhanced/multiframe CT is unsupported\")\n", - " require(ds.Modality == \"CT\", \"Non-CT instance in the CT series\")\n", - " for key in (\"PatientID\", \"StudyInstanceUID\", \"SeriesInstanceUID\", \"FrameOfReferenceUID\"):\n", - " require(str(ds.get(key, \"\")) == str(first.get(key, \"\")) != \"\", f\"Inconsistent CT {key}\")\n", - " require(\"ModalityLUTSequence\" not in ds, \"This example supports linear CT rescaling only\")\n", - " # DICOM PS3.3 C.8.2.1: original non-localizer, non-multi-energy CT is HU,\n", - " # and RescaleType is optional in that case. Do not infer units for derived CT.\n", - " image_type = list(ds.get(\"ImageType\", []))\n", - " implicit_hu = (len(image_type) >= 3 and image_type[0] == \"ORIGINAL\"\n", - " and image_type[2] != \"LOCALIZER\"\n", - " and str(ds.get(\"MultienergyCTAcquisition\", \"NO\")) == \"NO\")\n", - " units = str(ds.get(\"RescaleType\", \"HU\" if implicit_hu else \"\")).upper()\n", - " require(units == \"HU\", \"CT rescale units cannot be verified as HU\")\n", - " require(\"RescaleSlope\" in ds and \"RescaleIntercept\" in ds, \"CT rescale parameters are missing\")\n", - " slope, intercept = float(ds.RescaleSlope), float(ds.RescaleIntercept)\n", - " require(math.isfinite(slope) and math.isfinite(intercept) and slope != 0, \"Invalid CT rescale parameters\")\n", - " try:\n", - " volume = hd.get_volume_from_series(\n", - " datasets, dtype=np.float32, apply_real_world_transform=False,\n", - " apply_modality_transform=True, apply_voi_transform=False,\n", - " apply_presentation_lut=False, atol=0.01,\n", - " )\n", - " except (ValueError, RuntimeError) as error:\n", - " raise ValidationError(f\"highdicom CT volume validation failed: {error}\") from error\n", - " require(np.isfinite(volume.array).all(), \"Non-finite CT intensities after rescaling\")\n", - " # Use the library's affine to record provenance in the exported slice order.\n", - " positions = np.asarray([ds.ImagePositionPatient for ds in datasets], dtype=np.float64)\n", - " slice_indices = volume.map_reference_to_indices(positions)[:, 0]\n", - " ordered = [datasets[int(i)] for i in np.argsort(slice_indices, kind=\"stable\")]\n", - " return volume, ordered\n", - "\n", - "\n", - "def select_catalog_pairs(index_table, seg_table, geometry_table, collection, n_patients, max_dicom_mb):\n", - " \"\"\"Select referenced regular CT/SEG pairs before downloading any DICOM.\"\"\"\n", - " require(isinstance(n_patients, int) and n_patients >= 2, \"At least two patient pairs are required\")\n", - " require(math.isfinite(max_dicom_mb) and max_dicom_mb > 0, \"Invalid DICOM size cap\")\n", - " cohort = index_table.loc[index_table.collection_id.eq(collection)].copy()\n", - " seg_rows = cohort.loc[cohort.Modality.eq(\"SEG\") & cohort.instanceCount.eq(1)].merge(\n", - " seg_table[[\"SeriesInstanceUID\", \"segmented_SeriesInstanceUID\", \"SegmentationType\"]],\n", - " on=\"SeriesInstanceUID\", validate=\"one_to_one\")\n", - " seg_rows = seg_rows.loc[seg_rows.SegmentationType.eq(\"BINARY\")]\n", - " eligible_geometry = geometry_table.loc[\n", - " geometry_table.regularly_spaced_3d_volume.eq(True).fillna(False),\n", - " [\"SeriesInstanceUID\", \"regularly_spaced_3d_volume\"],\n", - " ]\n", - " ct_rows = cohort.loc[cohort.Modality.eq(\"CT\")].merge(\n", - " eligible_geometry, on=\"SeriesInstanceUID\", validate=\"one_to_one\")\n", - " pairs = seg_rows.merge(ct_rows, left_on=\"segmented_SeriesInstanceUID\", right_on=\"SeriesInstanceUID\",\n", - " suffixes=(\"_seg\", \"_ct\"), validate=\"many_to_one\")\n", - " pairs = pairs.loc[pairs.PatientID_seg.eq(pairs.PatientID_ct)\n", - " & pairs.StudyInstanceUID_seg.eq(pairs.StudyInstanceUID_ct)\n", - " & pairs.license_short_name_seg.eq(\"CC BY 3.0\")\n", - " & pairs.license_short_name_ct.eq(\"CC BY 3.0\")].copy()\n", - " pairs[\"estimated_mb\"] = pairs.series_size_MB_seg + pairs.series_size_MB_ct\n", - " pairs = pairs.sort_values([\"estimated_mb\", \"PatientID_ct\", \"SeriesInstanceUID_seg\"], kind=\"stable\")\n", - " pairs = pairs.drop_duplicates(\"PatientID_ct\").head(n_patients)\n", - " require(len(pairs) == n_patients, \"Not enough eligible patient pairs with regular CT geometry\")\n", - " require(np.isfinite(pairs.estimated_mb).all() and (pairs.estimated_mb > 0).all(), \"Invalid size estimates\")\n", - " require(float(pairs.estimated_mb.sum()) <= max_dicom_mb, \"Cohort exceeds MAX_DICOM_MB; reduce N_PATIENTS\")\n", - " return pairs" - ] - }, - { - "cell_type": "markdown", - "id": "715888b1", + "id": "1243f56f-3f8a-45e2-af6b-94bdf1268f9c", "metadata": {}, - "source": [ - "### Check every SEG frame against its referenced CT" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "01ae8c29", - "metadata": { - "tags": [ - "definitions" - ] - }, "outputs": [], "source": [ - "def functional_group(frame, shared, sequence_name):\n", - " sequence = frame.get(sequence_name, shared.get(sequence_name, None))\n", - " require(sequence is not None and len(sequence) == 1, f\"Missing/ambiguous SEG {sequence_name}\")\n", - " return sequence[0]\n", - "\n", - "\n", - "def validate_seg_geometry(seg, ct_datasets, allow_omitted_empty_frames: bool = False):\n", - " \"\"\"Verify source references and every SEG frame on the CT grid before decoding.\"\"\"\n", - " require(seg.Modality == \"SEG\" and str(seg.SegmentationType) == \"BINARY\", \"Only BINARY DICOM SEG is supported\")\n", - " ct = ct_datasets[0]\n", - " for key in (\"PatientID\", \"StudyInstanceUID\", \"FrameOfReferenceUID\"):\n", - " require(str(seg.get(key, \"\")) == str(ct.get(key, \"\")) != \"\", f\"CT/SEG {key} mismatch\")\n", - " require((seg.Rows, seg.Columns) == (ct.Rows, ct.Columns), \"CT/SEG dimensions differ\")\n", - " series_refs = seg.get(\"ReferencedSeriesSequence\", [])\n", - " require(len(series_refs) == 1, \"SEG must reference exactly one source series\")\n", - " require(str(series_refs[0].SeriesInstanceUID) == str(ct.SeriesInstanceUID), \"SEG references a different CT series\")\n", - " referenced_instances = series_refs[0].get(\"ReferencedInstanceSequence\", [])\n", - " referenced_uids = [str(item.ReferencedSOPInstanceUID) for item in referenced_instances]\n", - " require(len(referenced_uids) == len(set(referenced_uids)), \"Duplicate SEG source reference\")\n", - " ct_by_uid = {str(ds.SOPInstanceUID): ds for ds in ct_datasets}\n", - " require(set(referenced_uids) == set(ct_by_uid), \"SEG source references and complete CT series do not match\")\n", - " require(all(str(item.ReferencedSOPClassUID) == str(CTImageStorage) for item in referenced_instances), \"SEG references unsupported source classes\")\n", - " shared_sequence = seg.get(\"SharedFunctionalGroupsSequence\", [])\n", - " require(len(shared_sequence) == 1, \"SEG needs one Shared Functional Groups item\")\n", - " shared = shared_sequence[0]\n", - " frames = seg.get(\"PerFrameFunctionalGroupsSequence\", [])\n", - " require(len(frames) == int(seg.NumberOfFrames), \"SEG functional-group frame count mismatch\")\n", - " segment_numbers = [int(item.SegmentNumber) for item in seg.SegmentSequence]\n", - " require(segment_numbers and len(set(segment_numbers)) == len(segment_numbers), \"Invalid SEG segment numbering\")\n", - " seen = set()\n", - " preserved_flags = []\n", - " for frame in frames:\n", - " derivations = frame.get(\"DerivationImageSequence\", shared.get(\"DerivationImageSequence\", []))\n", - " sources = [source for derivation in derivations for source in derivation.get(\"SourceImageSequence\", [])]\n", - " require(len(sources) == 1, \"Each SEG frame must reference one single-frame CT instance\")\n", - " source = sources[0]\n", - " uid = str(source.ReferencedSOPInstanceUID)\n", - " require(uid in ct_by_uid, \"SEG frame references a CT instance that is missing\")\n", - " require(str(source.ReferencedSOPClassUID) == str(CTImageStorage), \"SEG frame source is not single-frame CT\")\n", - " require(\"ReferencedFrameNumber\" not in source, \"Multiframe source references are unsupported\")\n", - " preserved = str(source.get(\"SpatialLocationsPreserved\", \"\"))\n", - " require(preserved in (\"\", \"YES\"), \"SEG explicitly does not preserve spatial locations\")\n", - " preserved_flags.append(preserved)\n", - " source_ct = ct_by_uid[uid]\n", - " orientation = functional_group(frame, shared, \"PlaneOrientationSequence\").ImageOrientationPatient\n", - " position = functional_group(frame, shared, \"PlanePositionSequence\").ImagePositionPatient\n", - " pixel_spacing = functional_group(frame, shared, \"PixelMeasuresSequence\").PixelSpacing\n", - " require(np.allclose(numeric_vector(orientation, 6, \"SEG orientation\"), source_ct.ImageOrientationPatient, rtol=0, atol=1e-4), \"SEG orientation differs from its source CT\")\n", - " require(np.allclose(numeric_vector(position, 3, \"SEG position\"), source_ct.ImagePositionPatient, rtol=0, atol=0.01), \"SEG position differs from its source CT\")\n", - " require(np.allclose(numeric_vector(pixel_spacing, 2, \"SEG pixel spacing\"), source_ct.PixelSpacing, rtol=0, atol=1e-4), \"SEG pixel spacing differs from its source CT\")\n", - " number = int(functional_group(frame, shared, \"SegmentIdentificationSequence\").ReferencedSegmentNumber)\n", - " require(number in segment_numbers, \"SEG frame names an undefined segment\")\n", - " require((uid, number) not in seen, \"Duplicate SEG frame for one source instance and segment\")\n", - " seen.add((uid, number))\n", - " missing = {(uid, number) for uid in ct_by_uid for number in segment_numbers} - seen\n", - " require(not missing or allow_omitted_empty_frames, \"SEG omits source/segment frames; explicit empty-frame policy is required\")\n", - " return {\n", - " \"segment_numbers\": segment_numbers,\n", - " \"omitted_source_segment_frames\": len(missing),\n", - " \"spatial_locations_flag_absent\": any(flag == \"\" for flag in preserved_flags),\n", - " }" + "import idc_index\n", + "print(idc_index.__version__)\n" ] }, { "cell_type": "markdown", - "id": "615e8f2a", + "id": "21238b92", "metadata": {}, "source": [ - "### Preserve coded segment semantics" + "## 1. Find CT series that have a segmentation\n", + "\n", + "`idc-index` ships a table, `index`, with one row per DICOM series in IDC. It also provides additional tables\n", + "that can be fetched on demand. One of these, `seg_index`, has one row per SEG series, and its\n", + "`segmented_SeriesInstanceUID` column gives the `SeriesInstanceUID` of the image series that was segmented.\n", + "This means we can find CT/SEG pairs with a single SQL query, without opening any DICOM files.\n", + "\n", + "To keep the download small, we pick the three smallest pairs from C4KC-KiTS." ] }, { "cell_type": "code", "execution_count": null, - "id": "2751c843", - "metadata": { - "tags": [ - "definitions" - ] - }, + "id": "203205e1", + "metadata": {}, "outputs": [], "source": [ - "def code_metadata(code):\n", - " return {\n", - " \"value\": str(code.CodeValue),\n", - " \"scheme\": str(code.CodingSchemeDesignator),\n", - " \"meaning\": str(code.CodeMeaning),\n", - " }\n", - "\n", - "\n", - "def segment_metadata(seg):\n", - " result = []\n", - " for item in seg.SegmentSequence:\n", - " category = code_metadata(item.SegmentedPropertyCategoryCodeSequence[0])\n", - " property_type = code_metadata(item.SegmentedPropertyTypeCodeSequence[0])\n", - " modifiers = [code_metadata(c) for c in item.SegmentedPropertyTypeCodeSequence[0].get(\"SegmentedPropertyTypeModifierCodeSequence\", [])]\n", - " regions = []\n", - " for region in item.get(\"AnatomicRegionSequence\", []):\n", - " regions.append({\"code\": code_metadata(region), \"modifiers\": [code_metadata(c) for c in region.get(\"AnatomicRegionModifierSequence\", [])]})\n", - " # Include modifiers/anatomic region in the semantic key, not just the free-text label.\n", - " def identity(code):\n", - " return (code[\"scheme\"], code[\"value\"])\n", - "\n", - " semantic_key = json.dumps({\n", - " \"category\": identity(category), \"type\": identity(property_type),\n", - " \"type_modifiers\": sorted(identity(code) for code in modifiers),\n", - " \"regions\": sorted((identity(region[\"code\"]), sorted(identity(code) for code in region[\"modifiers\"])) for region in regions),\n", - " }, sort_keys=True)\n", - " result.append({\n", - " \"segment_number\": int(item.SegmentNumber),\n", - " \"label\": str(item.SegmentLabel),\n", - " \"description\": str(item.get(\"SegmentDescription\", \"\")),\n", - " \"algorithm_type\": str(item.SegmentAlgorithmType),\n", - " \"algorithm_name\": str(item.get(\"SegmentAlgorithmName\", \"\")),\n", - " \"category\": category, \"property_type\": property_type,\n", - " \"type_modifiers\": modifiers, \"anatomic_regions\": regions,\n", - " \"semantic_key\": semantic_key,\n", - " })\n", - " require(len({item[\"semantic_key\"] for item in result}) == len(result), \"Ambiguous duplicate semantic segment definitions\")\n", - " return sorted(result, key=lambda item: item[\"semantic_key\"])" + "client = IDCClient()\n", + "client.fetch_index(\"seg_index\")\n", + "\n", + "query = \"\"\"\n", + "SELECT\n", + " seg_series.PatientID,\n", + " seg_index.segmented_SeriesInstanceUID AS ct_series_uid,\n", + " seg_index.SeriesInstanceUID AS seg_series_uid,\n", + " ct_series.series_size_MB + seg_series.series_size_MB AS download_size_MB\n", + "FROM seg_index\n", + "JOIN index AS seg_series ON seg_index.SeriesInstanceUID = seg_series.SeriesInstanceUID\n", + "JOIN index AS ct_series ON seg_index.segmented_SeriesInstanceUID = ct_series.SeriesInstanceUID\n", + "WHERE seg_series.collection_id = 'c4kc_kits'\n", + "ORDER BY download_size_MB\n", + "LIMIT 3\n", + "\"\"\"\n", + "pairs = client.sql_query(query)\n", + "pairs" ] }, { "cell_type": "markdown", - "id": "10236b4e", + "id": "8080c031", "metadata": {}, "source": [ - "### Convert a validated pair" + "Now download both series of each pair. By default `idc-index` organizes files into a\n", + "`collection/patient/study/series` folder hierarchy. Here we use the `dirTemplate` argument to put each series\n", + "directly into a folder named by its `SeriesInstanceUID`, which makes the files easy to find afterwards." ] }, { "cell_type": "code", "execution_count": null, - "id": "cf41526b", - "metadata": { - "tags": [ - "definitions" - ] - }, + "id": "4ab233a5", + "metadata": {}, "outputs": [], "source": [ - "def validate_and_convert(ct_datasets, seg, expected_count: int, allow_omitted_empty_frames: bool = False):\n", - " \"\"\"Validate one real DICOM pair; return HU image, separate masks, and provenance.\"\"\"\n", - " ct_volume, ct_datasets = load_ct_volume(ct_datasets, expected_count)\n", - " checks = validate_seg_geometry(seg, ct_datasets, allow_omitted_empty_frames)\n", - " segments = segment_metadata(seg)\n", - " source_uids = [str(ds.SOPInstanceUID) for ds in ct_datasets]\n", - " try:\n", - " mask_volume = seg.get_volume(\n", - " segment_numbers=[item[\"segment_number\"] for item in segments],\n", - " combine_segments=False, dtype=np.uint8,\n", - " allow_missing_positions=allow_omitted_empty_frames, atol=0.01,\n", - " ).match_geometry(ct_volume, tol=1e-4)\n", - " except (ValueError, RuntimeError) as error:\n", - " raise ValidationError(f\"highdicom SEG grid alignment failed: {error}\") from error\n", - " masks_dhwk = mask_volume.array\n", - " require(masks_dhwk.shape == (*ct_volume.array.shape, len(segments)), \"Decoded SEG shape differs from CT\")\n", - " require(np.isin(masks_dhwk, (0, 1)).all(), \"Non-binary decoded mask values\")\n", - " image_hu = np.ascontiguousarray(ct_volume.array[np.newaxis], dtype=np.float32)\n", - " masks = np.ascontiguousarray(np.moveaxis(masks_dhwk, -1, 0))\n", - " affine = ct_volume.affine\n", - " checks[\"ct_volume_loader\"] = \"highdicom.get_volume_from_series\"\n", - " checks[\"hu_units_from_original_ct_rule\"] = any(\"RescaleType\" not in ds for ds in ct_datasets)\n", - " checks[\"seg_volume_loader\"] = \"highdicom.seg.Segmentation.get_volume + Volume.match_geometry\"\n", - " checks[\"grid_match_tolerance_voxels\"] = 1e-4\n", - " first = ct_datasets[0]\n", - " metadata = {\n", - " \"patient_id\": str(first.PatientID),\n", - " \"study_uid\": str(first.StudyInstanceUID),\n", - " \"ct_series_uid\": str(first.SeriesInstanceUID),\n", - " \"seg_series_uid\": str(seg.SeriesInstanceUID),\n", - " \"seg_sop_instance_uid\": str(seg.SOPInstanceUID),\n", - " \"frame_of_reference_uid\": str(first.FrameOfReferenceUID),\n", - " \"source_sop_instance_uids_in_slice_order\": source_uids,\n", - " \"array_axes\": [\"slice\", \"row\", \"column\"],\n", - " \"world_coordinate_system\": \"DICOM LPS, millimeters\",\n", - " \"affine_dhw_to_lps_mm\": affine.tolist(),\n", - " \"image_units\": \"HU\", \"segments\": segments, \"checks\": checks,\n", - " \"image_shape\": list(image_hu.shape), \"mask_shape\": list(masks.shape),\n", - " \"nonzero_voxels_by_channel\": [int(np.count_nonzero(mask)) for mask in masks],\n", - " }\n", - " return image_hu, masks, metadata" + "download_dir = Path(\"idc_data\")\n", + "\n", + "client.download_from_selection(\n", + " downloadDir=str(download_dir),\n", + " seriesInstanceUID=pairs.ct_series_uid.tolist() + pairs.seg_series_uid.tolist(),\n", + " dirTemplate=\"%SeriesInstanceUID\",\n", + ")" ] }, { "cell_type": "markdown", - "id": "1152e91f", + "id": "56e625d5", "metadata": {}, "source": [ - "### Patient grouping and source checksums" + "## 2. Load a CT series as a 3D volume\n", + "\n", + "A CT series is stored as a set of single-frame DICOM files, one per slice. To turn these into a 3D array, the\n", + "slices must be sorted by their position in space (the file names and `InstanceNumber` are not reliable for this),\n", + "and the stored pixel values must be converted to Hounsfield units using `RescaleSlope` and `RescaleIntercept`.\n", + "\n", + "`highdicom.get_volume_from_series()` does both of these things for us. It returns a `highdicom.Volume`, which\n", + "holds:\n", + "\n", + "- `array`: the voxel values, as a NumPy array with axes `(slices, rows, columns)`\n", + "- `affine`: a 4×4 matrix mapping voxel indices to positions in the patient coordinate system (in mm)\n", + "\n", + "We will work with the first pair for now." ] }, { "cell_type": "code", "execution_count": null, - "id": "286f15e7", - "metadata": { - "tags": [ - "definitions" - ] - }, + "id": "b058c984", + "metadata": {}, "outputs": [], "source": [ - "def patient_disjoint_split(records, seed: int = 2026, test_fraction: float = 0.25):\n", - " \"\"\"Assign all cases for each collection/patient group to the same partition.\"\"\"\n", - " require(0 < test_fraction < 1, \"test_fraction must be between 0 and 1\")\n", - " for record in records:\n", - " require(bool(record.get(\"collection_id\")) and bool(record.get(\"patient_id\")), \"Missing patient group identity\")\n", - " groups = sorted({(record[\"collection_id\"], record[\"patient_id\"]) for record in records})\n", - " require(len(groups) >= 2, \"At least two patient groups are needed\")\n", - " random.Random(seed).shuffle(groups)\n", - " n_test = min(len(groups) - 1, max(1, math.ceil(len(groups) * test_fraction)))\n", - " test_groups = set(groups[:n_test])\n", - " result = [dict(record, split=\"test\" if (record[\"collection_id\"], record[\"patient_id\"]) in test_groups else \"train\") for record in records]\n", - " train = {(r[\"collection_id\"], r[\"patient_id\"]) for r in result if r[\"split\"] == \"train\"}\n", - " test = {(r[\"collection_id\"], r[\"patient_id\"]) for r in result if r[\"split\"] == \"test\"}\n", - " require(train and test and train.isdisjoint(test), \"Patient split validation failed\")\n", - " return result\n", - "\n", - "\n", - "def file_sha256(path: Path) -> str:\n", - " digest = hashlib.sha256()\n", - " with path.open(\"rb\") as stream:\n", - " for chunk in iter(lambda: stream.read(1024 * 1024), b\"\"):\n", - " digest.update(chunk)\n", - " return digest.hexdigest()" + "case = pairs.iloc[0]\n", + "\n", + "ct_files = sorted((download_dir / case.ct_series_uid).glob(\"*.dcm\"))\n", + "ct_datasets = [pydicom.dcmread(f) for f in ct_files]\n", + "ct_volume = hd.get_volume_from_series(ct_datasets, dtype=np.float32)\n", + "\n", + "print(\"CT shape:\", ct_volume.shape)\n", + "print(\"Voxel spacing (mm):\", ct_volume.spacing)\n", + "print(\"HU range:\", ct_volume.array.min(), \"to\", ct_volume.array.max())\n", + "print(\"Affine:\\n\", ct_volume.affine)" ] }, { "cell_type": "markdown", - "id": "1a6c69cc", + "id": "e409bb10", "metadata": {}, "source": [ - "### 2. Load a pinned metadata catalog\n", - "\n", - "The main catalog has one row per series. `seg_index` supplies the referenced\n", - "CT series, and `volume_geometry_index` supplies precomputed regular-grid checks. We compare selected cloud object identifiers with the download\n", - "client's package index, so a catalog/client mismatch fails before downloading.\n", - "The Hugging Face option downloads three Parquet metadata files and records their\n", - "SHA256 checksums; this can be larger than the small DICOM example itself." + "## 3. Load the segmentation and align it with the CT\n", + "\n", + "A DICOM SEG is a single multi-frame file. that can contain several *segments* (one per structure).\n", + "There are two variants: in BINARY SEGs (like those in C4KC-KiTS), each frame stores a binary mask\n", + "for one segment on one slice, and for LABELMAP SEGs, segments are combined together in a single\n", + "frame for a given slice. Highdicom's `highdicom.seg.Segmentation` class works with both types and\n", + "gives consistent results regardless of how frames are stored, hiding this difference from you.\n", + "Let's read the file with `highdicom.seg.segread()` and see which segments it contains." ] }, { "cell_type": "code", "execution_count": null, - "id": "e7480611", + "id": "57f4d9b2", "metadata": {}, "outputs": [], "source": [ - "client = IDCClient()\n", - "require(client.get_idc_version() == \"v25\", \"This example is pinned to IDC v25\")\n", - "catalog_paths = {}\n", - "if CATALOG_SOURCE == \"huggingface\":\n", - " cache_dir = WORK_DIR / \"catalog\" / HF_REVISION\n", - " cache_dir.mkdir(parents=True, exist_ok=True)\n", - " for name in (\"idc_index\", \"seg_index\", \"volume_geometry_index\"):\n", - " path = cache_dir / f\"{name}.parquet\"\n", - " if not path.exists():\n", - " url = f\"https://huggingface.co/datasets/ImagingDataCommons/idc-index-data/resolve/{HF_REVISION}/{name}.parquet\"\n", - " temporary = path.with_suffix(\".partial\")\n", - " with requests.get(url, stream=True, timeout=(15, 90)) as response:\n", - " response.raise_for_status()\n", - " with temporary.open(\"wb\") as stream:\n", - " for chunk in response.iter_content(1024 * 1024):\n", - " stream.write(chunk)\n", - " temporary.replace(path)\n", - " catalog_paths[name] = path\n", - " index_table = pd.read_parquet(catalog_paths[\"idc_index\"])\n", - " seg_table = pd.read_parquet(catalog_paths[\"seg_index\"])\n", - " geometry_table = pd.read_parquet(catalog_paths[\"volume_geometry_index\"])\n", - "elif CATALOG_SOURCE == \"idc-index\":\n", - " client.fetch_index(\"seg_index\")\n", - " client.fetch_index(\"volume_geometry_index\")\n", - " require(getattr(client, \"seg_index\", None) is not None, \"SEG index download failed\")\n", - " require(getattr(client, \"volume_geometry_index\", None) is not None, \"Geometry index download failed\")\n", - " index_table, seg_table, geometry_table = client.index, client.seg_index, client.volume_geometry_index\n", - " catalog_paths = {\"idc_index\": Path(client.indices_overview[\"index\"][\"file_path\"]),\n", - " \"seg_index\": Path(client.indices_overview[\"seg_index\"][\"file_path\"]),\n", - " \"volume_geometry_index\": Path(client.indices_overview[\"volume_geometry_index\"][\"file_path\"])}\n", - "else:\n", - " raise ValidationError(\"Unknown CATALOG_SOURCE\")\n", - "\n", - "catalog_provenance = {}\n", - "for name, path in catalog_paths.items():\n", - " generated_version = pq.read_metadata(path).metadata[b\"idc_index_data_version\"].decode()\n", - " require(generated_version.split(\"-\")[0] == version(\"idc-index-data\"), \"Catalog release differs from the download client\")\n", - " catalog_provenance[name] = {\"sha256\": file_sha256(path), \"generated_version\": generated_version,\n", - " \"size_bytes\": path.stat().st_size}\n", - "print(f\"Catalog source: {CATALOG_SOURCE}; IDC release: {client.get_idc_version()}\")" + "seg_file = next((download_dir / case.seg_series_uid).glob(\"*.dcm\"))\n", + "seg = hd.seg.segread(seg_file)\n", + "\n", + "for segment in seg.SegmentSequence:\n", + " print(segment.SegmentNumber, segment.SegmentLabel)" ] }, { "cell_type": "markdown", - "id": "2a5f05b7", + "id": "72718f2c", "metadata": {}, "source": [ - "### 3. Select four referenced pairs and inspect the download estimate\n", - "\n", - "Restrict both modalities to this collection and the documented CC BY 3.0 license.\n", - "Require one SEG instance per series and a CT marked `regularly_spaced_3d_volume`\n", - "in the geometry index, retain one pair per patient, and show the\n", - "estimate before any DICOM download. Raise an error if the configured size cap\n", - "or patient count cannot be met. Check the [collection's usage and citation requirements](https://www.cancerimagingarchive.net/collection/c4kc-kits/)." + "Each segment also carries coded descriptions of what it represents (for example\n", + "`SegmentedPropertyTypeCodeSequence`). These coded values are more reliable than the free-text `SegmentLabel`\n", + "when combining segmentations from different sources.\n", + "\n", + "We now need the masks as an array with the same shape and voxel ordering as the CT. A SEG does not have to store\n", + "its frames in the same order as the CT slices, and it may omit frames that are empty. Two `highdicom` methods\n", + "handle this for us:\n", + "\n", + "- `seg.get_volume()` assembles the frames into a 3D volume. With `combine_segments=False`, the result has one\n", + " binary channel per segment in a trailing axis: `(slices, rows, columns, segments)`. (Use `combine_segments=True`\n", + " instead if you want a single label map in which each voxel holds a segment number.)\n", + "- `.match_geometry(ct_volume)` rearranges the mask volume so that it lies on exactly the same voxel grid as the\n", + " CT. This only pads, flips or permutes the array; it never interpolates, and it raises an error if the two\n", + " volumes cannot be aligned.\n", + "\n", + "We pass `tol=1e-4` (in units of voxels) to `match_geometry()` because the SEG files in this collection store\n", + "positions and pixel spacing with fewer decimal places than the CT files, so the two grids differ by a tiny\n", + "rounding error." ] }, { "cell_type": "code", "execution_count": null, - "id": "2e65c739", + "id": "54bb6587", "metadata": {}, "outputs": [], "source": [ - "cohort = index_table.loc[index_table.collection_id.eq(COLLECTION)].copy()\n", - "pairs = select_catalog_pairs(index_table, seg_table, geometry_table, COLLECTION, N_PATIENTS, MAX_DICOM_MB)\n", - "estimated_mb = float(pairs.estimated_mb.sum())\n", - "require(estimated_mb <= MAX_DICOM_MB, \"Cohort exceeds MAX_DICOM_MB; reduce N_PATIENTS\")\n", - "selected_uids = pairs.SeriesInstanceUID_ct.tolist() + pairs.SeriesInstanceUID_seg.tolist()\n", - "selected_catalog = cohort.loc[cohort.SeriesInstanceUID.isin(selected_uids)]\n", - "client_rows = client.index.set_index(\"SeriesInstanceUID\")\n", - "for row in selected_catalog.to_dict(\"records\"):\n", - " require(row[\"SeriesInstanceUID\"] in client_rows.index, \"Selected series missing from download client\")\n", - " download_row = client_rows.loc[row[\"SeriesInstanceUID\"]]\n", - " require(str(download_row.crdc_series_uuid) == str(row[\"crdc_series_uuid\"])\n", - " and int(download_row.instanceCount) == int(row[\"instanceCount\"]), \"Catalog/download client mismatch\")\n", - "selected_catalog.to_json(OUTPUT_DIR / \"cohort_catalog.json\", orient=\"records\", indent=2)\n", - "print(f\"DICOM estimate: {estimated_mb:.2f} MB for {len(pairs)} patients; cap {MAX_DICOM_MB:.0f} MB\")\n", - "display(pairs[[\"PatientID_ct\", \"instanceCount_ct\", \"estimated_mb\", \"license_short_name_ct\"]]\n", - " .rename(columns={\"PatientID_ct\": \"Patient\", \"instanceCount_ct\": \"CT slices\", \"estimated_mb\": \"CT + SEG (MB)\", \"license_short_name_ct\": \"License\"}))" + "mask_volume = seg.get_volume(combine_segments=False).match_geometry(ct_volume, tol=1e-4)\n", + "\n", + "print(\"CT shape: \", ct_volume.shape)\n", + "print(\"Mask shape:\", mask_volume.shape)" ] }, { "cell_type": "markdown", - "id": "2bb962c5", + "id": "8337802a", "metadata": {}, "source": [ - "### 4. Download, validate and export one patient at a time\n", - "\n", - "Every file's SHA256 is recorded. The catalog's CT instance count must match the\n", - "downloaded series and the SEG's global source references. All source/segment\n", - "frames must be present by default; zero-filled omitted frames require an explicit\n", - "dataset assumption, controlled by `ALLOW_OMITTED_EMPTY_FRAMES`.\n", - "\n", - "The preview uses a fixed HU display window and the slice with the largest total\n", - "segmented area. This changes only the display, not stored HU values. No patient\n", - "identifiers are drawn into the image. Each exported mask channel has its original\n", - "segment metadata; DICOM algorithm type is preserved without reinterpretation." + "Let's check the alignment by overlaying each segment on the CT slice where the segmentation is largest. We use a\n", + "soft-tissue display window (-200 to 300 HU) for viewing; this does not change the underlying values." ] }, { "cell_type": "code", "execution_count": null, - "id": "528ac35b", + "id": "305dff4c", "metadata": {}, "outputs": [], "source": [ - "records = []\n", - "canonical_channel_keys = None\n", - "for row in pairs.to_dict(\"records\"):\n", - " patient = row[\"PatientID_ct\"]\n", - " ct_uid, seg_uid = row[\"SeriesInstanceUID_ct\"], row[\"SeriesInstanceUID_seg\"]\n", - " client.download_from_selection(downloadDir=str(DICOM_DIR), seriesInstanceUID=[ct_uid, seg_uid],\n", - " show_progress_bar=False, use_s5cmd_sync=True)\n", - " study_dir = DICOM_DIR / COLLECTION / patient / row[\"StudyInstanceUID_ct\"]\n", - " ct_paths = sorted((study_dir / f\"CT_{ct_uid}\").glob(\"*.dcm\"))\n", - " seg_paths = sorted((study_dir / f\"SEG_{seg_uid}\").glob(\"*.dcm\"))\n", - " require(len(seg_paths) == 1, \"Expected exactly one downloaded SEG instance\")\n", - " ct_datasets = [pydicom.dcmread(path) for path in ct_paths]\n", - " seg = hd.seg.segread(seg_paths[0])\n", - " require(all(str(ds.SeriesInstanceUID) == ct_uid and str(ds.PatientID) == patient for ds in ct_datasets), \"Downloaded CT differs from selection\")\n", - " require(str(seg.SeriesInstanceUID) == seg_uid, \"Downloaded SEG differs from selection\")\n", - " image_hu, masks, metadata = validate_and_convert(ct_datasets, seg, int(row[\"instanceCount_ct\"]),\n", - " ALLOW_OMITTED_EMPTY_FRAMES)\n", - " channel_keys = [item[\"semantic_key\"] for item in metadata[\"segments\"]]\n", - " if canonical_channel_keys is None:\n", - " canonical_channel_keys = channel_keys\n", - " require(channel_keys == canonical_channel_keys, \"Cases have different semantic channel definitions\")\n", - " sample_id = hashlib.sha256(ct_uid.encode()).hexdigest()[:16]\n", - " tensor_file = f\"{sample_id}.pt\"\n", - " torch.save({\"image_hu\": torch.from_numpy(image_hu), \"masks\": torch.from_numpy(masks)}, OUTPUT_DIR / tensor_file)\n", - " metadata[\"collection_id\"] = COLLECTION\n", - " metadata[\"catalog_regularly_spaced_3d_volume\"] = bool(row[\"regularly_spaced_3d_volume\"])\n", - " metadata[\"sources\"] = [{\"path\": str(path.relative_to(DICOM_DIR)), \"sha256\": file_sha256(path),\n", - " \"size_bytes\": path.stat().st_size} for path in ct_paths + seg_paths]\n", - " metadata[\"license_ct\"] = row[\"license_short_name_ct\"]\n", - " metadata[\"license_seg\"] = row[\"license_short_name_seg\"]\n", - " metadata[\"source_dois\"] = sorted({row[\"source_DOI_ct\"], row[\"source_DOI_seg\"]})\n", - " metadata_file = f\"{sample_id}.json\"\n", - " (OUTPUT_DIR / metadata_file).write_text(json.dumps(metadata, indent=2))\n", - " records.append({\"collection_id\": COLLECTION, \"patient_id\": patient,\n", - " \"ct_series_uid\": ct_uid, \"seg_series_uid\": seg_uid,\n", - " \"tensor_file\": tensor_file, \"metadata_file\": metadata_file,\n", - " \"ct_slices\": image_hu.shape[1], \"image_shape\": list(image_hu.shape),\n", - " \"mask_shape\": list(masks.shape), \"nonzero_voxels\": metadata[\"nonzero_voxels_by_channel\"],\n", - " \"dicom_bytes\": sum(item[\"size_bytes\"] for item in metadata[\"sources\"])})\n", - " slice_index = int(np.argmax(masks.sum(axis=(0, 2, 3))))\n", - " fig, axes = plt.subplots(1, len(masks) + 1, figsize=(5 * (len(masks) + 1), 4))\n", - " axes[0].imshow(image_hu[0, slice_index], cmap=\"gray\", vmin=-200, vmax=300)\n", - " axes[0].set_title(f\"CT (HU); slice {slice_index}\")\n", - " for axis, mask, segment in zip(axes[1:], masks, metadata[\"segments\"]):\n", - " axis.imshow(image_hu[0, slice_index], cmap=\"gray\", vmin=-200, vmax=300)\n", - " axis.imshow(np.ma.masked_equal(mask[slice_index], 0), cmap=\"autumn\", alpha=0.45, vmin=0, vmax=1)\n", - " axis.set_title(f\"{segment['label']} — {segment['algorithm_type']}\")\n", - " for axis in axes:\n", - " axis.axis(\"off\")\n", - " fig.suptitle(f\"Validated CT/SEG pair · {sample_id}\")\n", - " fig.tight_layout()\n", - " fig.savefig(OUTPUT_DIR / f\"{sample_id}_overlay.png\", dpi=120)\n", - " plt.show()\n", - " plt.close(fig)\n", - " print(f\"Validated {patient}: {image_hu.shape} image; {masks.shape} masks; omitted frames {metadata['checks']['omitted_source_segment_frames']}\")\n", - " del ct_datasets, seg, image_hu, masks" + "masks = mask_volume.array\n", + "slice_index = masks.sum(axis=(1, 2, 3)).argmax()\n", + "\n", + "fig, axes = plt.subplots(1, masks.shape[-1], figsize=(10, 5))\n", + "for k, ax in enumerate(axes):\n", + " ax.imshow(ct_volume.array[slice_index], cmap=\"gray\", vmin=-200, vmax=300)\n", + " ax.imshow(np.ma.masked_equal(masks[slice_index, :, :, k], 0), cmap=\"autumn\", alpha=0.5)\n", + " ax.set_title(seg.SegmentSequence[k].SegmentLabel)\n", + " ax.axis(\"off\")\n", + "plt.show()" ] }, { "cell_type": "markdown", - "id": "19a5fc07", + "id": "2325769f", "metadata": {}, "source": [ - "### 5. Write and verify the patient-disjoint split\n", + "## 4. Convert to PyTorch tensors\n", "\n", - "Grouping uses `(collection_id, PatientID)`, since PatientID is scoped to a\n", - "collection. The fixed seed is applied to a sorted patient list, so input row\n", - "order does not affect assignments. All cases for a patient share one partition.\n", - "This grouping cannot identify aliases or duplicate people across collections." + "PyTorch 3D models (e.g. `torch.nn.Conv3d`) expect inputs laid out as `(channels, depth, height, width)`.\n", + "\n", + "- The CT Volume has no channel dimension, but we want a pytorch tensor with a single channel, as is the standard format for grayscale images, so we add a leading axis: `(1, D, H, W)`.\n", + "- The mask Volume has one channel per segment in the last axis, so we move that axis to the front: `(K, D, H, W)`. Note that if you had used `combine_segments=True` above, the mask Volume would have no channels, and we would instead treat it just like we treat the CT volume." ] }, { "cell_type": "code", "execution_count": null, - "id": "280de258", + "id": "86f8c8b6", "metadata": {}, "outputs": [], "source": [ - "split_records = patient_disjoint_split(records, SPLIT_SEED, TEST_FRACTION)\n", - "manifest = {\n", - " \"idc_release\": client.get_idc_version(), \"catalog_source\": CATALOG_SOURCE,\n", - " \"hf_revision\": HF_REVISION if CATALOG_SOURCE == \"huggingface\" else None,\n", - " \"catalog_provenance\": catalog_provenance,\n", - " \"packages\": {name: version(name) for name in (\"idc-index\", \"idc-index-data\", \"highdicom\", \"pydicom\", \"numpy\", \"pandas\", \"torch\", \"matplotlib\")},\n", - " \"split_seed\": SPLIT_SEED, \"test_fraction\": TEST_FRACTION,\n", - " \"sampling\": \"smallest referenced regular-grid CT/SEG pair per patient; tutorial convenience sample\",\n", - " \"geometry_preflight\": \"volume_geometry_index.regularly_spaced_3d_volume\",\n", - " \"allow_omitted_empty_frames\": ALLOW_OMITTED_EMPTY_FRAMES,\n", - " \"channel_semantic_keys\": canonical_channel_keys,\n", - " \"estimated_dicom_mb\": estimated_mb,\n", - " \"actual_dicom_bytes\": sum(record[\"dicom_bytes\"] for record in records),\n", - " \"runtime_seconds\": round(perf_counter() - run_started, 2),\n", - " \"records\": split_records,\n", - "}\n", - "(OUTPUT_DIR / \"split_manifest.json\").write_text(json.dumps(manifest, indent=2))\n", - "display(pd.DataFrame(split_records)[[\"patient_id\", \"split\", \"ct_slices\", \"image_shape\", \"mask_shape\", \"nonzero_voxels\"]])\n", - "print(f\"Actual selected DICOM: {manifest['actual_dicom_bytes'] / 1e6:.2f} MB; runtime {manifest['runtime_seconds']:.2f} seconds\")" + "image = torch.from_numpy(ct_volume.array).unsqueeze(0)\n", + "target = torch.from_numpy(masks).permute(3, 0, 1, 2)\n", + "\n", + "print(\"Image:\", tuple(image.shape), image.dtype)\n", + "print(\"Masks:\", tuple(target.shape), target.dtype)" ] }, { "cell_type": "markdown", - "id": "42271c87", + "id": "a97bbc63", "metadata": {}, "source": [ - "### 6. Consume saved tensors with PyTorch\n", + "## 5. Put it together in a `Dataset`\n", + "\n", + "Finally, we combine the steps above into a function and use it in a minimal `torch.utils.data.Dataset` that loads\n", + "one CT/SEG pair per item. We also scale the HU values to the range [0, 1] using the same window as above, which is\n", + "a simple and common choice of intensity normalization for CT.\n", "\n", - "This Dataset opens one case at a time. Batch size 1 accommodates different\n", - "volume depths. The optional fixed window `[-200, 300] HU -> [0, 1]` is an\n", - "illustrative input transform, not a recommended training configuration or a\n", - "transformation fitted on the test set. Stored tensors always retain HU.\n", - "Only load `.pt` files generated by this notebook, using `weights_only=True`." + "Since each CT has a different number of slices, we use `batch_size=1`. In a real training pipeline you would\n", + "typically add resampling, cropping or patch extraction here (for example using\n", + "[MONAI](https://monai.io/) transforms) so that cases can be batched together. For larger datasets, you may also\n", + "want to save the tensors to disk once (e.g. with `torch.save`) rather than parsing the DICOM files every epoch." ] }, { "cell_type": "code", "execution_count": null, - "id": "49a882e8", + "id": "de0ef3e9", "metadata": {}, "outputs": [], "source": [ + "def load_case(ct_dir: Path, seg_file: Path) -> tuple[torch.Tensor, torch.Tensor]:\n", + " \"\"\"Load a CT series and its SEG as (1, D, H, W) and (K, D, H, W) tensors.\"\"\"\n", + " ct_volume = hd.get_volume_from_series(\n", + " [pydicom.dcmread(f) for f in sorted(ct_dir.glob(\"*.dcm\"))],\n", + " dtype=np.float32,\n", + " )\n", + " seg = hd.seg.segread(seg_file)\n", + " mask_volume = seg.get_volume(combine_segments=False).match_geometry(ct_volume, tol=1e-4)\n", + "\n", + " image = torch.from_numpy(ct_volume.array).unsqueeze(0)\n", + " masks = torch.from_numpy(mask_volume.array).permute(3, 0, 1, 2)\n", + " return image, masks\n", + "\n", + "\n", "class CTSegDataset(torch.utils.data.Dataset):\n", - " def __init__(self, output_dir, partition, normalize=False):\n", - " require(partition in (\"train\", \"test\"), \"Unknown dataset partition\")\n", - " self.output_dir = Path(output_dir)\n", - " saved_manifest = json.loads((self.output_dir / \"split_manifest.json\").read_text())\n", - " self.records = [record for record in saved_manifest[\"records\"] if record[\"split\"] == partition]\n", - " self.normalize = normalize\n", + " def __init__(self, pairs, download_dir: Path):\n", + " self.pairs = pairs.reset_index(drop=True)\n", + " self.download_dir = download_dir\n", "\n", " def __len__(self):\n", - " return len(self.records)\n", - "\n", - " def __getitem__(self, index):\n", - " record = self.records[index]\n", - " tensors = torch.load(self.output_dir / record[\"tensor_file\"], map_location=\"cpu\", weights_only=True)\n", - " image = tensors[\"image_hu\"]\n", - " if self.normalize:\n", - " image = (image.clamp(-200, 300) + 200) / 500\n", - " return {\"image\": image, \"masks\": tensors[\"masks\"], \"patient_id\": record[\"patient_id\"]}\n", - "\n", - "train_dataset = CTSegDataset(OUTPUT_DIR, \"train\", normalize=True)\n", - "test_dataset = CTSegDataset(OUTPUT_DIR, \"test\", normalize=True)\n", - "loader = torch.utils.data.DataLoader(train_dataset, batch_size=1, shuffle=False, num_workers=0)\n", - "batch = next(iter(loader))\n", - "print(\"Train/test patients:\", len(train_dataset), len(test_dataset))\n", - "print(\"Batch shapes:\", tuple(batch[\"image\"].shape), tuple(batch[\"masks\"].shape))" - ] - }, - { - "cell_type": "markdown", - "id": "e23bf602", - "metadata": {}, - "source": [ - "## Checks\n", - "\n", - "Reload saved files and verify shapes, binary values, coded channel consistency,\n", - "source count, and patient separation. These checks demonstrate data preparation;\n", - "they do not establish annotation quality or predict segmentation accuracy." + " return len(self.pairs)\n", + "\n", + " def __getitem__(self, i):\n", + " row = self.pairs.iloc[i]\n", + " ct_dir = self.download_dir / row.ct_series_uid\n", + " seg_file = next((self.download_dir / row.seg_series_uid).glob(\"*.dcm\"))\n", + " image, masks = load_case(ct_dir, seg_file)\n", + " image = (image.clamp(-200, 300) + 200) / 500 # window and scale to [0, 1]\n", + " return image, masks" ] }, { "cell_type": "code", "execution_count": null, - "id": "4fb10c3c", + "id": "31c213a4", "metadata": {}, "outputs": [], "source": [ - "saved = json.loads((OUTPUT_DIR / \"split_manifest.json\").read_text())\n", - "train_groups = {(r[\"collection_id\"], r[\"patient_id\"]) for r in saved[\"records\"] if r[\"split\"] == \"train\"}\n", - "test_groups = {(r[\"collection_id\"], r[\"patient_id\"]) for r in saved[\"records\"] if r[\"split\"] == \"test\"}\n", - "require(train_groups and test_groups and train_groups.isdisjoint(test_groups), \"Saved split leaks patients\")\n", - "for record in saved[\"records\"]:\n", - " tensors = torch.load(OUTPUT_DIR / record[\"tensor_file\"], map_location=\"cpu\", weights_only=True)\n", - " metadata = json.loads((OUTPUT_DIR / record[\"metadata_file\"]).read_text())\n", - " require(list(tensors[\"image_hu\"].shape) == record[\"image_shape\"], \"Saved image shape mismatch\")\n", - " require(list(tensors[\"masks\"].shape) == record[\"mask_shape\"], \"Saved mask shape mismatch\")\n", - " require(tensors[\"image_hu\"].dtype == torch.float32 and torch.isfinite(tensors[\"image_hu\"]).all().item(), \"Invalid HU tensor\")\n", - " require(tensors[\"masks\"].dtype == torch.uint8 and torch.all((tensors[\"masks\"] == 0) | (tensors[\"masks\"] == 1)).item(), \"Invalid binary masks\")\n", - " require([segment[\"semantic_key\"] for segment in metadata[\"segments\"]] == saved[\"channel_semantic_keys\"], \"Saved channels differ between cases\")\n", - " require(len(metadata[\"source_sop_instance_uids_in_slice_order\"]) == record[\"ct_slices\"], \"Saved source count mismatch\")\n", - " require(metadata[\"catalog_regularly_spaced_3d_volume\"] is True, \"Missing catalog geometry preflight\")\n", - " del tensors\n", - "print(\"PASS: validated real CT/SEG pairs, consistent tensor channels, and disjoint patient groups\")" + "dataset = CTSegDataset(pairs, download_dir)\n", + "loader = torch.utils.data.DataLoader(dataset, batch_size=1, shuffle=True)\n", + "\n", + "for image, masks in loader:\n", + " print(\"Image batch:\", tuple(image.shape), \" Mask batch:\", tuple(masks.shape))" ] }, { "cell_type": "markdown", - "id": "3d0ac3a8", + "id": "1e4d45a8", "metadata": {}, "source": [ - "## Next Steps\n", + "Note that this assumes every SEG in the dataset contains the same segments in the same order (here, kidney then\n", + "mass), which is true for C4KC-KiTS. When combining data from different collections, check the segment\n", + "descriptions and reorder the channels so that each channel always means the same structure.\n", "\n", - "Use a larger, representative cohort before training or assessing a model.\n", - "Keep patient grouping when adding studies, slices or patches. Fit any learned\n", - "preprocessing on training patients only. Validate resampling and class semantics\n", - "explicitly before extending this example to other collections or modalities.\n", + "## Next steps\n", "\n", - "Outputs are under `idc_ct_seg_data/outputs`: `.pt` tensors, per-case JSON geometry\n", - "and source checksums, overlay PNGs, the selected catalog rows, and the split\n", - "manifest. The four-patient split is an educational artifact, not a benchmark.\n", + "- Select a larger cohort by changing the `LIMIT` in the query, or explore other collections with segmentations\n", + " using `seg_index`.\n", + "- If you will train and evaluate a model, split the data by patient (not by series) so that no patient\n", + " appears in both the training and test sets.\n", + "- See the [highdicom documentation](https://highdicom.readthedocs.io/en/latest/volume.html) for more on working\n", + " with volumes, including resampling and reorienting them.\n", "\n", "## Support\n", "\n", - "For questions, use the [IDC user forum](https://discourse.canceridc.dev) or open\n", - "an issue in [IDC-Tutorials](https://github.com/ImagingDataCommons/IDC-Tutorials/issues).\n", - "More examples are available in the [IDC-Tutorials repository](https://github.com/ImagingDataCommons/IDC-Tutorials).\n", + "For questions, use the [IDC user forum](https://discourse.canceridc.dev) or open an issue in\n", + "[IDC-Tutorials](https://github.com/ImagingDataCommons/IDC-Tutorials/issues).\n", "\n", "## Acknowledgments\n", "\n", - "Imaging Data Commons has been funded in whole or in part with Federal funds from\n", - "the National Cancer Institute, National Institutes of Health, under Task Order\n", - "No. HHSN26110071 under Contract No. HHSN261201500003I.\n", - "\n", - "If you use IDC in research, cite Fedorov et al., *National Cancer Institute Imaging\n", - "Data Commons: Toward Transparency, Reproducibility, and Scalability in Imaging\n", - "Artificial Intelligence*, RadioGraphics (2023), [doi:10.1148/rg.230180](https://doi.org/10.1148/rg.230180).\n", + "Imaging Data Commons has been funded in whole or in part with Federal funds from the National Cancer Institute,\n", + "National Institutes of Health, under Task Order No. HHSN26110071 under Contract No. HHSN261201500003I.\n", "\n", - "The C4KC-KiTS data are provided under CC BY 3.0. Cite Heller et al., *Data from\n", - "C4KC-KiTS*, The Cancer Imaging Archive (2019), [doi:10.7937/TCIA.2019.IX49E8NX](https://doi.org/10.7937/TCIA.2019.IX49E8NX),\n", - "and follow the [collection's full citation instructions](https://www.cancerimagingarchive.net/collection/c4kc-kits/),\n", - "including Heller et al., *The KiTS19 Challenge Data: 300 Kidney Tumor Cases with\n", - "Clinical Context, CT Semantic Segmentations, and Surgical Outcomes*,\n", - "[arXiv:1904.00445](https://arxiv.org/abs/1904.00445).\n", + "If you use IDC in research, cite Fedorov et al., *National Cancer Institute Imaging Data Commons: Toward\n", + "Transparency, Reproducibility, and Scalability in Imaging Artificial Intelligence*, RadioGraphics (2023),\n", + "[doi:10.1148/rg.230180](https://doi.org/10.1148/rg.230180).\n", "\n", - "Libraries used: [idc-index](https://github.com/ImagingDataCommons/idc-index),\n", - "[highdicom](https://highdicom.readthedocs.io/en/latest/seg.html),\n", - "[pydicom](https://pydicom.github.io/), and [PyTorch](https://pytorch.org/)." + "The C4KC-KiTS data are provided under CC BY 3.0. Cite Heller et al., *Data from C4KC-KiTS*, The Cancer Imaging\n", + "Archive (2019), [doi:10.7937/TCIA.2019.IX49E8NX](https://doi.org/10.7937/TCIA.2019.IX49E8NX), and follow the\n", + "[collection's full citation instructions](https://www.cancerimagingarchive.net/collection/c4kc-kits/)." ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", "name": "python", - "version": "3.12" + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.2" } }, "nbformat": 4, From ca095201328d2a106ad796d6a3763954bb0c3534 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Wed, 21 Oct 2026 16:50:21 -0400 Subject: [PATCH 4/5] Remove outdated test file --- test/src/test_ct_seg_to_tensors.py | 297 ----------------------------- 1 file changed, 297 deletions(-) delete mode 100644 test/src/test_ct_seg_to_tensors.py diff --git a/test/src/test_ct_seg_to_tensors.py b/test/src/test_ct_seg_to_tensors.py deleted file mode 100644 index e4ff02c..0000000 --- a/test/src/test_ct_seg_to_tensors.py +++ /dev/null @@ -1,297 +0,0 @@ -"""Behavior checks for the self-contained CT/SEG tutorial definitions. - -Run after installing the notebook's dependencies: - python -m pytest test/src/test_ct_seg_to_tensors.py -""" -import copy -from pathlib import Path - -import highdicom as hd -import nbformat -import numpy as np -import pandas as pd -import pydicom -import pytest -from pydicom.dataset import Dataset, FileMetaDataset -from pydicom.uid import CTImageStorage, ExplicitVRLittleEndian, generate_uid - - -NOTEBOOK = Path(__file__).resolve().parents[2] / "notebooks/advanced_topics/ct_seg_to_tensors.ipynb" -namespace = {} -for cell in nbformat.read(NOTEBOOK, as_version=4).cells: - if "definitions" in cell.metadata.get("tags", []): - exec(compile(cell.source, str(NOTEBOOK), "exec"), namespace) -convert = namespace["validate_and_convert"] -split_patients = namespace["patient_disjoint_split"] -ValidationError = namespace["ValidationError"] -select_pairs = namespace["select_catalog_pairs"] - - -@pytest.fixture -def pair(): - study, series, frame_of_reference = generate_uid(), generate_uid(), generate_uid() - images = [] - for z in (0, 5, 10): - image = Dataset() - image.file_meta = FileMetaDataset() - image.file_meta.TransferSyntaxUID = ExplicitVRLittleEndian - image.SOPClassUID = CTImageStorage - image.SOPInstanceUID = generate_uid() - image.PatientID = "tutorial-patient" - image.PatientName = "Synthetic^Patient" - image.PatientBirthDate = "" - image.PatientSex = "" - image.StudyInstanceUID = study - image.SeriesInstanceUID = series - image.FrameOfReferenceUID = frame_of_reference - image.StudyID = "1" - image.SeriesNumber = 1 - image.InstanceNumber = z + 1 - image.StudyDate = "20260101" - image.StudyTime = "120000" - image.AccessionNumber = "" - image.Modality = "CT" - image.ImageOrientationPatient = [1, 0, 0, 0, 1, 0] - image.ImagePositionPatient = [0, 0, z] - image.PixelSpacing = [2, 3] - image.SliceThickness = 5 - image.Rows, image.Columns = 4, 4 - image.SamplesPerPixel = 1 - image.PhotometricInterpretation = "MONOCHROME2" - image.BitsAllocated = 16 - image.BitsStored = 16 - image.HighBit = 15 - image.PixelRepresentation = 0 - image.RescaleSlope = 2 - image.RescaleIntercept = -1000 - image.RescaleType = "HU" - image.PixelData = np.full((4, 4), z + 1, dtype=np.uint16).tobytes() - images.append(image) - masks = np.zeros((3, 4, 4, 2), dtype=np.uint8) - masks[0, 1, 2, 0] = 1 - masks[2, 2, 1, 0] = 1 - masks[2, 2, 1, 1] = 1 # Intentional overlap must survive conversion. - descriptions = [ - hd.seg.SegmentDescription( - segment_number=number, segment_label=label, - segmented_property_category=pydicom.sr.codedict.codes.SCT.AnatomicalStructure, - segmented_property_type=code, - algorithm_type=hd.seg.SegmentAlgorithmTypeValues.MANUAL, - ) - for number, label, code in [ - (1, "Kidney", pydicom.sr.codedict.codes.SCT.Kidney), - (2, "Liver", pydicom.sr.codedict.codes.SCT.Liver), - ] - ] - seg = hd.seg.Segmentation( - source_images=images, pixel_array=masks, - segmentation_type=hd.seg.SegmentationTypeValues.BINARY, - segment_descriptions=descriptions, omit_empty_frames=False, - series_instance_uid=generate_uid(), series_number=2, - sop_instance_uid=generate_uid(), instance_number=1, - manufacturer="Tutorial", manufacturer_model_name="Synthetic", - software_versions="1", device_serial_number="1", - ) - return images, copy.deepcopy(seg) - - -def test_ct_order_hu_affine_and_overlapping_mask_channels(pair): - images, seg = pair - image, masks, metadata = convert([images[2], images[0], images[1]], seg, 3) - assert image.shape == (1, 3, 4, 4) - assert image.dtype == np.float32 - # highdicom's default volume order runs from superior to inferior. - assert image[0, :, 0, 0].tolist() == [-978, -988, -998] - np.testing.assert_allclose(metadata["affine_dhw_to_lps_mm"], - [[0, 0, 3, 0], [0, 2, 0, 0], [-5, 0, 0, 10], [0, 0, 0, 1]]) - assert metadata["source_sop_instance_uids_in_slice_order"] == [ - str(ds.SOPInstanceUID) for ds in reversed(images)] - assert masks.dtype == np.uint8 - assert masks[:, 0, 2, 1].tolist() == [1, 1] - assert sorted(np.count_nonzero(mask) for mask in masks) == [1, 2] - - -@pytest.mark.parametrize("change, message", [ - ("patient", "PatientID mismatch"), - ("series", "different CT series"), - ("frame_of_reference", "FrameOfReferenceUID mismatch"), - ("position", "SEG position differs"), - ("spacing", "SEG pixel spacing differs"), - ("orientation", "SEG orientation differs"), - ("spatial_no", "does not preserve"), -]) -def test_incompatible_pair_is_rejected(pair, change, message): - images, seg = pair - if change == "patient": - seg.PatientID = "wrong-patient" - elif change == "series": - seg.ReferencedSeriesSequence[0].SeriesInstanceUID = generate_uid() - elif change == "frame_of_reference": - seg.FrameOfReferenceUID = generate_uid() - elif change == "position": - seg.PerFrameFunctionalGroupsSequence[0].PlanePositionSequence[0].ImagePositionPatient = [1, 0, 0] - elif change == "spacing": - seg.SharedFunctionalGroupsSequence[0].PixelMeasuresSequence[0].PixelSpacing = [1, 1] - elif change == "orientation": - seg.SharedFunctionalGroupsSequence[0].PlaneOrientationSequence[0].ImageOrientationPatient = [0, 1, 0, 1, 0, 0] - else: - seg.PerFrameFunctionalGroupsSequence[0].DerivationImageSequence[0].SourceImageSequence[0].SpatialLocationsPreserved = "NO" - with pytest.raises(ValidationError, match=message): - convert(images, seg, 3) - - -@pytest.mark.parametrize("change, message", [ - ("missing", "instance count"), - ("duplicate_uid", "Duplicate CT SOP"), - ("irregular", "highdicom CT volume validation failed"), - ("gantry_tilt", "highdicom CT volume validation failed"), -]) -def test_incomplete_or_unsupported_ct_is_rejected(pair, change, message): - images, seg = pair - if change == "missing": - images = images[:-1] - elif change == "duplicate_uid": - images[1].SOPInstanceUID = images[0].SOPInstanceUID - elif change == "irregular": - images[2].ImagePositionPatient = [0, 0, 11] - else: - images[2].ImagePositionPatient = [1, 0, 10] - with pytest.raises(ValidationError, match=message): - convert(images, seg, 3) - - -def test_missing_spatial_flag_is_allowed_only_after_geometry_validation(pair): - images, seg = pair - for frame in seg.PerFrameFunctionalGroupsSequence: - source = frame.DerivationImageSequence[0].SourceImageSequence[0] - if "SpatialLocationsPreserved" in source: - del source.SpatialLocationsPreserved - _, masks, metadata = convert(images, seg, 3) - assert metadata["checks"]["spatial_locations_flag_absent"] is True - assert masks[:, 0, 2, 1].tolist() == [1, 1] - - -def test_hu_values_do_not_receive_display_or_presentation_transforms(pair): - images, seg = pair - for image in images: - image.PhotometricInterpretation = "MONOCHROME1" - image.PresentationLUTShape = "INVERSE" - image.WindowCenter = -500 - image.WindowWidth = 100 - image, _, _ = convert(images, seg, 3) - assert image[0, :, 0, 0].tolist() == [-978, -988, -998] - - -def test_unspecified_rescale_units_are_not_silently_labelled_hu(pair): - images, seg = pair - del images[0].RescaleType - with pytest.raises(ValidationError, match="cannot be verified as HU"): - convert(images, seg, 3) - - -def test_original_ct_can_omit_rescale_type_under_dicom_hu_rule(pair): - images, seg = pair - for image in images: - del image.RescaleType - image.ImageType = ["ORIGINAL", "PRIMARY", "AXIAL"] - image, _, metadata = convert(images, seg, 3) - assert image[0, :, 0, 0].tolist() == [-978, -988, -998] - assert metadata["checks"]["hu_units_from_original_ct_rule"] is True - - -@pytest.mark.parametrize("image_type,multi_energy", [ - (["DERIVED", "PRIMARY", "AXIAL"], "NO"), - (["ORIGINAL", "PRIMARY", "LOCALIZER"], "NO"), - (["ORIGINAL", "PRIMARY", "AXIAL"], "YES"), -]) -def test_missing_rescale_units_are_not_inferred_outside_original_ct_rule(pair, image_type, multi_energy): - images, seg = pair - for image in images: - del image.RescaleType - image.ImageType = image_type - image.MultienergyCTAcquisition = multi_energy - with pytest.raises(ValidationError, match="cannot be verified as HU"): - convert(images, seg, 3) - - -def test_omitted_empty_seg_frames_need_explicit_policy_and_align_to_full_ct(pair): - images, seg = pair - pixels = seg.get_pixels_by_source_instance( - [str(image.SOPInstanceUID) for image in images], combine_segments=False) - sparse_seg = hd.seg.Segmentation( - source_images=images, pixel_array=pixels, - segmentation_type=hd.seg.SegmentationTypeValues.BINARY, - segment_descriptions=list(seg.SegmentSequence), omit_empty_frames=True, - series_instance_uid=generate_uid(), series_number=2, - sop_instance_uid=generate_uid(), instance_number=1, - manufacturer="Tutorial", manufacturer_model_name="Synthetic", - software_versions="1", device_serial_number="1", - ) - with pytest.raises(ValidationError, match="explicit empty-frame policy"): - convert(images, sparse_seg, 3) - image, masks, metadata = convert(images, sparse_seg, 3, allow_omitted_empty_frames=True) - assert image.shape == (1, 3, 4, 4) - assert masks.shape == (2, 3, 4, 4) - assert not masks[:, 1].any() - assert masks[:, 0, 2, 1].tolist() == [1, 1] - assert metadata["checks"]["omitted_source_segment_frames"] > 0 - - -@pytest.fixture -def catalog(): - rows, segmentations, geometries = [], [], [] - for patient, size in [("A", 1.0), ("B", 2.0), ("C", 3.0)]: - for modality in ["CT", "SEG"]: - rows.append({ - "collection_id": "c4kc_kits", "PatientID": patient, - "StudyInstanceUID": f"study-{patient}", - "SeriesInstanceUID": f"{modality}-{patient}", "Modality": modality, - "instanceCount": 3 if modality == "CT" else 1, - "series_size_MB": size, "license_short_name": "CC BY 3.0", - }) - segmentations.append({"SeriesInstanceUID": f"SEG-{patient}", - "segmented_SeriesInstanceUID": f"CT-{patient}", - "SegmentationType": "BINARY"}) - geometries.append({"SeriesInstanceUID": f"CT-{patient}", - "regularly_spaced_3d_volume": True}) - return pd.DataFrame(rows), pd.DataFrame(segmentations), pd.DataFrame(geometries) - - -@pytest.mark.parametrize("geometry_value", [False, None, pd.NA]) -def test_catalog_geometry_filter_excludes_failed_or_unknown_ct(catalog, geometry_value): - index, segmentations, geometry = catalog - geometry["regularly_spaced_3d_volume"] = geometry.regularly_spaced_3d_volume.astype(object) - geometry.loc[0, "regularly_spaced_3d_volume"] = geometry_value - pairs = select_pairs(index, segmentations, geometry, "c4kc_kits", 2, 100) - assert pairs.PatientID_ct.tolist() == ["B", "C"] - - -def test_catalog_missing_geometry_and_wrong_reference_do_not_fall_back_to_patient(catalog): - index, segmentations, geometry = catalog - geometry = geometry.iloc[1:] - segmentations.loc[1, "segmented_SeriesInstanceUID"] = "CT-A" - with pytest.raises(ValidationError, match="Not enough eligible"): - select_pairs(index, segmentations, geometry, "c4kc_kits", 2, 100) - - -def test_catalog_size_cap_fails_before_download(catalog): - with pytest.raises(ValidationError, match="exceeds MAX_DICOM_MB"): - select_pairs(*catalog, "c4kc_kits", 2, 1) - - -def test_patient_split_is_grouped_and_independent_of_record_order(): - records = [{"collection_id": "c4kc_kits", "patient_id": patient, "study": study} - for patient, study in [("A", "1"), ("A", "2"), ("B", "3"), ("C", "4"), ("D", "5")]] - first = split_patients(records, seed=2026) - second = split_patients(list(reversed(records)), seed=2026) - assignment = {record["patient_id"]: record["split"] for record in first} - assert assignment == {record["patient_id"]: record["split"] for record in second} - assert len({record["split"] for record in first if record["patient_id"] == "A"}) == 1 - assert sum(partition == "test" for partition in assignment.values()) == 1 - assert {record["patient_id"] for record in first if record["split"] == "train"}.isdisjoint( - {record["patient_id"] for record in first if record["split"] == "test"}) - - -def test_one_patient_cannot_form_train_and_test(): - with pytest.raises(ValidationError, match="At least two patient groups"): - split_patients([{"collection_id": "c4kc_kits", "patient_id": "A"}]) From 34157218a1803fb4025fc98ba1319018d6ccbdb1 Mon Sep 17 00:00:00 2001 From: Chris Bridge Date: Thu, 22 Oct 2026 09:36:25 -0400 Subject: [PATCH 5/5] Update README, tweaks to explanations in notebook --- README.md | 2 +- notebooks/advanced_topics/ct_seg_to_tensors.ipynb | 8 +++----- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/README.md b/README.md index b79ab68..b4e352f 100644 --- a/README.md +++ b/README.md @@ -50,7 +50,7 @@ Zenodo permanent citeable record: https://zenodo.org/records/22691908. Notebooks in this folder focus on topics that will require understanding of the basics, and aim to address more narrow use cases of IDC usage. * [Searching DICOM private tags](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/dicom_private_tags_intro.ipynb): all of DICOM attributes for the imaging data in IDC are searchable using BigQuery. DICOM private tags often contain critical information, such as diffusion b-values, but are a bit more tricky to access from BigQuery. In this tutorial you will learn how to accomplish this. -* [Validated CT/DICOM SEG tensors and patient-disjoint splits](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/ct_seg_to_tensors.ipynb): select referenced regular-grid C4KC-KiTS pairs with seg_index and volume_geometry_index, load HU volumes and separate mask channels with highdicom, validate source correspondence, and create a reproducible patient-disjoint split for PyTorch. +* [Loading paired CT/DICOM SEG objects as Pytorch tensors](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/ct_seg_to_tensors.ipynb): select referenced regular-grid C4KC-KiTS pairs with idc-index's seg_index, load HU volumes and associated segmentation mask with highdicom, and convert to matched PyTorch tensors. * [Working with IDC clinical data](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/clinical_data_intro.ipynb): this tutorial will help you become familiar with how clinical and other tabular data accompanying images in IDC is organized, and how you can use it to build cohorts. * [Direct access to tags/frames from GCS/AWS buckets](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/gcs_aws_direct_access.ipynb): learn how to access individual frames or tags of large DICOM files from the bucket without having to download the entire file (this notebook accompanies documentation article here: https://learn.canceridc.dev/data/downloading-data/direct-loading) * [Using DICOMweb to access IDC data](https://github.com/ImagingDataCommons/IDC-Tutorials/blob/master/notebooks/advanced_topics/idc_dicomweb_access.ipynb): both metadata and pixel data can be access using DICOMweb, which is particularly important while working with digital pathology, as it enables granular access to the individual pyramid tiles (frames) diff --git a/notebooks/advanced_topics/ct_seg_to_tensors.ipynb b/notebooks/advanced_topics/ct_seg_to_tensors.ipynb index fe4880a..7af6a27 100644 --- a/notebooks/advanced_topics/ct_seg_to_tensors.ipynb +++ b/notebooks/advanced_topics/ct_seg_to_tensors.ipynb @@ -196,10 +196,8 @@ "## 3. Load the segmentation and align it with the CT\n", "\n", "A DICOM SEG is a single multi-frame file. that can contain several *segments* (one per structure).\n", - "There are two variants: in BINARY SEGs (like those in C4KC-KiTS), each frame stores a binary mask\n", - "for one segment on one slice, and for LABELMAP SEGs, segments are combined together in a single\n", - "frame for a given slice. Highdicom's `highdicom.seg.Segmentation` class works with both types and\n", - "gives consistent results regardless of how frames are stored, hiding this difference from you.\n", + "This example uses BINARY segmentations; DICOM also supports FRACTIONAL and LABELMAP representations.\n", + "In BINARY SEGs (like those in C4KC-KiTS), each frame stores a binary mask for one segment on one slice.\n", "Let's read the file with `highdicom.seg.segread()` and see which segments it contains." ] }, @@ -234,7 +232,7 @@ " binary channel per segment in a trailing axis: `(slices, rows, columns, segments)`. (Use `combine_segments=True`\n", " instead if you want a single label map in which each voxel holds a segment number.)\n", "- `.match_geometry(ct_volume)` rearranges the mask volume so that it lies on exactly the same voxel grid as the\n", - " CT. This only pads, flips or permutes the array; it never interpolates, and it raises an error if the two\n", + " CT. This only pads, crops, flips, and/or permutes the array; it never interpolates, and it raises an error if the two\n", " volumes cannot be aligned.\n", "\n", "We pass `tol=1e-4` (in units of voxels) to `match_geometry()` because the SEG files in this collection store\n",