09_segmentation_evaluation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\">Segmentation Evaluation</h1>\n",
    "\n",
    "**Summary:**\n",
    "\n",
    "1. SimpleITK supports two ways of combining expert segmentations to obtain a reference segmentation.\n",
    "2. A variety of criteria used for evaluating a segmentation result are readily available or implemented in SimpleITK.\n",
    "\n",
    "<u>Reference Segmentation</u>\n",
    "\n",
    "Evaluating segmentation algorithms is most often done using reference data to which you compare your results. In the medical domain reference data is commonly obtained via manual segmentation by an expert (don't forget to thank your clinical colleagues for their hard work). When you are resource limited, the reference data may be defined by a single expert. This is less than ideal. When multiple experts provide you with their input then you can potentially combine them to obtain reference data that is closer to the ever elusive \"ground truth\". In this notebook we show two approaches to combining input from multiple observers, majority vote and the Simultaneous Truth and Performance Level\n",
    "Estimation [(STAPLE)](https://www.ncbi.nlm.nih.gov/pubmed/15250643) algorithm.\n",
    "\n",
    "<u>Segmentation Evaluation</u>\n",
    "\n",
    "Once we have a reference, we compare the algorithm's performance using multiple criteria, as usually there is no single evaluation measure that conveys all of the relevant information. In this notebook we illustrate the use of the following evaluation criteria:\n",
    "* Overlap measures:\n",
    "  * Jaccard and Dice coefficients \n",
    "  * false negative and false positive errors\n",
    "* Surface distance measures:\n",
    "  * Hausdorff distance (symmetric)\n",
    "  * mean, median, max and standard deviation between surfaces\n",
    "* Volume measures:\n",
    "  * volume similarity $ \\frac{2*(v1-v2)}{v1+v2}$\n",
    "\n",
    "The relevant criteria are task dependent, so you need to ask yourself whether you are interested in detecting spurious errors or not (mean or max surface distance), whether over/under segmentation should be differentiated (volume similarity and Dice or just Dice), and what is the ratio between acceptable errors and the size of the segmented object (Dice coefficient may be too sensitive to small errors when the segmented object is small and not sensitive enough to large errors when the segmented object is large).\n",
    "\n",
    "In the context of segmentation challenges, algorithm rankings are often based on a weighted combination of these criteria. These ranking schemes are not necessarily robust, as discussed in \"[Why rankings of biomedical image analysis competitions should be interpreted with care](https://www.nature.com/articles/s41467-018-07619-7)\", L. Maier-Hein et al.\n",
    "\n",
    "The data we use in the notebook is a set of manually segmented liver tumors from a single clinical CT scan. The relevant publication is: T. Popa et al., \"Tumor Volume Measurement and Volume Measurement Comparison Plug-ins for VolView Using ITK\", SPIE Medical Imaging: Visualization, Image-Guided Procedures, and Display, 2006.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import SimpleITK as sitk\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from downloaddata import fetch_data as fdata\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import gui\n",
    "\n",
    "from ipywidgets import interact, fixed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utility method for display"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def display_with_overlay(\n",
    "    segmentation_number, slice_number, image, segs, window_min, window_max\n",
    "):\n",
    "    \"\"\"\n",
    "    Display a CT slice with segmented contours overlaid onto it. The contours are the edges of\n",
    "    the labeled regions.\n",
    "    \"\"\"\n",
    "    img = image[:, :, slice_number]\n",
    "    msk = segs[segmentation_number][:, :, slice_number]\n",
    "    overlay_img = sitk.LabelMapContourOverlay(\n",
    "        sitk.Cast(msk, sitk.sitkLabelUInt8),\n",
    "        sitk.Cast(\n",
    "            sitk.IntensityWindowing(\n",
    "                img, windowMinimum=window_min, windowMaximum=window_max\n",
    "            ),\n",
    "            sitk.sitkUInt8,\n",
    "        ),\n",
    "        opacity=1,\n",
    "        contourThickness=[2, 2],\n",
    "    )\n",
    "    # We assume the original slice is isotropic, otherwise the display would be distorted\n",
    "    plt.imshow(sitk.GetArrayViewFromImage(overlay_img))\n",
    "    plt.axis(\"off\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fetch the data\n",
    "\n",
    "Retrieve a single CT scan and three manual delineations of a liver tumor. Visual inspection of the data highlights the variability between experts. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image = sitk.ReadImage(fdata(\"liverTumorSegmentations/Patient01Homo.mha\"))\n",
    "segmentation_file_names = [\n",
    "    \"liverTumorSegmentations/Patient01Homo_Rad01.mha\",\n",
    "    \"liverTumorSegmentations/Patient01Homo_Rad02.mha\",\n",
    "    \"liverTumorSegmentations/Patient01Homo_Rad03.mha\",\n",
    "]\n",
    "\n",
    "segmentations = [\n",
    "    sitk.ReadImage(fdata(file_name), sitk.sitkUInt8)\n",
    "    for file_name in segmentation_file_names\n",
    "]\n",
    "\n",
    "interact(\n",
    "    display_with_overlay,\n",
    "    segmentation_number=(0, len(segmentations) - 1),\n",
    "    slice_number=(0, image.GetSize()[2] - 1),\n",
    "    image=fixed(image),\n",
    "    segs=fixed(segmentations),\n",
    "    window_min=fixed(-1024),\n",
    "    window_max=fixed(976),\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Derive a reference\n",
    "\n",
    "There are a variety of ways to derive a reference segmentation from multiple expert inputs (\"[A comparison of ground truth estimation methods](https://www.ncbi.nlm.nih.gov/pubmed/20033494)\", A. M. Biancardi, A. C. Jirapatnakul, A. P. Reeves).\n",
    "\n",
    "Two methods that are available in SimpleITK are [majority vote](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1LabelVotingImageFilter.html) and the STAPLE algorithm ([single label](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1STAPLEImageFilter.html) or [multi label](https://simpleitk.org/doxygen/latest/html/classitk_1_1simple_1_1MultiLabelSTAPLEImageFilter.html))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the STAPLE algorithm to obtain the reference segmentation. This implementation of the original algorithm\n",
    "# combines a single label from multiple segmentations, the label is user specified. The result of the\n",
    "# filter is the voxel's probability of belonging to the foreground. We then have to threshold the result to obtain\n",
    "# a reference binary segmentation.\n",
    "foregroundValue = 1\n",
    "threshold = 0.95\n",
    "reference_segmentation_STAPLE_probabilities = sitk.STAPLE(\n",
    "    segmentations, foregroundValue\n",
    ")\n",
    "# We use the overloaded operator to perform thresholding, another option is to use the BinaryThreshold function.\n",
    "reference_segmentation = reference_segmentation_STAPLE_probabilities > threshold\n",
    "\n",
    "manual_plus_staple = list(segmentations)\n",
    "# Append the reference segmentation to the list of manual segmentations\n",
    "manual_plus_staple.append(reference_segmentation)\n",
    "\n",
    "interact(\n",
    "    display_with_overlay,\n",
    "    segmentation_number=(0, len(manual_plus_staple) - 1),\n",
    "    slice_number=(0, image.GetSize()[2] - 1),\n",
    "    image=fixed(image),\n",
    "    segs=fixed(manual_plus_staple),\n",
    "    window_min=fixed(-1024),\n",
    "    window_max=fixed(976),\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate segmentations using the reference\n",
    "\n",
    "Once we derive a reference from our experts input we can compare segmentation results to it.\n",
    "\n",
    "Note that in this notebook we compare the expert segmentations to the reference derived from them. This is not relevant for algorithm evaluation, but it can potentially be used to rank your experts.\n",
    "\n",
    "In this specific implementation we take advantage of the fact that we have a binary segmentation with 1 for foreground and 0 for background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from enum import Enum\n",
    "\n",
    "\n",
    "# Use enumerations to represent the various evaluation measures\n",
    "class OverlapMeasures(Enum):\n",
    "    jaccard, dice, volume_similarity, false_negative, false_positive = range(5)\n",
    "\n",
    "\n",
    "class SurfaceDistanceMeasures(Enum):\n",
    "    (\n",
    "        hausdorff_distance,\n",
    "        mean_surface_distance,\n",
    "        median_surface_distance,\n",
    "        std_surface_distance,\n",
    "        max_surface_distance,\n",
    "    ) = range(5)\n",
    "\n",
    "\n",
    "# Empty numpy arrays to hold the results\n",
    "overlap_results = np.zeros(\n",
    "    (len(segmentations), len(OverlapMeasures.__members__.items()))\n",
    ")\n",
    "surface_distance_results = np.zeros(\n",
    "    (len(segmentations), len(SurfaceDistanceMeasures.__members__.items()))\n",
    ")\n",
    "\n",
    "# Compute the evaluation criteria\n",
    "\n",
    "# Note that for the overlap measures filter, because we are dealing with a single label we\n",
    "# use the combined, all labels, evaluation measures without passing a specific label to the methods.\n",
    "overlap_measures_filter = sitk.LabelOverlapMeasuresImageFilter()\n",
    "\n",
    "hausdorff_distance_filter = sitk.HausdorffDistanceImageFilter()\n",
    "\n",
    "# Use the absolute values of the distance map to compute the surface distances (distance map sign, outside or inside\n",
    "# relationship, is irrelevant)\n",
    "label = 1\n",
    "reference_distance_map = sitk.Abs(\n",
    "    sitk.SignedMaurerDistanceMap(\n",
    "        reference_segmentation, squaredDistance=False, useImageSpacing=True\n",
    "    )\n",
    ")\n",
    "reference_surface = sitk.LabelContour(reference_segmentation)\n",
    "\n",
    "statistics_image_filter = sitk.StatisticsImageFilter()\n",
    "# Get the number of pixels in the reference surface by counting all pixels that are 1.\n",
    "statistics_image_filter.Execute(reference_surface)\n",
    "num_reference_surface_pixels = int(statistics_image_filter.GetSum())\n",
    "\n",
    "for i, seg in enumerate(segmentations):\n",
    "    # Overlap measures\n",
    "    overlap_measures_filter.Execute(reference_segmentation, seg)\n",
    "    overlap_results[i, OverlapMeasures.jaccard.value] = (\n",
    "        overlap_measures_filter.GetJaccardCoefficient()\n",
    "    )\n",
    "    overlap_results[i, OverlapMeasures.dice.value] = (\n",
    "        overlap_measures_filter.GetDiceCoefficient()\n",
    "    )\n",
    "    overlap_results[i, OverlapMeasures.volume_similarity.value] = (\n",
    "        overlap_measures_filter.GetVolumeSimilarity()\n",
    "    )\n",
    "    overlap_results[i, OverlapMeasures.false_negative.value] = (\n",
    "        overlap_measures_filter.GetFalseNegativeError()\n",
    "    )\n",
    "    overlap_results[i, OverlapMeasures.false_positive.value] = (\n",
    "        overlap_measures_filter.GetFalsePositiveError()\n",
    "    )\n",
    "    # Hausdorff distance\n",
    "    hausdorff_distance_filter.Execute(reference_segmentation, seg)\n",
    "\n",
    "    surface_distance_results[i, SurfaceDistanceMeasures.hausdorff_distance.value] = (\n",
    "        hausdorff_distance_filter.GetHausdorffDistance()\n",
    "    )\n",
    "    # Symmetric surface distance measures\n",
    "    segmented_distance_map = sitk.Abs(\n",
    "        sitk.SignedMaurerDistanceMap(seg, squaredDistance=False, useImageSpacing=True)\n",
    "    )\n",
    "    segmented_surface = sitk.LabelContour(seg)\n",
    "\n",
    "    # Multiply the binary surface segmentations with the distance maps. The resulting distance\n",
    "    # maps contain non-zero values only on the surface (they can also contain zero on the surface)\n",
    "    seg2ref_distance_map = reference_distance_map * sitk.Cast(\n",
    "        segmented_surface, sitk.sitkFloat32\n",
    "    )\n",
    "    ref2seg_distance_map = segmented_distance_map * sitk.Cast(\n",
    "        reference_surface, sitk.sitkFloat32\n",
    "    )\n",
    "\n",
    "    # Get the number of pixels in the reference surface by counting all pixels that are 1.\n",
    "    statistics_image_filter.Execute(segmented_surface)\n",
    "    num_segmented_surface_pixels = int(statistics_image_filter.GetSum())\n",
    "\n",
    "    # Get all non-zero distances and then add zero distances if required.\n",
    "    seg2ref_distance_map_arr = sitk.GetArrayViewFromImage(seg2ref_distance_map)\n",
    "    seg2ref_distances = list(seg2ref_distance_map_arr[seg2ref_distance_map_arr != 0])\n",
    "    seg2ref_distances = seg2ref_distances + list(\n",
    "        np.zeros(num_segmented_surface_pixels - len(seg2ref_distances))\n",
    "    )\n",
    "    ref2seg_distance_map_arr = sitk.GetArrayViewFromImage(ref2seg_distance_map)\n",
    "    ref2seg_distances = list(ref2seg_distance_map_arr[ref2seg_distance_map_arr != 0])\n",
    "    ref2seg_distances = ref2seg_distances + list(\n",
    "        np.zeros(num_reference_surface_pixels - len(ref2seg_distances))\n",
    "    )\n",
    "\n",
    "    all_surface_distances = seg2ref_distances + ref2seg_distances\n",
    "\n",
    "    # The maximum of the symmetric surface distances is the Hausdorff distance between the surfaces. In\n",
    "    # general, it is not equal to the Hausdorff distance between all voxel/pixel points of the two\n",
    "    # segmentations, though in our case it is. More on this below.\n",
    "    surface_distance_results[i, SurfaceDistanceMeasures.mean_surface_distance.value] = (\n",
    "        np.mean(all_surface_distances)\n",
    "    )\n",
    "    surface_distance_results[\n",
    "        i, SurfaceDistanceMeasures.median_surface_distance.value\n",
    "    ] = np.median(all_surface_distances)\n",
    "    surface_distance_results[i, SurfaceDistanceMeasures.std_surface_distance.value] = (\n",
    "        np.std(all_surface_distances)\n",
    "    )\n",
    "    surface_distance_results[i, SurfaceDistanceMeasures.max_surface_distance.value] = (\n",
    "        np.max(all_surface_distances)\n",
    "    )\n",
    "\n",
    "# Print the matrices\n",
    "np.set_printoptions(precision=3)\n",
    "print(overlap_results)\n",
    "print(surface_distance_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Improved output\n",
    "\n",
    "Using the [pandas](http://pandas.pydata.org/) package we can easily produce high quality output. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from IPython.display import display, HTML\n",
    "\n",
    "# Graft our results matrix into pandas data frames\n",
    "overlap_results_df = pd.DataFrame(\n",
    "    data=overlap_results,\n",
    "    index=list(range(len(segmentations))),\n",
    "    columns=[name for name, _ in OverlapMeasures.__members__.items()],\n",
    ")\n",
    "surface_distance_results_df = pd.DataFrame(\n",
    "    data=surface_distance_results,\n",
    "    index=list(range(len(segmentations))),\n",
    "    columns=[name for name, _ in SurfaceDistanceMeasures.__members__.items()],\n",
    ")\n",
    "\n",
    "# Display the data as HTML tables and graphs\n",
    "display(HTML(overlap_results_df.to_html(float_format=lambda x: \"%.3f\" % x)))\n",
    "display(HTML(surface_distance_results_df.to_html(float_format=lambda x: \"%.3f\" % x)))\n",
    "overlap_results_df.plot(kind=\"bar\").legend(bbox_to_anchor=(1.6, 0.9))\n",
    "surface_distance_results_df.plot(kind=\"bar\").legend(bbox_to_anchor=(1.6, 0.9))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also export the data as a table for your LaTeX manuscript using the [to_latex](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.to_latex.html) function.\n",
    "<b>Note</b>: You will need to add the \\usepackage{booktabs} to your LaTeX document's preamble. \n",
    "\n",
    "To create the minimal LaTeX document which will allow you to see the difference between the tables below, copy paste:\n",
    "\n",
    "\\documentclass{article}\n",
    "\n",
    "\\usepackage{booktabs}\n",
    "\n",
    "\\begin{document}\n",
    "\n",
    "paste the tables here\n",
    "\n",
    "\\end{document}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The formatting of the table using the default settings is less than ideal\n",
    "print(overlap_results_df.to_latex())\n",
    "\n",
    "# We can improve on this by specifying the table's column format and the float format\n",
    "print(\n",
    "    overlap_results_df.to_latex(\n",
    "        column_format=\"ccccccc\", float_format=lambda x: \"%.3f\" % x\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visual Diff\n",
    "\n",
    "It is always nice to have a figure with a visual display of the difference between the segmentation and ground truth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "simpleitk_error_allowed": "Exception thrown in SimpleITK Show:"
   },
   "outputs": [],
   "source": [
    "# Use the first segmentation\n",
    "segmentation = segmentations[0]\n",
    "\n",
    "# Save ink, the differences will be in black and background is white\n",
    "segmentation_diff = (segmentation == reference_segmentation) * 255\n",
    "\n",
    "# Flatten for 2D presentation, create a montage from the volume\n",
    "num_slices = segmentation_diff.GetDepth()\n",
    "tile_w = int(np.sqrt(num_slices))\n",
    "tile_h = int(np.ceil(num_slices / tile_w))\n",
    "default_background_color = 255\n",
    "tile_image = sitk.Tile(\n",
    "    [segmentation_diff[:, :, i] for i in range(num_slices)],\n",
    "    (tile_w, tile_h),\n",
    "    default_background_color,\n",
    ")\n",
    "sitk.Show(tile_image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"10_results_visualization.ipynb\"><h2 align=right>Next &raquo;</h2></a>"
   ]
  }
 ],
 "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
}