10_results_visualization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\">Visualization of Segmentation and Registration Results</h1>\n",
    "\n",
    "In this notebook we illustrate various ways one can display the results of segmentation and registration algorithms so that they can be easily incorporated into a manuscript. For interactive data exploration we recommend using dedicated programs (e.g. 3D slicer). \n",
    "\n",
    "Two key points to remember when working with bio-medical images:\n",
    "\n",
    "1. Most often images have a high dynamic range. Thus, to write them to file in a format appropriate for use in a manuscript we will need to map the intensities to a low dynamic range (e.g. [0,255]). In SimpleITK this is readily done with the [IntensityWindowingImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1IntensityWindowingImageFilter.html).\n",
    "2. Images may have non-isotropic spacing between pixels. The file formats appropriate for use in a manuscript (e.g. png, jpg) assume isotropic pixel spacing. This requires that we resample the image before writing to disk. The function `make_isotropic` in the code cell bellow resolves this issue. \n",
    "\n",
    "The following filters and their procedural counterparts are useful for various image creation tasks, as illustrated in this notebook:\n",
    " * [CastImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1CastImageFilter.html)\n",
    " * [ResampleImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1ResampleImageFilter.html), one of the more important filters in your toolbox.\n",
    " * [TileImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1TileImageFilter.html)\n",
    " * [CheckerBoardImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1CheckerBoardImageFilter.html)\n",
    " * [ComposeImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1ComposeImageFilter.html)\n",
    " * [LabelToRGBImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1LabelToRGBImageFilter.html)\n",
    " * [ScalarToRGBColormapImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1ScalarToRGBColormapImageFilter.html)\n",
    " * [LabelOverlayImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1LabelOverlayImageFilter.html)\n",
    " * [LabelContourImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1LabelContourImageFilter.html)\n",
    " * [LabelMapContourOverlayImageFilter](https://simpleitk.org/doxygen/latest/html/namespaceitk_1_1simple.html#a4f6af69f85171e44bcff90d7860d456e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import SimpleITK as sitk\n",
    "\n",
    "# Utility method that either downloads data from the Girder repository or\n",
    "# if already downloaded returns the file name for reading from disk (cached data).\n",
    "from downloaddata import fetch_data as fdata\n",
    "\n",
    "import utilities\n",
    "import gui\n",
    "\n",
    "# Always write output to a separate directory, we don't want to pollute the source directory.\n",
    "import os\n",
    "\n",
    "OUTPUT_DIR = \"output\"\n",
    "\n",
    "\n",
    "def make_isotropic(image, interpolator=sitk.sitkLinear):\n",
    "    \"\"\"\n",
    "    Resample an image to isotropic pixels (using smallest spacing from original) and save to file. Many file formats\n",
    "    (jpg, png,...) expect the pixels to be isotropic. By default the function uses a linear interpolator. For\n",
    "    label images one should use the sitkNearestNeighbor interpolator so as not to introduce non-existant labels.\n",
    "    \"\"\"\n",
    "    original_spacing = image.GetSpacing()\n",
    "    # Image is already isotropic, just return a copy.\n",
    "    if all(spc == original_spacing[0] for spc in original_spacing):\n",
    "        return sitk.Image(image)\n",
    "    # Make image isotropic via resampling.\n",
    "    original_size = image.GetSize()\n",
    "    min_spacing = min(original_spacing)\n",
    "    new_spacing = [min_spacing] * image.GetDimension()\n",
    "    new_size = [\n",
    "        int(round(osz * ospc / min_spacing))\n",
    "        for osz, ospc in zip(original_size, original_spacing)\n",
    "    ]\n",
    "    return sitk.Resample(\n",
    "        image,\n",
    "        new_size,\n",
    "        sitk.Transform(),\n",
    "        interpolator,\n",
    "        image.GetOrigin(),\n",
    "        new_spacing,\n",
    "        image.GetDirection(),\n",
    "        0,\n",
    "        image.GetPixelID(),\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Combining two images\n",
    "\n",
    "There are a variety of ways we can overlay two (partially) overlapping images onto each other. The common approaches include:\n",
    "1. Use of alpha blending.\n",
    "2. Use of a checkerboard pattern with the pixel values in adjacent squares/boxes taken from each of the images.\n",
    "3. When the pixel values are scalars (gray scale images), combine the two images in different channels, resulting in a color image.\n",
    "\n",
    "We will start by loading two images whose content luckily overlaps in physical space. Before we can combine the two, we need to resample one of them so that they both occupy the same spatial region. In addition we should also rescale the intensities so that they occupy the same range. In our case we will map them to [0,255], based on the desired windowing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img1 = sitk.ReadImage(fdata(\"training_001_mr_T1.mha\"))\n",
    "img2_original = sitk.ReadImage(fdata(\"training_001_ct.mha\"))\n",
    "img2 = sitk.Resample(img2_original, img1)\n",
    "\n",
    "# Obtain foreground masks for the two images using Otsu thresholding, we use these later on.\n",
    "msk1 = sitk.OtsuThreshold(img1, 0, 1)\n",
    "msk2 = sitk.OtsuThreshold(img2, 0, 1)\n",
    "\n",
    "gui.MultiImageDisplay(\n",
    "    image_list=[img1, img2], title_list=[\"image1\", \"image2\"], figure_size=(9, 3)\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Having identified the desired intensity range for each of the\n",
    "# images using the GUI above, we use these values to perform intensity windowing and map the intensity values\n",
    "# to [0,255] and cast to 8-bit unsigned int\n",
    "img1_255 = sitk.Cast(\n",
    "    sitk.IntensityWindowing(\n",
    "        img1, windowMinimum=2, windowMaximum=657, outputMinimum=0.0, outputMaximum=255.0\n",
    "    ),\n",
    "    sitk.sitkUInt8,\n",
    ")\n",
    "img2_255 = sitk.Cast(\n",
    "    sitk.IntensityWindowing(\n",
    "        img2,\n",
    "        windowMinimum=-1018,\n",
    "        windowMaximum=1126,\n",
    "        outputMinimum=0.0,\n",
    "        outputMaximum=255.0,\n",
    "    ),\n",
    "    sitk.sitkUInt8,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alpha blending\n",
    "\n",
    "Alpha blending combines the pixels from the two images as follows:\n",
    "$$\n",
    "I_{output} = \\alpha I_1 + (1-\\alpha)I_2,\\;\\;\\; \\alpha \\in[0.0,1.0]\n",
    "$$\n",
    "\n",
    "When our images consist of a foreground and background we can use alpha blending in a manner that takes this into account. Instead of blending all of the pixels using the formula above, we use this formula only in the regions where the foregrounds overlap. In regions where the foreground from one image overlaps with the background of the other we simply copy the foreground. This improves visibility as we are not blending a region that contains information with an empty region.\n",
    "\n",
    "The code below allows us to experiment with various alpha blending strategies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mask_image_multiply(mask, image):\n",
    "    components_per_pixel = image.GetNumberOfComponentsPerPixel()\n",
    "    if components_per_pixel == 1:\n",
    "        return mask * image\n",
    "    else:\n",
    "        return sitk.Compose(\n",
    "            [\n",
    "                mask * sitk.VectorIndexSelectionCast(image, channel)\n",
    "                for channel in range(components_per_pixel)\n",
    "            ]\n",
    "        )\n",
    "\n",
    "\n",
    "def alpha_blend(image1, image2, alpha=0.5, mask1=None, mask2=None):\n",
    "    \"\"\"\n",
    "    Alaph blend two images, pixels can be scalars or vectors.\n",
    "    The region that is alpha blended is controled by the given masks.\n",
    "    \"\"\"\n",
    "\n",
    "    if not mask1:\n",
    "        mask1 = sitk.Image(image1.GetSize(), sitk.sitkFloat32) + 1.0\n",
    "        mask1.CopyInformation(image1)\n",
    "    else:\n",
    "        mask1 = sitk.Cast(mask1, sitk.sitkFloat32)\n",
    "    if not mask2:\n",
    "        mask2 = sitk.Image(image2.GetSize(), sitk.sitkFloat32) + 1\n",
    "        mask2.CopyInformation(image2)\n",
    "    else:\n",
    "        mask2 = sitk.Cast(mask2, sitk.sitkFloat32)\n",
    "\n",
    "    components_per_pixel = image1.GetNumberOfComponentsPerPixel()\n",
    "    if components_per_pixel > 1:\n",
    "        img1 = sitk.Cast(image1, sitk.sitkVectorFloat32)\n",
    "        img2 = sitk.Cast(image2, sitk.sitkVectorFloat32)\n",
    "    else:\n",
    "        img1 = sitk.Cast(image1, sitk.sitkFloat32)\n",
    "        img2 = sitk.Cast(image2, sitk.sitkFloat32)\n",
    "\n",
    "    intersection_mask = mask1 * mask2\n",
    "\n",
    "    intersection_image = mask_image_multiply(\n",
    "        alpha * intersection_mask, img1\n",
    "    ) + mask_image_multiply((1 - alpha) * intersection_mask, img2)\n",
    "    return (\n",
    "        intersection_image\n",
    "        + mask_image_multiply(mask2 - intersection_mask, img2)\n",
    "        + mask_image_multiply(mask1 - intersection_mask, img1)\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now create 3D images using all four combinations of alpha-blending and masks. As we are working with a 3D image and we want to save it as a figure for use in a manuscript, we will create a 2D montage image using the axial slices from the volumes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combine the two volumes\n",
    "images_list = [\n",
    "    (alpha_blend(img1_255, img2_255), \"alpha_blend_standard\"),\n",
    "    (alpha_blend(img1_255, img2_255, mask1=msk1), \"alpha_blend_mask1\"),\n",
    "    (alpha_blend(img1_255, img2_255, mask2=msk2), \"alpha_blend_mask2\"),\n",
    "    (\n",
    "        alpha_blend(img1_255, img2_255, mask1=msk1, mask2=msk2),\n",
    "        \"alpha_blend_mask1_mask2\",\n",
    "    ),\n",
    "]\n",
    "\n",
    "# Tile the volumes using the x-y plane (axial slices)\n",
    "all_montages = []\n",
    "for img, img_name in images_list:\n",
    "    num_slices = img.GetDepth()\n",
    "    tile_w = int(np.sqrt(num_slices))\n",
    "    tile_h = int(np.ceil(num_slices / tile_w))\n",
    "    tile_image = sitk.Tile([img[:, :, i] for i in range(num_slices)], (tile_w, tile_h))\n",
    "    sitk.WriteImage(\n",
    "        sitk.Cast(tile_image, sitk.sitkUInt8),\n",
    "        os.path.join(OUTPUT_DIR, img_name + \".png\"),\n",
    "    )\n",
    "    all_montages.append(tile_image)\n",
    "\n",
    "# Display all montages by combining them into a faux volume. Notice that scrolling through this\n",
    "# volume creates the illusion of motion due to the change in intensities (the interested\n",
    "# reader is referred to \"Visual dissociations of movement, position, and stereo depth: Some phenomenal\n",
    "# phenomena\", R. L. Gregory, P. F. Heard).\n",
    "gui.MultiImageDisplay(\n",
    "    image_list=[sitk.JoinSeries(all_montages)],\n",
    "    title_list=[\"Montages With Different Alpha Blending Strategies\"],\n",
    "    figure_size=(8, 8),\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checkerboard\n",
    "\n",
    "Combine the original and the intensity windowed images using the checkerboard pattern. This illustrates the need to map both images to the same intensity range. This is particularly visible in the background region, where both images contain air. \n",
    "\n",
    "You can specify the number of checkerboard tiles per dimension as illustrated below. The actual number of checkerboard tiles depends on the number of pixels per dimension and the specified number of tiles. You may get more tiles than specified, for example try specifying [4,4,7] below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gui.MultiImageDisplay(\n",
    "    image_list=[\n",
    "        sitk.CheckerBoard(img1, img2, [4, 4, 4]),\n",
    "        sitk.CheckerBoard(img1_255, img2_255, (10, 10, 4)),\n",
    "    ],\n",
    "    title_list=[\"original intensities\", \"rescaled intensities\"],\n",
    "    figure_size=(9, 3),\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combine scalar images to create color image\n",
    "\n",
    "There are a variety of ways we can combine the scalar images to create a color image. Some of these combinations should be avoided as they are not discernible by a significant portion of the population (i.e. red-green channel encoding). For additional details see:\n",
    "\n",
    "M. Geissbuehler, T. Lasser, \"How to display data by color schemes compatible with red-green color perception deficiencies\", Opt Express., 21(8):9862-74, 2013.  \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zeros = sitk.Image(img1_255.GetSize(), img1_255.GetPixelID())\n",
    "zeros.CopyInformation(img1_255)\n",
    "\n",
    "gui.MultiImageDisplay(\n",
    "    image_list=[\n",
    "        sitk.Cast(sitk.Compose(img1_255, img2_255, zeros), sitk.sitkVectorUInt8),\n",
    "        sitk.Cast(sitk.Compose(img1_255, img2_255, img1_255), sitk.sitkVectorUInt8),\n",
    "        sitk.Cast(\n",
    "            sitk.Compose(img1_255, 0.5 * img1_255 + 0.5 * img2_255, img2_255),\n",
    "            sitk.sitkVectorUInt8,\n",
    "        ),\n",
    "    ],\n",
    "    title_list=[\"avoid red-green\", \"use magenta-green\", \"use orange-blue\"],\n",
    "    figure_size=(9, 3),\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overlay scalar image onto another via color map\n",
    "\n",
    "In some situations we have an underlying image (spatial structures) and we want to overlay a scalar based pseudo-color image on top of it.\n",
    "\n",
    "This is relevant for presenting co-registered PET/CT data, with the PET providing functional information and the CT providing the underlying spatial structures. A similar use case in the context of deep learning is to display activation maps illustrating the regions in an image on which the network is focusing its attention for the particular class.\n",
    "\n",
    "The two main decisions we make are:\n",
    "1. Selection of pseudo-color scheme using the [ScalarToRGBColormapImageFilter](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1ScalarToRGBColormapImageFilter.html) which supports a variety of color maps.\n",
    "2. Alpha blending approach (alpha value and usage of masks as done above).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a signed distance map which we will overlay onto the original image using\n",
    "# pseudo-coloring. We are only interested in locations that are at a distance of [0,512] from the object\n",
    "distance_map = sitk.SignedMaurerDistanceMap(msk1)\n",
    "# Get binary region of interest mask.\n",
    "roi = sitk.Cast(distance_map > 0.0, sitk.sitkFloat32) * sitk.Cast(\n",
    "    distance_map < 512.0, sitk.sitkFloat32\n",
    ")\n",
    "roi_distance_map = roi * distance_map\n",
    "\n",
    "overlay_color_img = sitk.ScalarToRGBColormap(\n",
    "    roi_distance_map, sitk.ScalarToRGBColormapImageFilter.Jet\n",
    ")\n",
    "\n",
    "# Combine the color overlay volume with the spatial structure volume using alpha blending\n",
    "# and cast to a three component vector 8 bit unsigned int. We can readily incorporate a\n",
    "# mask into the blend (pun intended). By adding mask2=roi we can limit the overlay to\n",
    "# the region of interest. We can also use a per pixel alpha blending by normalizing the\n",
    "# distance map in our ROI.\n",
    "#\n",
    "# To see the various options, re-run this cell and change the blending_approach (basic, mask, per_pixel_alpha)\n",
    "blending_approach = \"basic\"\n",
    "\n",
    "if blending_approach == \"basic\":\n",
    "    combined_volume = sitk.Cast(\n",
    "        alpha_blend(sitk.Compose(img1_255, img1_255, img1_255), overlay_color_img),\n",
    "        sitk.sitkVectorUInt8,\n",
    "    )\n",
    "elif blending_approach == \"mask\":\n",
    "    combined_volume = sitk.Cast(\n",
    "        alpha_blend(\n",
    "            sitk.Compose(img1_255, img1_255, img1_255), overlay_color_img, mask2=roi\n",
    "        ),\n",
    "        sitk.sitkVectorUInt8,\n",
    "    )\n",
    "elif blending_approach == \"per_pixel_alpha\":\n",
    "    mmFilter = sitk.MinimumMaximumImageFilter()\n",
    "    mmFilter.Execute(roi_distance_map)\n",
    "    alpha = sitk.Cast(\n",
    "        1.0\n",
    "        - (roi_distance_map - mmFilter.GetMinimum())\n",
    "        / (mmFilter.GetMaximum() - mmFilter.GetMinimum()),\n",
    "        sitk.sitkFloat32,\n",
    "    )\n",
    "    combined_volume = sitk.Cast(\n",
    "        alpha_blend(\n",
    "            sitk.Compose(img1_255, img1_255, img1_255), overlay_color_img, alpha=alpha\n",
    "        ),\n",
    "        sitk.sitkVectorUInt8,\n",
    "    )\n",
    "\n",
    "# Given a volume we can either create a montage as above or we can take a representative\n",
    "# slice (axial/sagittal/coronal). As image formats used in manuscripts assume isotropic\n",
    "# pixels we need to ensure this before we write to disk.\n",
    "all_central_slices = [\n",
    "    combined_volume[:, :, int(combined_volume.GetDepth() / 2.0 + 0.5)],\n",
    "    combined_volume[:, int(combined_volume.GetHeight() / 2.0 + 0.5), :],\n",
    "    combined_volume[int(combined_volume.GetWidth() / 2.0 + 0.5), :, :],\n",
    "]\n",
    "\n",
    "# Resample to isotropic pixels and write to file.\n",
    "for i, img in enumerate(all_central_slices):\n",
    "    all_central_slices[i] = make_isotropic(img)\n",
    "    sitk.WriteImage(\n",
    "        all_central_slices[i],\n",
    "        os.path.join(OUTPUT_DIR, \"color_overlay{0}.png\".format(i)),\n",
    "    )\n",
    "\n",
    "gui.multi_image_display2D(\n",
    "    [sitk.Tile(all_central_slices, (1, 3))], figure_size=(4, 4), horizontal=False\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Combining an image and segmentation\n",
    "\n",
    "To display the results of segmentation in context, we need to overlay them onto the original image. There are two common options for doing this:\n",
    "  1. Map the segmentation labels to a color image and alpha blend onto the original image.\n",
    "  2. Overlay the segmentation boundaries onto the original image.\n",
    "\n",
    "We illustrate both approaches below.\n",
    "\n",
    "For this example we use the Point-validated Pixel-based Breathing Thorax Model (POPI) model. The model is provided by the Léon Bérard Cancer Center & CREATIS Laboratory, Lyon, France. The relevant publication is:\n",
    "\n",
    "J. Vandemeulebroucke, D. Sarrut, P. Clarysse, \"The POPI-model, a point-validated pixel-based breathing thorax model\", Proc. XVth International Conference on the Use of Computers in Radiation Therapy (ICCR), Toronto, Canada, 2007."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "img = sitk.ReadImage(fdata(\"POPI/meta/00-P.mhd\"))\n",
    "segmentation = sitk.ReadImage(fdata(\"POPI/masks/00-air-body-lungs.mhd\"))\n",
    "\n",
    "gui.MultiImageDisplay(\n",
    "    image_list=[img, segmentation, sitk.LabelToRGB(segmentation)],\n",
    "    title_list=[\"image\", \"raw segmentation labels\", \"segmentation labels in color\"],\n",
    "    figure_size=(9, 3),\n",
    "    shared_slider=True,\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will work with the central coronal slice from this dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Identify the desired intensity range for our image using the GUI above, mapping the high dynamic range\n",
    "# image to a low dynamic range, [0,255], extract the central coronal slice and we flip it for display purposes.\n",
    "coronal_255 = sitk.Cast(\n",
    "    sitk.IntensityWindowing(\n",
    "        img[:, int(img.GetHeight() / 2), :][:, ::-1],\n",
    "        windowMinimum=-1000,\n",
    "        windowMaximum=170,\n",
    "        outputMinimum=0.0,\n",
    "        outputMaximum=255.0,\n",
    "    ),\n",
    "    sitk.sitkUInt8,\n",
    ")\n",
    "coronal_255_isotropic = make_isotropic(coronal_255)\n",
    "\n",
    "coronal_segmentation = segmentation[:, int(segmentation.GetHeight() / 2), :][:, ::-1]\n",
    "# Use nearest neighbor interpolation for a label image.\n",
    "coronal_segmentation_isotropic = make_isotropic(\n",
    "    coronal_segmentation, sitk.sitkNearestNeighbor\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In many cases the values in a label image are not appropriate for direct display. For instance the values in our\n",
    "label image are 0,1,2. You can see the values if you hover your cursor over the raw segmentation label image above (figure's bottom right corner). \n",
    "\n",
    "In theory we could map these intensities to [0,255] and save the image. In practice we may have more than 256 labels and therefor it is preferable to map the labels to colors and save the color image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the default color map when mapping labels to colors and write the image.\n",
    "sitk.WriteImage(\n",
    "    sitk.LabelToRGB(coronal_segmentation_isotropic),\n",
    "    os.path.join(OUTPUT_DIR, \"coronal_segmentation.png\"),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overlay segmentation labels onto original image\n",
    "\n",
    "SimpleITK allows you to overlay the segmentation labels onto the original image using a color map and alpha blending. You can specify the value for alpha blending, the color map (there is a default color map), and the background label value which will not be overlaid with a label.\n",
    "\n",
    "The color map in SimpleITK is a set of values in the RGB color space strung together. For example [255, 0, 0, 0, 255, 0] is a two entry color map with red and green. To create a human readable color map, use lists to represent the colors and string them together, as done below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Overlay the segmentation using default color map and an alpha value of 0.5\n",
    "coronal_combined1 = sitk.LabelOverlay(\n",
    "    image=coronal_255_isotropic,\n",
    "    labelImage=coronal_segmentation_isotropic,\n",
    "    opacity=0.5,\n",
    "    backgroundValue=utilities.popi_air_label,\n",
    ")\n",
    "\n",
    "# Create an \"interesting\" color map and specify backgroundValue to a non existent label\n",
    "# so that the background label is also overlaid.\n",
    "pink = [255, 105, 180]\n",
    "green = [0, 255, 0]\n",
    "gold = [255, 215, 0]\n",
    "coronal_combined2 = sitk.LabelOverlay(\n",
    "    image=coronal_255_isotropic,\n",
    "    labelImage=coronal_segmentation_isotropic,\n",
    "    opacity=0.5,\n",
    "    backgroundValue=-1.0,\n",
    "    colormap=pink + green + gold,\n",
    ")\n",
    "\n",
    "# Display the two images as a faux volume, JoinSeries, approach.\n",
    "gui.MultiImageDisplay(\n",
    "    image_list=[sitk.JoinSeries([coronal_combined1, coronal_combined2])],\n",
    "    title_list=[\"overlaid labels\"],\n",
    "    figure_size=(9, 3),\n",
    "    shared_slider=True,\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overlay segmentation boundaries onto original image\n",
    "\n",
    "We can readily obtain the segmentation boundaries from the raw segmentation. We then either just save the contours as an image or we can directly overly them onto the image. \n",
    "\n",
    "Some points to note:\n",
    " 1. When working with 3D images and segmentations, our boundaries are surfaces. When these surfaces are intersected  with 2D planes they may define a region and not a contour, which is what we usually expect (e.g. slice 24 in the results displayed by the following code cell).\n",
    " 2. When the labels are next to each other, they share a boundary. As a consequence, drawing the boundaries may result in contours overwriting each other or in contour crossings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "red = [255, 0, 0]\n",
    "green = [0, 255, 0]\n",
    "blue = [0, 0, 255]\n",
    "\n",
    "# red goes to the first label, green to second, blue to third\n",
    "# utilities.popi_body_label=0, utilities.popi_air_label=1, utilities.popi_lung_label=2\n",
    "contour_image = sitk.LabelToRGB(\n",
    "    sitk.LabelContour(segmentation, fullyConnected=True, backgroundValue=255),\n",
    "    colormap=red + green + blue,\n",
    "    backgroundValue=255,\n",
    ")\n",
    "gui.MultiImageDisplay(image_list=[contour_image], figure_size=(9, 3));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following code cell we overlay the contours onto the original image. We use a filter that is intended for usage with a label-map image as opposed to a label image. The former is a more efficient data structure for representing segmentations (run length encoded). We therefor need to cast the label image to a label-map image.\n",
    "We can also set several visualization related parameters such as overlay opacity, contour thickness, priority of overlay (which label overwrites which if they overlap) etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "contour_overlaid_image = sitk.LabelMapContourOverlay(\n",
    "    sitk.Cast(coronal_segmentation_isotropic, sitk.sitkLabelUInt8),\n",
    "    coronal_255_isotropic,\n",
    "    opacity=1,\n",
    "    contourThickness=[4, 4],\n",
    "    dilationRadius=[3, 3],\n",
    "    colormap=red + green + blue,\n",
    ")\n",
    "gui.multi_image_display2D([contour_overlaid_image], figure_size=(6, 3));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparing two segmentations\n",
    "\n",
    "In this section we show how to create a binary image illustrating all the locations where two segmentations differ. This is a trivial one liner in SimpleITK.\n",
    "\n",
    "The following cell modifies our original coronal segmentation by dilating the body region in the top half of the image and dilating the lung region in the bottom half of the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "binary_dilate_filter = sitk.BinaryDilateImageFilter()\n",
    "binary_dilate_filter.SetKernelRadius(2)\n",
    "mid_size = int(coronal_segmentation_isotropic.GetHeight() / 2)\n",
    "\n",
    "# Over-segment the body region on the upper image region.\n",
    "binary_dilate_filter.SetForegroundValue(utilities.popi_body_label)\n",
    "top_segmentation = binary_dilate_filter.Execute(\n",
    "    coronal_segmentation_isotropic[:, 0:mid_size]\n",
    ")\n",
    "# Over-segment the lung region on the lower image region.\n",
    "binary_dilate_filter.SetForegroundValue(utilities.popi_lung_label)\n",
    "bottom_segmentation = binary_dilate_filter.Execute(\n",
    "    coronal_segmentation_isotropic[:, mid_size:]\n",
    ")\n",
    "\n",
    "modified_segmentation = sitk.Tile(top_segmentation, bottom_segmentation, (1, 2))\n",
    "modified_segmentation.CopyInformation(coronal_segmentation_isotropic)\n",
    "# Faux volume which allows us to visually compare the two segmentations by switching back and\n",
    "# forth between them.\n",
    "gui.MultiImageDisplay(\n",
    "    image_list=[sitk.JoinSeries(coronal_segmentation_isotropic, modified_segmentation)],\n",
    "    figure_size=(6, 3),\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see where the two segmentations differ, we directly compare them. If we don't want to waste ink, we can invert the result so that black pixels are the foreground and white the background. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff_image = coronal_segmentation_isotropic != modified_segmentation\n",
    "sitk.WriteImage(\n",
    "    diff_image * 255, os.path.join(OUTPUT_DIR, \"segmentation_differences.jpg\")\n",
    ")\n",
    "sitk.WriteImage(\n",
    "    (diff_image != 1) * 255,\n",
    "    os.path.join(OUTPUT_DIR, \"segmentation_differences_inverted.jpg\"),\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advanced visualization\n",
    "\n",
    "For more advanced visualization options we need to turn to other packages beyond SimpleITK. We highly recommend the [itkwidgets package](https://github.com/InsightSoftwareConsortium/itkwidgets) for various 3D visualizations. If you need finer control over the visualization there's always the visualization toolkit ([VTK](https://github.com/Kitware/VTK)). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itkwidgets\n",
    "\n",
    "itkwidgets.compare(\n",
    "    img,\n",
    "    segmentation == utilities.popi_lung_label,\n",
    "    cmap=itkwidgets.cm.BrBG,\n",
    "    annotations=False,\n",
    "    ui_collapsed=True,\n",
    "    link_cmap=False,\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}