sessions/tcia-idc/RSNA_2021_IDC_and_TCIA.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "RSNA 2021: TCIA and IDC",
      "private_outputs": true,
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/RSNA/AI-Deep-Learning-Lab-2021/blob/main/sessions/tcia-idc/RSNA_2021_IDC_and_TCIA.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gy9vfOYmSgiz"
      },
      "source": [
        "# RSNA 2021: Working with public datasets: TCIA and IDC\n",
        "\n",
        "**Short link to this notebook: https://tinyurl.com/RSNA21-IDC-TCIA**\n",
        "\n",
        "---\n",
        "\n",
        "_This notebook is part of the RSNA 2021 Deep Learning Lab tutorial series, see all notebooks here: https://github.com/RSNA/AI-Deep-Learning-Lab-2021_\n",
        "\n",
        "---\n",
        "\n",
        "_You will find more notebooks demonstrating how to work with the data in NCI Imaging Data Commons here: https://github.com/ImagingDataCommons/IDC-Examples/tree/master/notebooks_\n",
        "\n",
        "---\n",
        "\n",
        "The goal of this session is to introduce you to the two data repositories supported by the US National Cancer Institute:\n",
        "\n",
        "* [The Cancer Imaging Archive (TCIA)](https://www.cancerimagingarchive.net/)\n",
        "* [Imaging Data Commons (IDC)](https://imaging.datacommons.cancer.gov/), which is the imaging repository within [NCI Cancer Research Data Commons (CRDC)](https://datacommons.cancer.gov/)\n",
        "\n",
        "**Learning Objectives:**\n",
        "1. Understand basic capabilities of TCIA and IDC, and the differences between the two repositories.\n",
        "2. Explore relevant functionality of TCIA and IDC to support data exploration, cohort definition, and retrieval of the data.\n",
        "3. Learn how to analyze and visualize the data retrieved from TCIA/IDC on an example of segmentation of abdominal organs at risk.\n",
        "\n",
        "This notebook will guide you thought the complete process of identifying a relevant DICOM dataset from a public repository, retrieving it to the Colab VM, preparing it for processing by the specific analysis tool, installing the tool and applying it to the dataset, and visualizing the segmentation results produced by the tool.\n",
        "\n",
        "Note that it is not the purpose of this tutorial to promote a specific tool, or assess its robustness. We aim to provide an example of how a tool can be used for analyzing a sample dataset from TCIA/IDC. We hope that after completing this tutorial you will be empowered and motivated to experiment with more tools and apply them to more datasets in TCIA/IDC!\n",
        "\n",
        "**Session Authors**\n",
        "\n",
        "* Andrey Fedorov - Department of Radiology, Brigham and Women's Hospital, Boston, USA (contact: andrey.fedorov AT gmail dot com)\n",
        "* Justin Kirby - Frederick National Lab, Frederick, USA\n",
        "* Dennis Bontempi - Artificial Intelligence in Medicine Program, Brigham and Women's Hospital, Boston, USA\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fY1nlTbvZg2b"
      },
      "source": [
        "## Outline\n",
        "\n",
        "There is a growing number of tools being developed for analyzing medical images. More and more of those are shared openly by the authors to support dissemination of findings and reproducibility of scientific studies. \n",
        "\n",
        "However, getting such tools to work is not always straightforward. Some of the challenges include deployment of the tool, identifying datasets that are suitable for analysis using a specific tool, preprocessing of the data.\n",
        "\n",
        "In this notebook we will guide you through the process of deploying one such tool on a Colab VM, and demonstrate how you can utilize publicly available repositories of cancer imaging data to find relevant datasets, how to preprocess them for analysis by a specific tool, and how to visualize results of image segmentation produced by the tool. \n",
        "\n",
        "We will also demonstrate how existing segmentation results, available as either DICOM RTSTRUCT or DICOM SEG, can be loaded and visualized in the notebook. \n",
        "\n",
        "In this tutorial we will work with the nnU-Net segmentation tool developed by Isensee et al in the following publication:\n",
        "\n",
        "> Isensee, F., Jaeger, P. F., Kohl, S. A. A., Petersen, J. & Maier-Hein, K. H. nnU-Net: a self-configuring method for deep learning-based biomedical image segmentation. Nat. Methods 18, 203–211 (2021).\n",
        "\n",
        "and shared in https://github.com/MIC-DKFZ/nnUNet.\n",
        "\n",
        "Specifically, we will utilize the network pretrained to segment 4 abdominal organs at risk (AORs): heart, aorta, trachea and esophagus. The network was trained using the [SegTHOR (Segmentation of THoracic Organs at Risk) dataset](https://arxiv.org/abs/1912.05950) and is shared in this Zenodo entry (see [Task055_SegTHOR.zip](https://zenodo.org/record/4485926/files/Task055_SegTHOR.zip?download=1)):\n",
        "\n",
        "> Isensee, Fabian, Jäger, Paul F., Kohl, Simon A. A., Petersen, Jens, & Maier-Hein, Klaus H. (2021). pretrained models for 3D semantic image segmentation with nnU-Net (2.1). Zenodo. https://doi.org/10.5281/zenodo.4485926\n",
        "\n",
        "We warn you that for many this may not really be a beginner-level notebook! We hope this material will be helpful to start with the exploration of the DICOM data available in IDC, and experimentation with analyzing this data! \n",
        "\n",
        "Please contact us on the [Imaging Data Commons user forum](https://discourse.canceridc.dev/) if you have any questions related to this notebook.\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wMDzWEcX90lX"
      },
      "source": [
        "## Prerequisites\n",
        "\n",
        "* To use Colab, and to access data in IDC, you will need a [Google Account](https://support.google.com/accounts/answer/27441?hl=en)\n",
        "* Make sure your Colab instance has a GPU! For this check \"Runtime > Change runtime type\" and make sure to choose the GPU runtime.\n",
        "* To perform queries against IDC BigQuery tables you will need a cloud project. You can get started with Google Cloud free project with the following steps (they are also illustrated in [this short video](https://youtu.be/i08S0KJLnyw)):\n",
        "  1. Go to https://console.cloud.google.com/, and accept Terms and conditions.\n",
        "  2. Click \"Select a project\" button in the upper left corner of the screen, and then click \"New project\".\n",
        "  3. Open the console menu by clicking the ☰ menu icon in the upper left corner, and select \"Dashboard\". You will see information about your project, including your Project ID. Insert that project ID in the cell below in place of `REPLACE_ME_WITH_YOUR_PROJECT_ID`.\n",
        "\n",
        "  "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-dhXZF07n8Le"
      },
      "source": [
        "# initialize this variable with your Google Cloud Project ID!\n",
        "my_ProjectID = \"REPLACE_ME_WITH_YOUR_PROJECT_ID\"\n",
        "\n",
        "import os\n",
        "os.environ[\"GCP_PROJECT_ID\"] = my_ProjectID"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kGobegw-Bvz5"
      },
      "source": [
        "## nnU-Net model setup\n",
        "\n",
        "**Download of the pretrained network is the most time-consuming steps in this tutorial, this is why we start the download this early!**\n",
        "\n",
        "While the nnU-Net framework should take care of the model download (from Zenodo), some of the zip files containing the pre-trained weights are particularly large, so the download can take a lot of time, get stuck, or produce errors (as [reported by other users](https://github.com/MIC-DKFZ/nnUNet/issues/358#issue-726410474) and in the [repository FAQ](https://github.com/MIC-DKFZ/nnUNet/blob/master/documentation/common_problems_and_solutions.md#downloading-pretrained-models-unzip-cannot-find-zipfile-directory-in-one-of-homeisenseennunetdownload_16031094034174126)) .\n",
        "\n",
        "For this reason, and for the purpose of speeding up this tutorial, we decided to copy the relevant model weights in a shared Dropbox folder. In the following cells, we use Linux `wget` to pull such files from the folder - and use the nnU-Net framework command `nnUNet_install_pretrained_model_from_zip` to unpack and install the pre-trained model."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9sipNQgeB7fo"
      },
      "source": [
        "# create the directory tree\n",
        "!mkdir -p tutorial \n",
        "!mkdir -p tutorial/models tutorial/data tutorial/output\n",
        "!mkdir -p tutorial/data/dicom tutorial/data/processed tutorial/data/nnUNet_raw_data"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "AzPWnwbqB-Bt"
      },
      "source": [
        "# this will usually take between one and five minutes (but can sometimes take up to eight)\n",
        "seg_model_url = \"https://www.dropbox.com/s/m7es2ojn8h0ybhv/Task055_SegTHOR.zip?dl=0\"\n",
        "output_path = \"tutorial/models/Task055_SegTHOR.zip\"\n",
        "\n",
        "!wget -O $output_path $seg_model_url"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ySyKkGA4b0R6"
      },
      "source": [
        "## The Cancer Imaging Archive (TCIA)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FhYD_mvRdKTd"
      },
      "source": [
        "This course assumes you have some basic familiarity with The Cancer Imaging Archive.  If you have never used TCIA you can [watch this presentation from RSNA 2020](https://vimeo.com/595989800) in order to understand the mission of TCIA and services it provides to the research community.  Options for accessing data from TCIA are summarized at https://www.cancerimagingarchive.net/access-data/. The most relevant data access methods for this course are briefly summarized below."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PVyKi2oBfZxX"
      },
      "source": [
        "#### Browsing Collections & Analysis Results\n",
        "\n",
        "The most basic way to find data on TCIA is to [Browse Collections](https://www.cancerimagingarchive.net/collections) and [Browse Analysis Results](https://www.cancerimagingarchive.net/tcia-analysis-results/). Using the information in the table you can identify potential datasets of interest. Clicking on a given dataset takes you to a page which provides a description, data usage policy and citation guidelines, and links to download the data.  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pOwf2B0Pir0C"
      },
      "source": [
        "#### Different approaches for downloading TCIA data\n",
        "\n",
        "**NOTE: the below is an example of how to use the NBIA Data Retriever to download the data defined by TCIA manifest. You do not need to run this example for the purposes of this tutorial!**\n",
        "\n",
        "TCIA hosts a variety of image types and other related files, but the majority of its data are radiology images stored in DICOM format. When downloading DICOM images the download link will save a *.TCIA \"manifest\" file rather than the actual images. These manifest files must be opened with a helper application called the [NBIA Data Retriever](https://wiki.cancerimagingarchive.net/x/egOnAg). The Data Retriever can be installed on Windows, Mac and Linux operating systems. The Linux version also supports a command-line interface option which can be used on Google Colab.  Let's install that now using the following code:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Qsm1MeREmdYm"
      },
      "source": [
        "# install NBIA Data Retriever for downloading images \n",
        "# documentation available at https://wiki.cancerimagingarchive.net/display/NBIA/Downloading+TCIA+Images\n",
        "\n",
        "!mkdir /usr/share/desktop-directories/\n",
        "!wget -P /content/NBIA-Data-Retriever https://cbiit-download.nci.nih.gov/nbia/releases/ForTCIA/NBIADataRetriever_4.2/nbia-data-retriever-4.2.deb\n",
        "!dpkg -i /content/NBIA-Data-Retriever/nbia-data-retriever-4.2.deb\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tqKjpGdain6z"
      },
      "source": [
        "Now let's try it out by downloading a sample dataset from TCIA.  From the previously mentioned [Browse Collections](https://www.cancerimagingarchive.net/collections) page you can use the filter box (top right of the table of datasets) to filter out datasets of interest.  Try typing \"lung CT\" to reduce the table of datasets to those that we might want to analyze with nnU-Net.  This should reduce the table to 33 results.  \n",
        "\n",
        "In the interest of time, let's download a very small dataset such as [APOLLO-5-ESCA](https://doi.org/10.7937/n69a-7a26).  Clicking on the link to this dataset in the table will open its summary page. After reviewing the page to learn more about this dataset, scroll down to the bottom \"Data Access\" section.  Click the blue \"Download\" button for the Images to save the associated manifest file.  You can then upload this file to Colab, and open it using the NBIA Data Retriever.\n",
        "\n",
        "Alternatively, we can right click the blue \"Download\" button on the APOLLO-5-ESCA page and copy the URL of the manifest file to our clipboard.  Then we can use wget to download and save the manifest file to Colab.  After the manifest file is in place, we can download the image data using NBIA Data Retriever.  These steps are demostrated in the code below.\n",
        "\n",
        "**NOTE: You will need to accept the usage agreement terms before the download will start!**\n",
        "\n",
        "The cell below can take as long as 7 minutes - no need to run it now, but you can keep it handy when you want to use NBIA Data Retriever to download data from TCIA."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YJqE_XZCrHjB"
      },
      "source": [
        "# NBIA Data Retriever's Linux CLI documentation is at: https://wiki.cancerimagingarchive.net/display/NBIA/NBIA+Data+Retriever+Command+Line+Interface \n",
        "\n",
        "# wget the manifest file by copying the URL of the Download button instead of performing the manual download/upload\n",
        "!mkdir -p /content/TCIA_download_example/\n",
        "!wget -O /content/TCIA_download_example/apollo-5-esca.tcia https://wiki.cancerimagingarchive.net/download/attachments/96338104/APOLLO-5-ESCA-manifest.tcia?api=v2\n",
        "\n",
        "# execute the download of the manifest - uncomment the line below if you really want to download the data!\n",
        "# download may take about 7 minutes, and is not needed for the subsequent steps of the notebook\n",
        "#!/opt/nbia-data-retriever/nbia-data-retriever --cli /content/TCIA_downloads/apollo-5-esca.tcia -d /content/TCIA_download_example/"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FdJhtWIAcOj5"
      },
      "source": [
        "## Imaging Data Commons (IDC)\n",
        "\n",
        "The National Cancer Institute (NCI) [Cancer Research Data Commons (CRDC) ](https://datacommons.cancer.gov/) aims to establish a national cloud-based datascience infrastructure. [Imaging Data Commons (IDC)](https://imaging.datacommons.cancer.gov/) is a new component of CRDC supported by the Cancer Moonshot. The goal of IDC is to enable a broad spectrum of cancer researchers, with and without imaging expertise, to easily access and explore the value of deidentified imaging data and to support integrated analyses with nonimaging data. We achieve this goal by colocating versatile imaging collections with cloud-based computing resources and data exploration, visualization, and analysis tools. \n",
        "\n",
        "IDC provides access to curated imaging collections, accompanied by documentation, a user forum, and a growing number of analysis use cases that aim todemonstrate the value of a data commons framework applied to cancer imaging research.\n",
        "\n",
        "Key resources maintained by IDC are the following:\n",
        "* [public cancer imaging data](https://imaging.datacommons.cancer.gov/collections/) stored in Google Storage buckets and [public metadata tables](https://console.cloud.google.com/bigquery?p=bigquery-public-data&d=idc_current&page=dataset) that contain all of the DICOM metadata for IDC-hosted images\n",
        "  * these are now available in [Google Public Dataset Program](https://console.cloud.google.com/marketplace/product/gcp-public-data-idc/nci-idc-data), which supports free egress of IDC data out of the cloud\n",
        "* [radiology](https://viewer.imaging.datacommons.cancer.gov/viewer/1.3.6.1.4.1.14519.5.2.1.6279.6001.224985459390356936417021464571) and [pathology](https://viewer.imaging.datacommons.cancer.gov/slim/studies/1.3.6.1.4.1.5962.99.1.2463087261.2121647220.1625960757917.3.0/series/1.3.6.1.4.1.5962.99.1.2463087261.2121647220.1625960757917.2.0) zero-footprint viewers that can be used to visualize any of the data hosted by IDC in your browser\n",
        "* radiology and pathology use cases: reproducible analysis workflows that operate on IDC data, [available as Colab notebooks](https://github.com/ImagingDataCommons/IDC-Examples/tree/master/notebooks)\n",
        "* [user portal](https://imaging.datacommons.cancer.gov/) that can be used to explore the data available in IDC, visualize images and annotations, and build cohorts\n",
        "* [API](https://api.imaging.datacommons.cancer.gov/v1/swagger) that can be used for programmatic operations with IDC cohorts\n",
        "\n",
        "At the moment, most of the data you will find in IDC has been replicated from TCIA (exception to this are DICOM-converted digital pathology collections). In the future, IDC will host cancer imaging data from sources other than TCIA."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eRtmCas7dT6D"
      },
      "source": [
        "## Setup of the Colab VM\n",
        "\n",
        "\n",
        "\n",
        "In the following cells we will confirm you have a GPU before doing anything else, and will install and import all the Python dependencies. \n",
        "\n",
        "The main python packages we need to install are:\n",
        "* `nnunet` - which is the [codebase for the nn-UNet framework](https://github.com/MIC-DKFZ/nnUNet) we are going to be using for the segmentation step;\n",
        "* `pydicom`, a Python [package](https://github.com/pydicom/pydicom) that lets the use read, modify, and write DICOM data in an easy \"pythonic\" way - that we are going to use to distinguish different DICOM objects from each other."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zLvysANUArnm"
      },
      "source": [
        "### GPU checks"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Pf2j172mddvK"
      },
      "source": [
        "# check wether the Colab Instance was correctly initialized with a GPU instance\n",
        "gpu_list = !nvidia-smi --list-gpus\n",
        "\n",
        "has_gpu = False if \"failed\" in gpu_list[0] else True\n",
        "\n",
        "if not has_gpu:\n",
        "  print(\"Your Colab VM does not have a GPU - check \\\"Runtime > Change runtime type\\\"\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kL6xCvo3eKQ0"
      },
      "source": [
        "# check which model of GPU the notebook is equipped with - a Tesla K80 or T4\n",
        "# T4 is the best performing on the two - and can about half the GPU processing time\n",
        "\n",
        "!nvidia-smi"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZJhwBz2ABT_p"
      },
      "source": [
        "### Environment Setup\n",
        "\n",
        "Here we will configure the Linux environment variables needed to run the nnU-Net pipeline. \n",
        "\n",
        "Three main variables are needed by default to run the nnU-Net segmentation pipelines:\n",
        "* `nnUNet_raw_data_base` is the path to the folder where the segmentation pipeline expects to find the data to process;\n",
        "* `nnUNet_preprocessed` is the path to the folder where the preprocessed data are saved;\n",
        "* `RESULTS_FOLDER` is the path to the folder storing by default the model weights and, in our case, for simplicity, the segmentation masks produced by the pipeline.\n",
        "\n",
        "We will use the additional variable `PATH_TO_MODEL_FILE` to point to the location where the pre-trained model weights for the chosen model will be stored (more on this later).\n",
        "\n",
        "Please notice that these variables need to be set using `os.environ[]` in Google Colab - as `!export` is not sufficient to guarantee the variables are kept from one cell to the other. For more in-depth information regarding what the nnU-Net framework uses these folders for, please visit [the dedicated nnU-Net documentation page](https://github.com/MIC-DKFZ/nnUNet/blob/master/documentation/setting_up_paths.md)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gjiHm26cBdKm"
      },
      "source": [
        "# set env variables for the bash process\n",
        "import os\n",
        "os.environ['nnUNet_raw_data_base'] = \"/content/tutorial/data/nnUNet_raw_data/\"\n",
        "os.environ['nnUNet_preprocessed'] = \"/content/tutorial/data/processed/\"\n",
        "\n",
        "os.environ[\"RESULTS_FOLDER\"] = \"/content/tutorial/output/\"\n",
        "os.environ[\"PATH_TO_MODEL_FILE\"] = \"/content/tutorial/models/Task055_SegTHOR.zip\"\n",
        "\n",
        "dicom_sorted_dir = \"/content/tutorial/data/dicom\""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "40E3HnI5A0SX"
      },
      "source": [
        "### Install command-line tools\n",
        "\n",
        "\n",
        "[Plastimatch](https://plastimatch.org/index.html) is considered to be the swiss army knife of medical images processing: we will use it to convert DICOM (CT, RTSTRUCT) series to NRRD files - but it can be used for a multitude of other tasks, such as registration, resampling, cropping, and computing statistics to name a few. Plastimatch is also available as a 3DSlicer plug-in and can be used directly from the Slicer GUI.\n",
        "\n",
        "For the sake of clarity and simplicity, we will call Plastimatch from a very simple [Python wrapper](https://github.com/denbonte/pyplastimatch) written for the occasion (unfortunately, Plastimatch does not provide an official one)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZZL6-ByHA7XY"
      },
      "source": [
        "%%capture\n",
        "!sudo apt update\n",
        "\n",
        "!sudo apt install plastimatch"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LMovHnKgBEfC"
      },
      "source": [
        "!echo $(plastimatch --version)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "K_8zlhmdo0HD"
      },
      "source": [
        "[dcmqi](https://github.com/QIICR/dcmqi) is an open source library that can help with the conversion between imaging research formats and the standard DICOM representation for image analysis results. More specifically, you can use dcmqi convert DICOM Segmentation objects (DICOM SEG) into research formats, such as NIfTI and NRRD."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KAkmz4jApGh8"
      },
      "source": [
        "%%capture\n",
        "!wget https://github.com/QIICR/dcmqi/releases/download/v1.2.4/dcmqi-1.2.4-linux.tar.gz\n",
        "!tar zxvf dcmqi-1.2.4-linux.tar.gz\n",
        "!cp dcmqi-1.2.4-linux/bin/* /usr/local/bin/\n",
        "!segimage2itkimage --help"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fFRWendbBH72"
      },
      "source": [
        "Finally, we are going to install [Subversion](https://subversion.apache.org/), a tool that will allow us to clone GitHub repositories only partially (to save time and space)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YBEJRe-2BKah"
      },
      "source": [
        "%%capture\n",
        "\n",
        "!sudo apt install subversion"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WI1_BixWBN1O"
      },
      "source": [
        "!echo $(svn --version | head -n 2)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0Ev-JlpCAuMs"
      },
      "source": [
        "### Install Python packages"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7XMzq-8nRvKd"
      },
      "source": [
        "%%capture\n",
        "!pip install nnunet\n",
        "!pip install pydicom"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "62yhEkkjITzn"
      },
      "source": [
        "Unpack and install model we downloaded earlier (under `PATH_TO_MODEL_FILE`). This step can take about 1-2 minutes."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QOVsPuzxIOXX"
      },
      "source": [
        "%%capture\n",
        "!nnUNet_install_pretrained_model_from_zip $PATH_TO_MODEL_FILE"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BP1tahM9tzlm"
      },
      "source": [
        "Next we set up few things to help with visualization of the segmentations later."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cURjj8rzAa2L"
      },
      "source": [
        "import os\n",
        "import sys\n",
        "import shutil\n",
        "\n",
        "import time\n",
        "import gdown\n",
        "\n",
        "import json\n",
        "import pprint\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "\n",
        "import pydicom\n",
        "import nibabel as nib\n",
        "import SimpleITK as sitk\n",
        "\n",
        "from medpy.metric.binary import dc as dice_coef\n",
        "from medpy.metric.binary import hd as hausdorff_distance\n",
        "from medpy.metric.binary import asd as avg_surf_distance\n",
        "\n",
        "from medpy.filter.binary import largest_connected_component\n",
        "\n",
        "# use the \"tensorflow_version\" magic to make sure TF 1.x is imported\n",
        "%tensorflow_version 1.x\n",
        "import tensorflow as tf\n",
        "import keras\n",
        "\n",
        "print(\"Python version               : \", sys.version.split('\\n')[0])\n",
        "print(\"Numpy version                : \", np.__version__)\n",
        "print(\"TensorFlow version           : \", tf.__version__)\n",
        "print(\"Keras (stand-alone) version  : \", keras.__version__)\n",
        "\n",
        "print(\"\\nThis Colab instance is equipped with a GPU.\")\n",
        "\n",
        "# ----------------------------------------\n",
        "\n",
        "#everything that has to do with plotting goes here below\n",
        "import matplotlib\n",
        "matplotlib.use(\"agg\")\n",
        "\n",
        "import matplotlib.pyplot as plt\n",
        "from matplotlib.colors import ListedColormap\n",
        "\n",
        "%matplotlib inline\n",
        "%config InlineBackend.figure_format = \"png\"\n",
        "\n",
        "import ipywidgets as ipyw\n",
        "\n",
        "## ----------------------------------------\n",
        "\n",
        "# create new colormap appending the alpha channel to the selected one\n",
        "# (so that we don't get a \\\"color overlay\\\" when plotting the segmask superimposed to the CT)\n",
        "cmap = plt.cm.Reds\n",
        "my_reds = cmap(np.arange(cmap.N))\n",
        "my_reds[:,-1] = np.linspace(0, 1, cmap.N)\n",
        "my_reds = ListedColormap(my_reds)\n",
        "\n",
        "cmap = plt.cm.Greens\n",
        "my_greens = cmap(np.arange(cmap.N))\n",
        "my_greens[:,-1] = np.linspace(0, 1, cmap.N)\n",
        "my_greens = ListedColormap(my_greens)\n",
        "\n",
        "cmap = plt.cm.Blues\n",
        "my_blues = cmap(np.arange(cmap.N))\n",
        "my_blues[:,-1] = np.linspace(0, 1, cmap.N)\n",
        "my_blues = ListedColormap(my_blues)\n",
        "\n",
        "cmap = plt.cm.spring\n",
        "my_spring = cmap(np.arange(cmap.N))\n",
        "my_spring[:,-1] = np.linspace(0, 1, cmap.N)\n",
        "my_spring = ListedColormap(my_spring)\n",
        "## ----------------------------------------\n",
        "\n",
        "import seaborn as sns"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UrAh2kjqb2Be"
      },
      "source": [
        "# PyPlastimatch - python wrapper for Plastimatch (and interactive notebook visualisation)\n",
        "!svn checkout https://github.com/AIM-Harvard/pyplastimatch/trunk/pyplastimatch tutorial/pyplastimatch"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GsePTNBmM9sw"
      },
      "source": [
        "# dicomsort is the pythong package that can sort DICOM files into\n",
        "# folder organization based on user-specified DICOM attributes\n",
        "!git clone https://github.com/pieper/dicomsort.git"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Duil1tviGWaG"
      },
      "source": [
        "## Data selection\n",
        "\n",
        "The model is trained to segment organs in chest CT. In the following cell we initialize variable that points to the specific CT image (more precisely, specific DICOM CT image series) using DICOM `SeriesInstanceUID` attribute that we will use with the segmentation tool. Utilizing this unique identifier you can retrieve the imaging series from either TCIA or IDC.\n",
        "\n",
        "Once you are done with this example, we will give you instructions how to find more chest CT image series that you can use to experiment with the nnU-Net segmentation tool!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vOBcXkDoHph1"
      },
      "source": [
        "my_SeriesInstanceUID = \"1.3.6.1.4.1.32722.99.99.232988001551799080335895423941323261228\""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "azGjFr-qG0us"
      },
      "source": [
        "## Visualization and download of data from IDC\n",
        "\n",
        "In order to work with Google Cloud, you will need to have a GCP project. You should have done this in the prerequisites. In the following cell you will be authenticated with your Google account. Make sure you use the same Google account you used in the Prerequisites section to create your Google Cloud project!\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZDHCwK-TG4n2"
      },
      "source": [
        "from google.colab import auth\n",
        "auth.authenticate_user()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "naJ5OAKZnZpG"
      },
      "source": [
        "In order to use data hosted by IDC effectively, you will need to utilize metadata to navigate what data is available and to select specific files that are relevant in your analysis. The main metadata table you will need for this purpose is the [`bigquery-public-data.idc_current.dicom_all`](https://console.cloud.google.com/bigquery?p=bigquery-public-data&d=idc_current&t=dicom_all&page=table) table.\n",
        "\n",
        "This query has one row per file hosted by IDC. All of IDC data is in DICOM format, and each of the rows in this table will have all of the DICOM attributes extracted from a given file. It will also have various columns containing non-DICOM metadata, such as the name of the collection where the file is included, size of the file, and URL that can be used to retrieve that file.\n",
        "\n",
        "To query IDC BigQuery tables, you can use one of the following approaches:\n",
        "1. `%%bigquery` magic will allow you to define your query in plain SQL, and load the result of the query into a Pandas dataframe.\n",
        "2. [BigQuery Python API](https://googleapis.dev/python/bigquery/latest/index.html) is more flexible in allowing you to parameterize your query.\n",
        "3. [Google Cloud BigQuery console](https://console.cloud.google.com/bigquery) is very convenient for interactive query exploration of tables.\n",
        "4. [`gcloud bq`](https://cloud.google.com/bigquery/docs/bq-command-line-tool) is the command line tool that comes as part of [Cloud SDK](https://cloud.google.com/sdk) and is convenient for scripting interactions from the shell. Cloud SDK is preinstalled on Colab.\n",
        "\n",
        "In the following cells we will utilize `%%bigquery`, Python BigQuery SDK and BigQuery console for working with IDC BigQuery tables.\n",
        "\n",
        "First, to verify that you are authenticated, and your project ID is working, let's run a test query against IDC BigQuery table to get the summary statistics about the  of data available in IDC.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "C2tx4sBInsFd"
      },
      "source": [
        "%%bigquery --project=$my_ProjectID idc_summary_df\n",
        "\n",
        "SELECT\n",
        "  COUNT(DISTINCT(PatientID)) AS total_patients,\n",
        "  COUNT(DISTINCT(StudyInstanceUID)) AS total_studies,\n",
        "  COUNT(DISTINCT(SeriesInstanceUID)) AS total_series,\n",
        "  COUNT(DISTINCT(SOPInstanceUID)) AS total_instances,\n",
        "  COUNT(DISTINCT(collection_id)) AS total_collections,\n",
        "  SUM(instance_size)/1099511627776 AS total_size_TB\n",
        "FROM\n",
        "  `bigquery-public-data.idc_current.dicom_all`"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "B8n3NTNvy_aQ"
      },
      "source": [
        "idc_summary_df"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "t9YKMGVvVZX_"
      },
      "source": [
        "Given `SeriesInstanceUID` value identifying the image series, we can query the IDC metadata table to get the list of files (defined by the Google Storage URLs) corresponding to this series.\n",
        "\n",
        "All of the DICOM metadata for each of the DICOM files is available in the BigQuery table we will be querying. We will get not just the `gcs_url`, but also identifiers for the Study, Series and Instance, to better understand organization of data, and since `StudyInstanceUID` will be handy later when we get to the visualization of the data."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "joRbZ81GWNz1"
      },
      "source": [
        "from google.cloud import bigquery\n",
        "bq_client = bigquery.Client(my_ProjectID)\n",
        "\n",
        "selection_query = f\"\\\n",
        "  SELECT  \\\n",
        "    StudyInstanceUID, \\\n",
        "    SeriesInstanceUID, \\\n",
        "    SOPInstanceUID, \\\n",
        "    gcs_url \\\n",
        "  FROM \\\n",
        "    `bigquery-public-data.idc_current.dicom_all` \\\n",
        "  WHERE \\\n",
        "    SeriesInstanceUID = \\\"{my_SeriesInstanceUID}\\\"\"\n",
        "\n",
        "selection_result = bq_client.query(selection_query)\n",
        "selection_df = selection_result.result().to_dataframe()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KgY51LXoYWzO"
      },
      "source": [
        "Let's look at the resulting table. Each row corresponds to a single DICOM file, which can be downloaded using the `gcs_url` URL. The values of `StudyInstanceUID` and `SeriesInstanceUID` are identical for all files, since they belong to the same study and series, but `SOPInstanceUID` values uniquely identify the specific DICOM file (instance)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tWf6sYpOYZ4S"
      },
      "source": [
        "import pandas as pd\n",
        "pd.set_option('display.max_colwidth', None)\n",
        "\n",
        "selection_df"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MFgSARVk9bUA"
      },
      "source": [
        "We will initiallize `StudyInstanceUID` as we will need it for several steps in the following cells."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UALQQLrT9mTq"
      },
      "source": [
        "import pandas as pd\n",
        "my_StudyInstanceUIDs = selection_df['StudyInstanceUID'].unique()\n",
        "my_StudyInstanceUIDs.sort()\n",
        "my_StudyInstanceUID = my_StudyInstanceUIDs[0]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mG2USrvUd7yl"
      },
      "source": [
        "We can use the IDC radiology image viewer to visualize the series and study of interest. This can be done using the identifiers we have in the dataframe we obtained in the earlier query."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HawHLzymeGcb"
      },
      "source": [
        "def get_idc_viewer_url(studyUID, seriesUID=None):\n",
        "  url = \"https://viewer.imaging.datacommons.cancer.gov/viewer/\"+studyUID\n",
        "  if seriesUID is not None:\n",
        "    url = url+\"?seriesInstanceUID=\"+seriesUID\n",
        "  return url\n",
        "\n",
        "print(\"URL to view the entire study:\")\n",
        "print(get_idc_viewer_url(my_StudyInstanceUID))\n",
        "print()\n",
        "print(\"URL to view the specific series:\")\n",
        "print(get_idc_viewer_url(my_StudyInstanceUID, my_SeriesInstanceUID))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NRfoM4AXcmCX"
      },
      "source": [
        "Now that we confirmed the series we identified is indeed a CT of the chest, let's download the DICOM files for this series to the Colab VM. We will use the `gsutil` command line tool to fetch each of the files defined by the list of GCS URLs saved in a plain text file.\n",
        "\n",
        "If you want to download large number of files from IDC, make sure to check out [this documentation article](https://learn.canceridc.dev/data/downloading-data) to learn about performance optimizations of the download!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "X4YXrk6PIkJe"
      },
      "source": [
        "# save the list of GCS URLs into a file\n",
        "import os\n",
        "idc_download_folder = \"/content/IDC_downloads\"\n",
        "if not os.path.exists(idc_download_folder):\n",
        "  os.mkdir(idc_download_folder)\n",
        "\n",
        "selection_manifest = os.path.join(idc_download_folder, \"idc_manifest.txt\")\n",
        "selection_df[\"gcs_url\"].to_csv(selection_manifest, header=False, index=False)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mIKSKvW8mdsN"
      },
      "source": [
        "# confirm the resulting manifest has as many lines as the number of rows in the\n",
        "# dataframe we initialized earlier\n",
        "!cat /content/IDC_downloads/idc_manifest.txt |wc"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "KdZSVftdjFPl"
      },
      "source": [
        "# let's make sure the download folder is clean, in case you ran this cell earlier\n",
        "# for a different dataset\n",
        "!rm -rf /content/IDC_downloads/*.dcm"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6ydwu1tSYScJ"
      },
      "source": [
        "# download is this simple!\n",
        "%%capture\n",
        "\n",
        "!cat /content/IDC_downloads/idc_manifest.txt | gsutil -m cp -I /content/IDC_downloads"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5WqUQnULEbPN"
      },
      "source": [
        "\n",
        "Now you have the data ready for the next steps of the processing. If you are interested how to download the same series using TCIA API, continue to the next section. Otherwise you can skip to **Sort the DICOM files**."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hxJzr-GDIlAV"
      },
      "source": [
        "## Download of data from TCIA\n",
        "\n",
        "To demonstrate how data can be retrieved from TCIA, we will utilize the [NBIA REST API](https://wiki.cancerimagingarchive.net/x/fILTB) instead of the NBIA Data Retriever to download the data. The reason for using NBIA REST API is because we do not have a manifest file to pass to the Data Retriever, and also to demonstrate the alternative approach to TCIA data download.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DOYNkKApoNFr"
      },
      "source": [
        "!mkdir -p /content/TCIA_downloads && rm -rf /content/TCIA_downloads/*"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fASAVcz2GG3j"
      },
      "source": [
        "# download zip file with the series instances\n",
        "import requests, os, zipfile\n",
        "\n",
        "tcia_download_folder = \"/content/TCIA_downloads\"\n",
        "params = {\"SeriesInstanceUID\":my_SeriesInstanceUID}\n",
        "image_request = requests.get(\" https://services.cancerimagingarchive.net/nbia-api/services/v1/getImage\", params=params, stream=True)\n",
        "print(f\"Completed request: {image_request.url}\")\n",
        "if image_request.status_code == 200:\n",
        "  series_zip_name = os.path.join(tcia_download_folder, f\"{my_SeriesInstanceUID}.zip\")\n",
        "  with open(series_zip_name, \"wb\") as f:\n",
        "    for chunk in image_request.iter_content(chunk_size=1024):\n",
        "      f.write(chunk)\n",
        "\n",
        "  # extract individual instances from the series zip file\n",
        "  series_folder_name = os.path.join(tcia_download_folder, my_SeriesInstanceUID)\n",
        "  if not os.path.exists(series_folder_name):\n",
        "    os.mkdir(series_folder_name)\n",
        "  with zipfile.ZipFile(series_zip_name, 'r') as zip_ref:\n",
        "    zip_ref.extractall(series_folder_name)\n",
        "else:\n",
        "  print(f\"Failed with {r.status_code}\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZvlSBfPRdRpz"
      },
      "source": [
        "## Sort the DICOM files\n",
        "\n",
        "Let's check the content downloaded from the two repositories is identical, just in case. This section will apply only if you downloaded the data both from IDC and TCIA."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "q9aVI7qAbaT0"
      },
      "source": [
        "!mkdir -p IDC_sorted\n",
        "!python dicomsort/dicomsort.py -k -u IDC_downloads IDC_sorted/%PatientID/%StudyInstanceUID/%SeriesInstanceUID/%SOPInstanceUID.dcm"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ewcZsa1rzKHi"
      },
      "source": [
        "Move the sorted data into the right place"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "DZ_LAvUGEGxj"
      },
      "source": [
        "!mkdir -p tutorial/data/dicom && rm -rf tutorial/data/dicom/* && mv IDC_sorted/* tutorial/data/dicom"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5jerJXiNbcgo"
      },
      "source": [
        "## Data Pre-processing\n",
        "\n",
        "In order to run the AI segmentation pipeline, we need to convert the DICOM data in a format required by nnU-Net.\n",
        "\n",
        "Using the simple Plastimatch wrapper, let's convert the DICOM CT series in both NRRD (very flexible, simple handling with SimpleITK) and NIfTI (as required by the nnU-Net pipeline) format.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "vI7JuZhZbiXP"
      },
      "source": [
        "from tutorial.pyplastimatch import pyplastimatch as pypla\n",
        "from tutorial.pyplastimatch.utils import viz as viz_utils\n",
        "from tutorial.pyplastimatch.utils import data as data_utils"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7gLUnOsJDmaU"
      },
      "source": [
        "pat = os.listdir(dicom_sorted_dir)[0]\n",
        "\n",
        "# study_uid was initialized earlier, when we decided which study to download\n",
        "\n",
        "# directories storing NRRD and NIfTI files\n",
        "base_preproc_path = \"tutorial/data/processed\"\n",
        "\n",
        "pat_dir_path_nrrd = os.path.join(base_preproc_path, \"nrrd\", pat)\n",
        "pat_dir_path_nii = os.path.join(base_preproc_path, \"nii\", pat)\n",
        "  \n",
        "# patient subfolder where all the preprocessed NRRDs will be stored\n",
        "if not os.path.exists(pat_dir_path_nrrd): os.makedirs(pat_dir_path_nrrd)\n",
        "  \n",
        "# patient subfolder where all the preprocessed NIfTIs will be stored\n",
        "if not os.path.exists(pat_dir_path_nii): os.makedirs(pat_dir_path_nii)\n",
        "\n",
        "# path to the directory where the DICOM CT file is stored\n",
        "path_to_ct_dir = os.path.join(\"tutorial/data/dicom\", pat,\n",
        "                              my_StudyInstanceUID, my_SeriesInstanceUID)\n",
        "\n",
        "# path to the files where the NRRD and NIfTI CTs will be stored\n",
        "ct_nrrd_path = os.path.join(pat_dir_path_nrrd, pat + \"_ct.nrrd\")\n",
        "ct_nii_path = os.path.join(pat_dir_path_nii, pat + \"_ct.nii.gz\")\n",
        "\n",
        "verbose = True\n",
        "\n",
        "# logfile for the plastimatch conversion\n",
        "log_file_path_nrrd = os.path.join(pat_dir_path_nrrd, pat + '_pypla.log')\n",
        "log_file_path_nii = os.path.join(pat_dir_path_nii, pat + '_pypla.log')\n",
        "  \n",
        "# DICOM CT to NRRD conversion (if the file doesn't exist yet)\n",
        "if not os.path.exists(ct_nrrd_path):\n",
        "  convert_args_ct = {\"input\" : path_to_ct_dir,\n",
        "                     \"output-img\" : ct_nrrd_path}\n",
        "  \n",
        "  # clean old log file if it exist\n",
        "  if os.path.exists(log_file_path_nrrd): os.remove(log_file_path_nrrd)\n",
        "  \n",
        "  pypla.convert(verbose = verbose, path_to_log_file = log_file_path_nrrd, **convert_args_ct)\n",
        "\n",
        "# DICOM CT to NIfTI conversion (if the file doesn't exist yet)\n",
        "if not os.path.exists(ct_nii_path):\n",
        "  convert_args_nii = {\"input\" : path_to_ct_dir, \n",
        "                      \"output-img\" : ct_nii_path}\n",
        "  \n",
        "  # clean old log file if it exist\n",
        "  if os.path.exists(log_file_path_nii): os.remove(log_file_path_nii)\n",
        "  \n",
        "  pypla.convert(verbose = verbose, path_to_log_file = log_file_path_nii, **convert_args_nii)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CcYxd4NAFNps"
      },
      "source": [
        "As the final step before running the segmentation pipeline, we need to make sure the folder storing the data follows the structure required by the nnU-Net framework, described at the [dedicated documentation page](https://github.com/MIC-DKFZ/nnUNet/blob/master/documentation/dataset_conversion.md)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3CFhNcscFOtA"
      },
      "source": [
        "# create a folder (random task name) for nnU-Net inference\n",
        "proc_folder_path = os.path.join(os.environ[\"nnUNet_raw_data_base\"],\n",
        "                                \"segthor\", \"imagesTs\")\n",
        "\n",
        "!mkdir -p $proc_folder_path\n",
        "\n",
        "# populate the folder following the nnU-Net naming conventions\n",
        "copy_path = os.path.join(proc_folder_path, pat + \"_0000.nii.gz\")\n",
        "\n",
        "# copy NIfTI to the right dir for nnU-Net processing\n",
        "if not os.path.exists(copy_path):\n",
        "  shutil.copy(src = ct_nii_path, dst = copy_path)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "O2ci1vFuJTlG"
      },
      "source": [
        "## Segmentation of thoracic structures from CT series"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xvazHm6AFXVx"
      },
      "source": [
        "### Inference \n",
        "\n",
        "In order to run the segmentation pipeline, we can follow the [\"run inference\" section of the nnU-Net documentation](https://github.com/MIC-DKFZ/nnUNet#how-to-run-inference-with-pretrained-models), specifying the path to the input and output folders defined in the sections above, and the pretrained model we want to use (i.e., the one we downloaded earlier).\n",
        "\n",
        "For the purpose of this notebook, to make the processing faster, we are not going to use an ensemble of different U-Net configurations for inference or test time augmentation (TTA). You are invited to explore these options later - and if you decide to do so, you can read [this example](https://github.com/MIC-DKFZ/nnUNet/blob/master/documentation/inference_example_Prostate.md) from the nnU-Net documentation to learn how this can be achieved.\n",
        "\n",
        "To learn more about all the arguments that can be specified to the `nnUNet_predict` command, run `nnUNet_predict --help`.\n",
        "\n",
        "The following step will take a couple of minutes."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "cAEBKIbQFgLw"
      },
      "source": [
        "# run the inference phase\n",
        "# accepted options for --model are: 2d, 3d_lowres, 3d_fullres or 3d_cascade_fullres\n",
        "!nnUNet_predict --input_folder \"tutorial/data/nnUNet_raw_data/segthor/imagesTs\" \\\n",
        "                --output_folder $RESULTS_FOLDER \\\n",
        "                --task_name \"Task055_SegTHOR\" --model 2d --disable_tta "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wu7UIDM6JY2Z"
      },
      "source": [
        "### Post-processing of inference results\n",
        "\n",
        "After the inference is finished, we can convert the segmentation masks back to NRRD for visualisation purposes and for easier handling."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Lq1fP6XHJfoS"
      },
      "source": [
        "pred_nii_path = os.path.join(os.environ[\"RESULTS_FOLDER\"], pat + \".nii.gz\")\n",
        "\n",
        "sitk_ct = sitk.ReadImage(ct_nrrd_path)\n",
        "\n",
        "nrrd_spacing = sitk_ct.GetSpacing()\n",
        "nrrd_dim = sitk_ct.GetSize()\n",
        "\n",
        "nii_spacing = tuple(nib.load(pred_nii_path).header['pixdim'][1:4])\n",
        "nii_dim = tuple(nib.load(pred_nii_path).get_fdata().shape)\n",
        "\n",
        "assert (nrrd_spacing == nii_spacing) & (nrrd_dim == nii_dim)\n",
        "\n",
        "## ----------------------------------------\n",
        "# NIfTI TO NRRD CONVERSION\n",
        "\n",
        "# path to the output NRRD file (inferred segmasks)\n",
        "pred_nrrd_path = os.path.join(pat_dir_path_nrrd, pat + \"_pred_segthor.nrrd\")\n",
        "log_file_path = os.path.join(pat_dir_path_nrrd, pat + \"_pypla.log\")\n",
        "\n",
        "# Inferred NIfTI segmask to NRRD\n",
        "convert_args_pred = {\"input\" : pred_nii_path, \n",
        "                     \"output-img\" : pred_nrrd_path}\n",
        "\n",
        "pypla.convert(path_to_log_file = log_file_path, **convert_args_pred)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ShIjT6jEJkK-"
      },
      "source": [
        "### Visualising segmentation results\n",
        "\n",
        "We can visualise the raw AI-inferred segmentation mask (heart, aorta, esophagus, amd treachea - in green, yellow, red, and blue, respectively) and compare the heart (and esophagus, if available for the randomly selected patient) segmentation to the manual delineation."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "p__wFByuJjrL"
      },
      "source": [
        "# load NRRD volumes\n",
        "ct_nrrd = sitk.GetArrayFromImage(sitk_ct)\n",
        "\n",
        "# inferred segmask\n",
        "pred_nrrd_segthor = sitk.GetArrayFromImage(sitk.ReadImage(pred_nrrd_path))\n",
        "\n",
        "pred_nrrd_esophagus = np.copy(pred_nrrd_segthor)\n",
        "pred_nrrd_heart = np.copy(pred_nrrd_segthor)\n",
        "pred_nrrd_trachea = np.copy(pred_nrrd_segthor)\n",
        "pred_nrrd_aorta = np.copy(pred_nrrd_segthor)\n",
        "  \n",
        "# zero every segmask other than the esophagus and make the mask binary (0/1)\n",
        "pred_nrrd_esophagus[pred_nrrd_segthor != 1] = 0\n",
        "pred_nrrd_esophagus[pred_nrrd_esophagus != 0] = 1\n",
        "  \n",
        "# zero every segmask other than the heart and make the mask binary (0/1)\n",
        "pred_nrrd_heart[pred_nrrd_segthor != 2] = 0\n",
        "pred_nrrd_heart[pred_nrrd_heart != 0] = 1\n",
        "  \n",
        "# zero every segmask other than the trachea and make the mask binary (0/1)\n",
        "pred_nrrd_trachea[pred_nrrd_segthor != 3] = 0\n",
        "pred_nrrd_trachea[pred_nrrd_trachea != 0] = 1\n",
        "  \n",
        "# zero every segmask other than the aorta and make the mask binary (0/1)\n",
        "pred_nrrd_aorta[pred_nrrd_segthor != 4] = 0\n",
        "pred_nrrd_aorta[pred_nrrd_aorta != 0] = 1\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_TiIfDIvQ317"
      },
      "source": [
        "_ = viz_utils.AxialSliceSegmaskViz(ct_volume = ct_nrrd,\n",
        "                                          segmask_dict = {\"Heart\" : pred_nrrd_heart,\n",
        "                                                             \"Aorta\" : pred_nrrd_aorta,\n",
        "                                                             \"Trachea\" : pred_nrrd_trachea,\n",
        "                                                             \"Esophagus\" : pred_nrrd_esophagus},\n",
        "                                          segmask_cmap_dict = {\"Heart\" : my_greens,\n",
        "                                                               \"Aorta\" : my_spring,\n",
        "                                                               \"Esophagus\" : my_reds,\n",
        "                                                               \"Trachea\" : my_blues},\n",
        "                                          dpi = 100)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "063LqMlMiRdl"
      },
      "source": [
        "## Comparison of segmentations with existing annotations\n",
        "\n",
        "**NOTE: The cells in this section are applicable for the specific study/series we use as the example (`SeriesInstanceUID = 1.3.6.1.4.1.32722.99.99.203715003805996641695765332389135385095`). If you experiment with other series, this will either not be applicable, or you will need to adjust the code.**\n",
        "\n",
        "Many of the collections contain annotations alongside images, which can be used to train new algorithms, or to evaluate performance of algorithms. We can view the entire imaging study, as we did earlier, to see what annotations are available for the series we segmented. As you can see, this study contains segmentations of organs saved in RTSTRUCT and SEG series."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Wvu-rxiQnDxp"
      },
      "source": [
        "print(get_idc_viewer_url(my_StudyInstanceUID))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MzXiHlbunkYh"
      },
      "source": [
        "### Conversion and visualization of RTSTRUCT annotations\n",
        "\n",
        "We first get the names of the structures in the RTSTRUCT series that are available in the DICOM study containing our CT series using the following query:\n",
        "\n",
        "```sql\n",
        "SELECT\n",
        "  SeriesInstanceUID,\n",
        "  structureSetROISequence.ROIName AS ROIName\n",
        "FROM\n",
        " `bigquery-public-data.idc_current.dicom_all`\n",
        "CROSS JOIN\n",
        " UNNEST (StructureSetROISequence) AS structureSetROISequence\n",
        "WHERE\n",
        " Modality = \"RTSTRUCT\" AND StudyInstanceUID = <my_StudyInstanceUID>\n",
        " ```"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CPSF-l5siWKk"
      },
      "source": [
        "rt_selection_query = f\"\\\n",
        "  SELECT \\\n",
        "    SeriesInstanceUID,\\\n",
        "    structureSetROISequence.ROIName AS ROIName, \\\n",
        "    gcs_url \\\n",
        "  FROM \\\n",
        "    `bigquery-public-data.idc_current.dicom_all` \\\n",
        "  CROSS JOIN \\\n",
        "    UNNEST (StructureSetROISequence) AS structureSetROISequence \\\n",
        "  WHERE \\\n",
        "    Modality = \\\"RTSTRUCT\\\" AND StudyInstanceUID = \\\"{my_StudyInstanceUID}\\\"\"\n",
        "\n",
        "rt_selection_result = bq_client.query(rt_selection_query)\n",
        "rt_selection_df = rt_selection_result.result().to_dataframe()\n",
        "\n",
        "rt_selection_df"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oYG-4UAXq4Eb"
      },
      "source": [
        "As you can see, in this specific case, there is a single RTSTRUCT series that has segmentations of the structures that we also segmented in the previous step using nnU-Net. Let's download the corresponding DICOM file, and convert RTSTRUCT into a representation that we can use for visualization."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_NhspiuVspVD"
      },
      "source": [
        "!mkdir -p /content/IDC_downloads/RTSTRUCT\n",
        "!rm -rf /content/IDC_downloads/RTSTRUCT/*.dcm"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZBAlYZY5scjR"
      },
      "source": [
        "path_to_rt_dir = \"/content/IDC_downloads/RTSTRUCT\"\n",
        "rt_selection_manifest = \"/content/IDC_downloads/rt_manifest.txt\"\n",
        "rt_selection_df[\"gcs_url\"].to_csv(rt_selection_manifest, header=False, index=False)\n",
        "\n",
        "!cat /content/IDC_downloads/rt_manifest.txt | gsutil -m cp -I /content/IDC_downloads/RTSTRUCT"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Da_0aVnyiW3C"
      },
      "source": [
        "# path to the files where the NRRD RTSTRUCTs will be stored\n",
        "rt_folder = os.path.join(pat_dir_path_nrrd, \"RTSTRUCT\")\n",
        "rt_list_path = os.path.join(pat_dir_path_nrrd, \"RTSTRUCT_content\")\n",
        "\n",
        "# DICOM RTSTRUCT to NRRD conversion (if the file doesn't exist yet)\n",
        "if not os.path.exists(rt_folder):\n",
        "  convert_args_rt = {\"input\" : path_to_rt_dir, \n",
        "                     \"referenced-ct\" : path_to_ct_dir,\n",
        "                     \"output-prefix\" : rt_folder,\n",
        "                     \"prefix-format\" : 'nrrd',\n",
        "                     \"output-ss-list\" : rt_list_path}\n",
        "  \n",
        "  # clean old log file if it exist\n",
        "  if os.path.exists(log_file_path_nrrd): os.remove(log_file_path_nrrd)\n",
        "  \n",
        "  pypla.convert(verbose = verbose, path_to_log_file = log_file_path_nrrd, **convert_args_rt)\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "02Kn4GcmiVro"
      },
      "source": [
        "# manual segmask (from the RTSTRUCT)\n",
        "rt_segmask_heart = os.path.join(pat_dir_path_nrrd, \"RTSTRUCT\", \"Heart.nrrd\")\n",
        "rt_nrrd_heart = sitk.GetArrayFromImage(sitk.ReadImage(rt_segmask_heart))\n",
        "\n",
        "try:\n",
        "  rt_segmask_esophagus = os.path.join(pat_dir_path_nrrd, \"RTSTRUCT\", \"Esophagus.nrrd\")\n",
        "  rt_nrrd_esophagus = sitk.GetArrayFromImage(sitk.ReadImage(rt_segmask_esophagus))\n",
        "except:\n",
        "  # for the sake of simplicity, fill the volume with zeros\n",
        "  # (so that we can keep the code that comes after the same)\n",
        "  rt_nrrd_esophagus = np.zeros(rt_nrrd_heart.shape)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "I2Ll9UABJsFg"
      },
      "source": [
        "_ = viz_utils.AxialSliceSegmaskComparison(ct_volume = ct_nrrd,\n",
        "                                          segmask_ai_dict = {\"Heart\" : pred_nrrd_heart,\n",
        "                                                             \"Aorta\" : pred_nrrd_aorta,\n",
        "                                                             \"Trachea\" : pred_nrrd_trachea,\n",
        "                                                             \"Esophagus\" : pred_nrrd_esophagus},\n",
        "                                          segmask_manual_dict = {\"Heart\" : rt_nrrd_heart,\n",
        "                                                                 \"Esophagus\" : rt_nrrd_esophagus},\n",
        "                                          segmask_cmap_dict = {\"Heart\" : my_greens,\n",
        "                                                               \"Aorta\" : my_spring,\n",
        "                                                               \"Esophagus\" : my_reds,\n",
        "                                                               \"Trachea\" : my_blues},\n",
        "                                          dpi = 100)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BmYCHg1EqWJS"
      },
      "source": [
        "### Conversion and visualization of SEG annotations\n",
        "\n",
        "DICOM Segmentation object (DICOM SEG in short) organizes segmentations in _segments_, where each segment corresponds to a set of labeled image pixels, accompanied by metadata describing what is segmented in that segment. Note that DICOM SEG belongs to the `enhanced multiframe` family of DICOM objects, which means that the DICOM SEG series consists of a single instance / single file, and that single file contains all segments and all frames of each segment.\n",
        "\n",
        "Metadata describing the segments is buried in the nested sequences within the DICOM SEG object. To simplify interaction with the DICOM SEG metadata, IDC maintains a separate metadata table: [`bigquery-open-data.idc_current_segmentations`](https://console.cloud.google.com/bigquery?p=bigquery-public-data&d=idc_current&t=segmentations&page=table).\n",
        "\n",
        "Each of the rows in that table corresponds to a segment. Columns of the table contain coded dictionary values describing the content of the segment, and also include metadata that can be used to link the segment to the patient/study/series it is coming from. All segments within a DICOM SEG are numbered consecutively, and thus `SeriesInstanceUID` + `SegmentNumber` can be used to uniquely identify the segment.\n",
        "\n",
        "We will first get `SeriesInstanceUID` corresponding to the segmentation series.\n",
        "\n",
        "Let's query the segmentations table to get information about the segments that are available."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "435yDsVer4f-"
      },
      "source": [
        "series_selection_query = f\"\\\n",
        "  SELECT \\\n",
        "    DISTINCT(SeriesInstanceUID),\\\n",
        "    SOPInstanceUID, \\\n",
        "    Modality \\\n",
        "  FROM \\\n",
        "    `bigquery-public-data.idc_current.dicom_all` \\\n",
        "  WHERE \\\n",
        "    Modality = \\\"SEG\\\" AND StudyInstanceUID = \\\"{my_StudyInstanceUID}\\\"\"\n",
        "\n",
        "series_selection_result = bq_client.query(series_selection_query)\n",
        "series_selection_df = series_selection_result.result().to_dataframe()\n",
        "\n",
        "import pandas as pd\n",
        "pd.set_option('display.max_colwidth', None)\n",
        "series_selection_df"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4dENCbvwtYjk"
      },
      "source": [
        "From the query result above we can get the `SOPInstanceUID` identifier for the DICOM SEG series, and using that identifier query the table containing segmentation segments metadata to understand what is inside the segmentation series."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nhi_CxPDuPht"
      },
      "source": [
        "seg_SOPInstanceUID = \"1.2.276.0.7230010.3.1.4.2323910823.11504.1597260515.422\"\n",
        "\n",
        "segments_selection_query = f\"\\\n",
        "  SELECT \\\n",
        "    SOPInstanceUID, \\\n",
        "    SegmentNumber, \\\n",
        "    SegmentedPropertyType \\\n",
        "  FROM \\\n",
        "    `bigquery-public-data.idc_current.segmentations` \\\n",
        "  WHERE \\\n",
        "    SOPInstanceUID = \\\"{seg_SOPInstanceUID}\\\"\"\n",
        "\n",
        "segments_selection_result = bq_client.query(segments_selection_query)\n",
        "segments_selection_df = segments_selection_result.result().to_dataframe()\n",
        "\n",
        "segments_selection_df"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eZbunBaNwZQM"
      },
      "source": [
        "As you can see from the result of the query above, instead of a plain text description of the segmentation that we saw for RTSTRUCT, there is a _coded tuple_ (CodeValue, CodingSchemeDesignator, CodeMeaning) describing each segment. Without getting into the details of explaining the content of this tuple, you can use `CodeMeaning` value to establish that Esophagus segmentation corresponds to `SegmentNumber` 1, and `Heart` is in segment number 3.\n",
        "\n",
        "Next let's retrieve the file corresponding to this DICOM SEG, and convert segments contained in it into a NRRD format. To retrieve the file, we will need to look up `gcs_url`, given `SOPInstanceUID` (or `SeriesInstanceUID`, since there is only a single file in the DICOM SEG series), and then retrieve that file from the storage bucket."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BoXldzc1xQ1u"
      },
      "source": [
        "segmentation_selection_query = f\"\\\n",
        "  SELECT \\\n",
        "    gcs_url \\\n",
        "  FROM \\\n",
        "    `bigquery-public-data.idc_current.dicom_all` \\\n",
        "  WHERE \\\n",
        "    SOPInstanceUID = \\\"{seg_SOPInstanceUID}\\\"\"\n",
        "\n",
        "segmentation_selection_result = bq_client.query(segmentation_selection_query)\n",
        "segmention_selection_df = segmentation_selection_result.result().to_dataframe()\n",
        "\n",
        "print(f\"GCS URL for the DICOM SEG file: {segmention_selection_df['gcs_url'][0]}\")"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CIxJF7e9jdIo"
      },
      "source": [
        "# path to the files where the converted NRRD segments will be stored\n",
        "seg_folder = os.path.join(pat_dir_path_nrrd, \"SEG\")\n",
        "os.environ[\"SEG_FOLDER\"] = seg_folder"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "w7Yzy-FSxhlg"
      },
      "source": [
        "!mkdir -p $SEG_FOLDER/segments && rm -rf ${SEG_FOLDER}/*.dcm\n",
        "!gsutil cp gs://idc-open-cr/eff917af-8a2a-42fe-9e12-22bceaac5da9.dcm $SEG_FOLDER\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q_UVG2Puce7b"
      },
      "source": [
        "We will use the [`segimage2itkimage`](https://qiicr.gitbook.io/dcmqi-guide/opening/cmd_tools/seg/segimage2itkimage) tool from [dcmqi](https://github.com/QIICR/dcmqi) to extract individual segments from DICOM SEG series we just downloaded into NRRD format, saving each segment into a separate file. The reason the tool saves those into separate files is because in the general case, DICOM SEG segments may overlap, and NRRD 3D image format does not support overlapping segments."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HgtonbCcdtkR"
      },
      "source": [
        "!segimage2itkimage --inputDICOM ${SEG_FOLDER}/eff917af-8a2a-42fe-9e12-22bceaac5da9.dcm --outputDirectory ${SEG_FOLDER}/segments"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4XQegjKgkv7I"
      },
      "source": [
        "!echo $SEG_FOLDER"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lxuclJU8eQjv"
      },
      "source": [
        "If you now take a look at the content of the `$SEG_FOLDER` folder above, you will see NRRD files named 1 through 6, corresponding to the individual segments stored in DICOM SEG. `meta.json` in that same folder has the segment-level metadata in `segmentAttributes` describing what kind of segmentation is stored in each of the NRRD files. `labelID` attribute corresponds to the segment number, and the numbering of the NRRD files. For example, this part of the JSON file can be used to conclude that segment 3, which corresponds to the file `3.nrrd` contains segmentation of the heart (which is in agreement with the content of the segment-level table we examined earlier):\n",
        "\n",
        "```json\n",
        "         {\n",
        "            \"SegmentAlgorithmType\" : \"MANUAL\",\n",
        "            \"SegmentDescription\" : \"Heart\",\n",
        "            \"SegmentLabel\" : \"Heart\",\n",
        "            \"SegmentedPropertyCategoryCodeSequence\" : {\n",
        "               \"CodeMeaning\" : \"Anatomical Structure\",\n",
        "               \"CodeValue\" : \"123037004\",\n",
        "               \"CodingSchemeDesignator\" : \"SCT\"\n",
        "            },\n",
        "            \"SegmentedPropertyTypeCodeSequence\" : {\n",
        "               \"CodeMeaning\" : \"Heart\",\n",
        "               \"CodeValue\" : \"80891009\",\n",
        "               \"CodingSchemeDesignator\" : \"SCT\"\n",
        "            },\n",
        "            \"labelID\" : 3,\n",
        "            \"recommendedDisplayRGBValue\" : [ 206, 110, 84 ]\n",
        "         }\n",
        "```\n",
        "\n",
        "Similarly, we can establish that `1.nrrd` contains segmentation of the esophagus. We can next load those two segmentations into python arrays and do the comparison visualization. In this case, RTSTRUCT and SEG segmentations are pretty much the same - they are just alternative means to encode the same segmentation result. In other situations, you most likely will encounter segmentations available in only one of those formats."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UZM_bqhmiRqT"
      },
      "source": [
        "# manual segmask (from the RTSTRUCT)\n",
        "seg_segmask_heart = os.path.join(seg_folder, \"segments\", \"3.nrrd\")\n",
        "seg_nrrd_heart = sitk.GetArrayFromImage(sitk.ReadImage(seg_segmask_heart))\n",
        "\n",
        "try:\n",
        "  seg_segmask_esophagus = os.path.join(seg_folder, \"segments\", \"1.nrrd\")\n",
        "  seg_nrrd_esophagus = sitk.GetArrayFromImage(sitk.ReadImage(seg_segmask_esophagus))\n",
        "except:\n",
        "  # for the sake of simplicity, fill the volume with zeros\n",
        "  # (so that we can keep the code that comes after the same)\n",
        "  seg_nrrd_esophagus = np.zeros(seg_nrrd_heart.shape)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "N2gX3VwhlZ4N"
      },
      "source": [
        "_ = viz_utils.AxialSliceSegmaskComparison(ct_volume = ct_nrrd,\n",
        "                                          segmask_ai_dict = {\"Heart\" : pred_nrrd_heart,\n",
        "                                                             \"Aorta\" : pred_nrrd_aorta,\n",
        "                                                             \"Trachea\" : pred_nrrd_trachea,\n",
        "                                                             \"Esophagus\" : pred_nrrd_esophagus},\n",
        "                                          segmask_manual_dict = {\"Heart\" : seg_nrrd_heart,\n",
        "                                                                 \"Esophagus\" : seg_nrrd_esophagus},\n",
        "                                          segmask_cmap_dict = {\"Heart\" : my_greens,\n",
        "                                                               \"Aorta\" : my_spring,\n",
        "                                                               \"Esophagus\" : my_reds,\n",
        "                                                               \"Trachea\" : my_blues},\n",
        "                                          dpi = 100)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MhfNqPxNJ-xz"
      },
      "source": [
        "### Quantitative assessment of the results\n",
        "\n",
        "Let's start by defining a function to compute the center of mass (CoM) of the segmentation masks. Before computing the common segmentation metrics, the CoM can give us a rough idea of how different the 3D delineations are and if there are any major labelling errors (which we could correct, e.g., with a largest connected component analysis).\n",
        "\n",
        "We will base our function on the [implementation](https://github.com/AIM-Harvard/pyradiomics/blob/master/radiomics/generalinfo.py) found in the open source [PyRadiomics library](https://github.com/AIM-Harvard/pyradiomics)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tYF05ymwKMhE"
      },
      "source": [
        "def getCenterOfMassIndexValue(input_mask):\n",
        "    \n",
        "    \"\"\"\n",
        "    Returns z, y and x coordinates of the center of mass of the ROI in terms of\n",
        "    the image coordinate space (continuous index).\n",
        "\n",
        "    Calculation is based on the original (non-resampled) mask.\n",
        "    Because this represents the continuous index, the order of x, y and z is reversed,\n",
        "    i.e. the first element is the z index, the second the y index and the last element is the x index.\n",
        "\n",
        "    @params:\n",
        "      input_mask - required : numpy (binary) volume storing the segmentation mask.\n",
        "\n",
        "    \"\"\"\n",
        "\n",
        "    if input_mask is not None:\n",
        "      mask_coordinates = np.array(np.where(input_mask == 1))\n",
        "      center_index = np.mean(mask_coordinates, axis = 1)\n",
        "      return tuple(center_index)\n",
        "    else:\n",
        "      return None"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Xl_qnwbuKOn8"
      },
      "source": [
        "com_manual_heart = np.array(getCenterOfMassIndexValue(rt_nrrd_heart))\n",
        "com_manual_heart_int = np.ceil(com_manual_heart).astype(dtype = np.uint16)\n",
        "\n",
        "com_raw_heart = np.array(getCenterOfMassIndexValue(pred_nrrd_heart))\n",
        "com_raw_heart_int = np.ceil(com_raw_heart).astype(dtype = np.uint16)\n",
        "\n",
        "print(\"Heart Center of Mass (raw AI segmentation) \\t:\", com_raw_heart_int)\n",
        "print(\"Heart Center of Mass (manual segmentation) \\t:\", com_manual_heart_int)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jnN020q0KS_e"
      },
      "source": [
        "com_manual_heart = np.array(getCenterOfMassIndexValue(rt_nrrd_heart))\n",
        "com_manual_heart_int = np.ceil(com_manual_heart).astype(dtype = np.uint16)\n",
        "\n",
        "com_raw_heart = np.array(getCenterOfMassIndexValue(pred_nrrd_heart))\n",
        "com_raw_heart_int = np.ceil(com_raw_heart).astype(dtype = np.uint16)\n",
        "\n",
        "print(\"Heart Center of Mass (raw AI segmentation) \\t:\", com_raw_heart_int)\n",
        "print(\"Heart Center of Mass (manual segmentation) \\t:\", com_manual_heart_int)\n",
        "\n",
        "## ----------------------------------------\n",
        "\n",
        "# run this if and only if a manual esophagus segmentation mask is available\n",
        "if np.sum(rt_nrrd_esophagus):\n",
        "\n",
        "  com_manual_esophagus = np.array(getCenterOfMassIndexValue(rt_nrrd_esophagus))\n",
        "  com_manual_esophagus_int = np.ceil(com_manual_esophagus).astype(dtype = np.uint16)\n",
        "\n",
        "  com_raw_esophagus = np.array(getCenterOfMassIndexValue(pred_nrrd_esophagus))\n",
        "  com_raw_esophagus_int = np.ceil(com_raw_esophagus).astype(dtype = np.uint16)\n",
        "\n",
        "  print(\"\\nEsophagus Center of Mass (raw AI segmentation) \\t:\", com_raw_esophagus_int)\n",
        "  print(\"Esophagus Center of Mass (manual segmentation) \\t:\", com_manual_esophagus_int)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rYBAo8SzKW31"
      },
      "source": [
        "Another common way to evaluate the quality of the segmentation is computing the Dice Coefficient between the AI segmentation and the manual one. To do so, we will use [MedPy's implementation of the Dice coefficient](https://loli.github.io/medpy/generated/medpy.metric.binary.dc.html#medpy-metric-binary-dc) (for binary masks).\n",
        "\n",
        "We can use other MedPy's functions to compute the Hausdorff distance and the average surface distance as well*.\n",
        "\n",
        "_*in most cases, the Hausdorff Distance will be quite high for both the heart segmentation and, if available with the randomly selected patient, the esophagus one. This is not a clear indication the model performance is poor: rather, it could also be the segmentation guidelines of the two datasets (the one the nnU-Net model was trained on and the external and independent validation dataset pulled from IDC) differ significantly._"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "k1xANs3rKWTo"
      },
      "source": [
        "pred_nrrd_path = os.path.join(pat_dir_path_nrrd, pat + \"_pred_segthor.nrrd\")\n",
        "\n",
        "\n",
        "voxel_spacing = list(sitk_ct.GetSpacing())\n",
        "\n",
        "dc_heart = dice_coef(pred_nrrd_heart, rt_nrrd_heart)\n",
        "hd_heart = hausdorff_distance(pred_nrrd_heart, rt_nrrd_heart, voxelspacing = voxel_spacing)\n",
        "asd_heart = avg_surf_distance(pred_nrrd_heart, rt_nrrd_heart, voxelspacing = voxel_spacing)\n",
        "\n",
        "print(\"Heart Dice Coefficient (raw segmentation) :\", dc_heart)\n",
        "print(\"Heart Hausdorff Distance (raw segmentation) [mm]:\", hd_heart)\n",
        "print(\"Heart Average Surface Distance (raw segmentation) [mm]:\", asd_heart)\n",
        "\n",
        "\n",
        "# run this if and only if a manual esophagus segmentation mask is available\n",
        "if np.sum(rt_nrrd_esophagus):\n",
        "  dc_esophagus = dice_coef(pred_nrrd_esophagus, rt_nrrd_esophagus)\n",
        "  hd_esophagus = hausdorff_distance(pred_nrrd_esophagus, rt_nrrd_esophagus, voxelspacing = voxel_spacing)\n",
        "  asd_esophagus = avg_surf_distance(pred_nrrd_esophagus, rt_nrrd_esophagus, voxelspacing = voxel_spacing)\n",
        "\n",
        "  print(\"\\nEsophagus Dice Coefficient (raw segmentation) :\", dc_esophagus)\n",
        "  print(\"Esophagus Hausdorff Distance (raw segmentation) [mm]:\", hd_esophagus)\n",
        "  print(\"Esophagus Average Surface Distance (raw segmentation) [mm]:\", asd_esophagus)\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6ehkQblEkgxe"
      },
      "source": [
        "## Finding more relevant images to test\n",
        "\n",
        "Let's get some better idea of the CT series we segmented. We can query IDC BigQuery DICOM metadata table to retrieve some relevant attributes.\n",
        "\n",
        "```sql\n",
        "SELECT\n",
        "  DISTINCT(Manufacturer),\n",
        "  ManufacturerModelName\n",
        "FROM\n",
        " `bigquery-public-data.idc_current.dicom_all`\n",
        "WHERE\n",
        " SeriesInstanceUID = <my_SeriesInstanceUID>\n",
        "```"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "r7-hRJBd05It"
      },
      "source": [
        "details_selection_query = f\"\\\n",
        "  SELECT \\\n",
        "    DISTINCT(Manufacturer), \\\n",
        "    ManufacturerModelName, \\\n",
        "    BodyPartExamined \\\n",
        "  FROM \\\n",
        "    `bigquery-public-data.idc_current.dicom_all` \\\n",
        "  WHERE \\\n",
        "    SeriesInstanceUID = \\\"{my_SeriesInstanceUID}\\\"\"\n",
        "\n",
        "details_selection_result = bq_client.query(details_selection_query)\n",
        "details_selection_df = details_selection_result.result().to_dataframe()\n",
        "\n",
        "details_selection_df"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SYvs4Q6A1SsF"
      },
      "source": [
        "We can easily check what other combinations of `Manufacturer` and `ManufacturerModelName` we have in store, to test generalizability of the segmentation network we have been using. In the query below, we will get all unique combinations of `Manufacturer`/`ManufacturerModelName`, and count how many CT studies that have `BodyPartExamined = LUNG` are available for that specific combination.\n",
        "\n",
        "In a similar fashion, you can utilize any of the DICOM attributes to select representative data to use in testing.\n",
        "\n",
        "Note that such explorations can also be done using [IDC portal](https://imaging.datacommons.cancer.gov/), [IDC DataStudio dashboard](https://datastudio.google.com/c/reporting/04cf5976-4ea0-4fee-a749-8bfd162f2e87), or [TCIA radiology portal](https://nbia.cancerimagingarchive.net/nbia-search/). Exploratory SQL queries to IDC tables are best done in the [BigQuery console](https://console.cloud.google.com/bigquery).\n",
        "\n",
        "```sql\n",
        "SELECT\n",
        "  Manufacturer,\n",
        "  ARRAY_TO_STRING(ARRAY_AGG(DISTINCT(ManufacturerModelName)),',') AS ManufacturerModelNames,\n",
        "  ARRAY_TO_STRING(ARRAY_AGG(DISTINCT(ARRAY_TO_STRING(ConvolutionKernel,','))),',') AS ConvolutionKernels,\n",
        "  COUNT(DISTINCT(StudyInstanceUID)) AS number_of_studies\n",
        "FROM\n",
        "  `bigquery-public-data.idc_current.dicom_all`\n",
        "WHERE\n",
        "  Modality = \"CT\"\n",
        "  AND BodyPartExamined = \"LUNG\"\n",
        "GROUP BY\n",
        "  Manufacturer\n",
        "ORDER BY\n",
        "  number_of_studies DESC\n",
        "```"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "My6L8D1k2Aoq"
      },
      "source": [
        "sample_selection_query = f\"\\\n",
        "SELECT \\\n",
        "  Manufacturer, \\\n",
        "  ARRAY_TO_STRING(ARRAY_AGG(DISTINCT(ManufacturerModelName)),',') as ManufacturerModelNames, \\\n",
        "  COUNT(DISTINCT(StudyInstanceUID)) AS number_of_studies \\\n",
        "FROM \\\n",
        "  `bigquery-public-data.idc_current.dicom_all` \\\n",
        "WHERE \\\n",
        "  Modality = \\\"CT\\\" \\\n",
        "  AND BodyPartExamined = \\\"LUNG\\\" \\\n",
        "GROUP BY \\\n",
        "  Manufacturer \\\n",
        "ORDER BY \\\n",
        "  number_of_studies DESC\"\n",
        "\n",
        "sample_selection_result = bq_client.query(sample_selection_query)\n",
        "sample_selection_df = sample_selection_result.result().to_dataframe()\n",
        "\n",
        "sample_selection_df"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BOMn1YcH2iu5"
      },
      "source": [
        "Next we can query for a study that has `Manufacturer`/`ManufacturerModelName` combination we want to test, and visualize that study. \n",
        "\n",
        "Use [BigQuery console](https://console.cloud.google.com/bigquery) to run this queries!\n",
        "\n",
        "\n",
        "```sql\n",
        "SELECT\n",
        "  SeriesInstanceUID,\n",
        "  ARRAY_TO_STRING(ARRAY_AGG(DISTINCT(SeriesDescription)),',') AS SeriesDescriptions,\n",
        "  ARRAY_TO_STRING(ARRAY_AGG(DISTINCT(StudyInstanceUID)),',') AS StudyInstanceUIDs,  \n",
        "  ARRAY_TO_STRING(ARRAY_AGG(DISTINCT(ARRAY_TO_STRING(ConvolutionKernel,','))),',') AS ConvolutionKernels,\n",
        "  COUNT(DISTINCT(SOPInstanceUID)) AS numberOfInstances\n",
        "FROM\n",
        "  `bigquery-public-data.idc_current.dicom_all`\n",
        "WHERE\n",
        "  Modality = \"CT\"\n",
        "  AND BodyPartExamined = \"LUNG\"\n",
        "  AND Manufacturer = \"GE MEDICAL SYSTEMS\"\n",
        "  AND ManufacturerModelName = \"LightSpeed Xtra\"\n",
        "GROUP BY\n",
        "  SeriesInstanceUID\n",
        "```\n",
        "\n",
        "This query returns three `SeriesInstanceUID`s.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-Gd4RSiT53zV"
      },
      "source": [
        "print(get_idc_viewer_url(\"1.3.6.1.4.1.14519.5.2.1.3023.4012.146358232297157912733174589662\"))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jKky12H46oaU"
      },
      "source": [
        "Let's re-initialize `SeriesInstanceUID` with the series from this study, and re-run the segmentation steps by returning to the **Visualization and download of data from IDC** section of the notebook after running the cell below!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QxYJGUvC6rFI"
      },
      "source": [
        "my_SeriesInstanceUID = \"1.3.6.1.4.1.14519.5.2.1.3023.4012.162275549801143329076803363880\""
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rUVGRNxZQerC"
      },
      "source": [
        "## I want to train my network, not run inference - what do I do?\n",
        "\n",
        "Google Colab has a number of limitations on the resources available and time you can use this free resource. If you plan to train deep learning models on volumetric medical imaging data using cloud VMs, you may consider GCP [Vertex AI Notebooks](https://cloud.google.com/vertex-ai-workbench) as an alternative to Google Colab. AI Notebooks is a paid resource that you can use to define highly configurable AI development environments accessible via JupyterLab interface. You can [apply for a sponsored project from IDC](https://learn.canceridc.dev/introduction/requesting-gcp-cloud-credits) to experiment with those capabilities at no cost to you."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XHQSJ-jFQnFD"
      },
      "source": [
        "## References\n",
        "\n",
        "Clark, K., Vendt, B., Smith, K., Freymann, J., Kirby, J., Koppel, P., Moore, S., Phillips, S., Maffitt, D., Pringle, M., Tarbox, L. & Prior, F. The Cancer Imaging Archive (TCIA): maintaining and operating a public information repository. J. Digit. Imaging 26, 1045–1057 (2013). http://dx.doi.org/10.1007/s10278-013-9622-7\n",
        "\n",
        "Fedorov, A., Longabaugh, W. J. R., Pot, D., Clunie, D. A., Pieper, S., Aerts, H. J. W. L., Homeyer, A., Lewis, R., Akbarzadeh, A., Bontempi, D., Clifford, W., Herrmann, M. D., Höfener, H., Octaviano, I., Osborne, C., Paquette, S., Petts, J., Punzo, D., Reyes, M., Schacherer, D. P., Tian, M., White, G., Ziegler, E., Shmulevich, I., Pihl, T., Wagner, U., Farahani, K. & Kikinis, R. NCI Imaging Data Commons. Cancer Res. 81, 4188–4193 (2021). http://dx.doi.org/10.1158/0008-5472.CAN-21-0950\n",
        "* You can also review the [demonstration videos](https://www.youtube.com/playlist?list=PLhawVWNiPvwb2H0D9UTOIL23bc5DuDJRu) accompanying the manuscript above for a quick summary of some of the key features of the platform."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8nSlEEDVQLq4"
      },
      "source": [
        "## Acknowledgments\n",
        "\n",
        "The Cancer Imaging Archive has been funded in whole or in part with federal funds from the National Cancer Institute, National Insti-tutes of Health, under Contract No. HHSN261200800001E.\n",
        "\n",
        "NCI Imaging Data Commons has been funded in whole or in part with Federal funds from the National Cancer Institute, National Institutes of Health, under Task Order No.HHSN26110071 under Contract No. HHSN261201500003l."
      ]
    }
  ]
}