diff --git a/README.md b/README.md index b9b6fed..b4e352f 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. +* [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 new file mode 100644 index 0000000..7af6a27 --- /dev/null +++ b/notebooks/advanced_topics/ct_seg_to_tensors.ipynb @@ -0,0 +1,439 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "951549ba", + "metadata": {}, + "source": [ + "\"Open\n", + "\n", + "# Loading CT images and DICOM segmentations into PyTorch tensors\n", + "\n", + "Updated: October 2026\n", + "\n", + "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": "239b5820", + "metadata": {}, + "source": [ + "## Setup\n", + "\n", + "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": "091bbfb2", + "metadata": {}, + "outputs": [], + "source": [ + "%pip install -q \"idc-index>=0.13\" \"highdicom>=0.28\" torch matplotlib" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "51d509ad", + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "\n", + "import highdicom as hd\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import pydicom\n", + "import torch\n", + "from idc_index import IDCClient" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1243f56f-3f8a-45e2-af6b-94bdf1268f9c", + "metadata": {}, + "outputs": [], + "source": [ + "import idc_index\n", + "print(idc_index.__version__)\n" + ] + }, + { + "cell_type": "markdown", + "id": "21238b92", + "metadata": {}, + "source": [ + "## 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": "203205e1", + "metadata": {}, + "outputs": [], + "source": [ + "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": "8080c031", + "metadata": {}, + "source": [ + "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": "4ab233a5", + "metadata": {}, + "outputs": [], + "source": [ + "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": "56e625d5", + "metadata": {}, + "source": [ + "## 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": "b058c984", + "metadata": {}, + "outputs": [], + "source": [ + "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": "e409bb10", + "metadata": {}, + "source": [ + "## 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", + "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." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "57f4d9b2", + "metadata": {}, + "outputs": [], + "source": [ + "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": "72718f2c", + "metadata": {}, + "source": [ + "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, 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", + "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": "54bb6587", + "metadata": {}, + "outputs": [], + "source": [ + "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": "8337802a", + "metadata": {}, + "source": [ + "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": "305dff4c", + "metadata": {}, + "outputs": [], + "source": [ + "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": "2325769f", + "metadata": {}, + "source": [ + "## 4. Convert to PyTorch tensors\n", + "\n", + "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": "86f8c8b6", + "metadata": {}, + "outputs": [], + "source": [ + "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": "a97bbc63", + "metadata": {}, + "source": [ + "## 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", + "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": "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, 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.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": "31c213a4", + "metadata": {}, + "outputs": [], + "source": [ + "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": "1e4d45a8", + "metadata": {}, + "source": [ + "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", + "## Next steps\n", + "\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 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 the National Cancer Institute,\n", + "National Institutes of Health, under Task Order No. HHSN26110071 under Contract No. HHSN261201500003I.\n", + "\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", + "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 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}