WIP: ENH: Add watershed with distance map example.

Add watershed with distance map example.

Resolves #48.

src/Segmentation/Watersheds/CMakeLists.txt

@@ -1,4 +1,12 @@
 
+add_example(SegmentWithWatershedAndDistanceMap)
+compare_to_baseline(EXAMPLE_NAME SegmentWithWatershedAndDistanceMap
+  TEST_NAME SegmentWithWatershedAndDistanceMapTest01BaselineComparison
+  BASELINE_PREFIX SegmentWithWatershedAndDistanceMapTest01Baseline
+  DEPENDS SegmentWithWatershedAndDistanceMapTest01
+  TEST_IMAGE SegmentWithWatershedAndDistanceMapTest01.png
+  )
+
 add_example(SegmentWithWatershedImageFilter)
 
 foreach( i RANGE 1 5 )

src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/BubbleVolume.csv.sha512

@@ -0,0 +1 @@
+0d5403fda7b3c94d3229bf5e42b523a4a5bb2cd84a45973011520daf5b3001192e74ff3d290b43c0feb0ebd84c17a328f4646a7887408c1f63e52927676f2360

src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/CMakeLists.txt

@@ -0,0 +1,59 @@
+cmake_minimum_required(VERSION 3.10.2)
+
+project(SegmentWithWatershedAndDistanceMap)
+
+find_package(ITK REQUIRED)
+include(${ITK_USE_FILE})
+
+
+add_executable(SegmentWithWatershedAndDistanceMap Code.cxx)
+target_link_libraries(SegmentWithWatershedAndDistanceMap ${ITK_LIBRARIES})
+
+install(TARGETS SegmentWithWatershedAndDistanceMap
+  DESTINATION bin/ITKExamples/Segmentation/Watersheds
+  COMPONENT Runtime
+  )
+
+install(FILES Code.cxx CMakeLists.txt Code.py
+  DESTINATION share/ITKExamples/Code/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/
+  COMPONENT Code
+  )
+
+
+enable_testing()
+set(input_image ${CMAKE_CURRENT_BINARY_DIR}/PlateauBorder.tif)
+set(binarizingRadius01 2)
+set(majorityThreshold01 2)
+set(watershedThreshold01 0.01)
+set(watershedLevel01 0.5)
+set(cleaningStructuringElementRadius01 3)
+
+add_test(NAME SegmentWithWatershedAndDistanceMapTest01
+  COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/SegmentWithWatershedAndDistanceMap
+    ${input_image}
+    ReversedInputImageTest01.tif
+    DistanceMapImageTest01.tif
+    WatershedImageTest01.tif
+    SegmentWithWatershedAndDistanceMapTest01.tif
+    ${binarizingRadius01}
+    ${majorityThreshold01}
+    ${watershedThreshold01}
+    ${watershedLevel01}
+    ${cleaningStructuringElementRadius01}
+  )
+
+if(ITK_WRAP_PYTHON)
+  add_test(NAME SegmentWithWatershedAndDistanceMapTest01Python
+    COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/Code.py
+      ${input_image}
+      ReversedInputImageTest01.tif
+      DistanceMapImageTest01.tif
+      WatershedImageTest01.tif
+      SegmentWithWatershedAndDistanceMapTest01Python.tif
+      ${binarizingRadius01}
+      ${majorityThreshold01}
+      ${watershedThreshold01}
+      ${watershedLevel01}
+      ${cleaningStructuringElementRadius01}
+    )
+endif()

src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/Code.cxx

@@ -0,0 +1,177 @@
+/*=========================================================================
+ *
+ *  Copyright NumFOCUS
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *         http://www.apache.org/licenses/LICENSE-2.0.txt
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ *
+ *=========================================================================*/
+
+#include "itkBinaryBallStructuringElement.h"
+#include "itkBinaryMorphologicalOpeningImageFilter.h"
+#include "itkImageFileReader.h"
+#include "itkImageFileWriter.h"
+#include "itkScalarToRGBColormapImageFilter.h"
+#include "itkSignedMaurerDistanceMapImageFilter.h"
+#include "itkVotingBinaryIterativeHoleFillingImageFilter.h"
+#include "itkWatershedImageFilter.h"
+
+// Run with:
+// ./SegmentWithWatershedAndDistanceMap <inputImageFile>
+// <reversedInputImageFile> <distanceMapOutputImageFile>
+// <watershedOutputFileName> <segmentationResultOutputImageFile>
+// binarizingRadius majorityThreshold watershedThreshold watershedLevel
+// cleaningStructuringElementRadius
+// e.g.
+// ./SegmentWithWatershedAndDistanceMap PlateauBorder.tif
+// reversedInputImage.tif distanceMap.tif watershed.tif segmentationResult.tif
+// 2 2 0.01 0.5 3
+// (A rule of thumb is to set the threshold to be about 1 / 100 of the level.)
+
+int
+main(int argc, char * argv[])
+{
+  if (argc < 10)
+  {
+    std::cerr << "Missing parameters." << std::endl;
+    std::cerr << "Usage: " << argv[0] << " inputImageFile"
+              << " reversedInputImageFile"
+              << " distanceMapOutputImageFile"
+              << " watershedOutputFileName"
+              << " segmentationResultOutputImageFile"
+              << " binarizingRadius"
+              << " majorityThreshold"
+              << " watershedThreshold"
+              << " watershedLevel"
+              << " cleaningStructuringElementRadius" << std::endl;
+    return EXIT_FAILURE;
+  }
+
+  constexpr unsigned int Dimension = 3;
+
+  using UnsignedCharPixelType = unsigned char;
+  using FloatPixelType = float;
+
+  using InputImageType = itk::Image<UnsignedCharPixelType, Dimension>;
+  using FloatImageType = itk::Image<FloatPixelType, Dimension>;
+  using RGBPixelType = itk::RGBPixel<UnsignedCharPixelType>;
+  using RGBImageType = itk::Image<RGBPixelType, Dimension>;
+  using LabeledImageType = itk::Image<itk::IdentifierType, Dimension>;
+
+
+  using FileReaderType = itk::ImageFileReader<InputImageType>;
+  FileReaderType::Pointer reader = FileReaderType::New();
+  reader->SetFileName(argv[1]);
+  reader->Update();
+
+
+  // Create bubble image: get a binarized version of the input image
+  using VotingBinaryIterativeHoleFillingImageFilterType =
+    itk::VotingBinaryIterativeHoleFillingImageFilter<InputImageType>;
+  VotingBinaryIterativeHoleFillingImageFilterType::Pointer votingBinaryHoleFillingImageFilter =
+    VotingBinaryIterativeHoleFillingImageFilterType::New();
+  votingBinaryHoleFillingImageFilter->SetInput(reader->GetOutput());
+
+  const unsigned int binarizingRadius = std::stoi(argv[6]);
+
+  InputImageType::SizeType indexRadius;
+  indexRadius.Fill(binarizingRadius);
+
+  votingBinaryHoleFillingImageFilter->SetRadius(indexRadius);
+
+  votingBinaryHoleFillingImageFilter->SetBackgroundValue(0);
+  votingBinaryHoleFillingImageFilter->SetForegroundValue(255);
+
+  const unsigned int majorityThreshold = std::stoi(argv[7]);
+  votingBinaryHoleFillingImageFilter->SetMajorityThreshold(majorityThreshold);
+
+  votingBinaryHoleFillingImageFilter->Update();
+
+  using FileWriterType = itk::ImageFileWriter<InputImageType>;
+  FileWriterType::Pointer reversedImageWriter = FileWriterType::New();
+  reversedImageWriter->SetFileName(argv[2]);
+  reversedImageWriter->SetInput(votingBinaryHoleFillingImageFilter->GetOutput());
+  reversedImageWriter->Update();
+
+
+  // Get the distance map of the input image
+  using SignedMaurerDistanceMapImageFilterType =
+    itk::SignedMaurerDistanceMapImageFilter<InputImageType, FloatImageType>;
+  SignedMaurerDistanceMapImageFilterType::Pointer distanceMapImageFilter =
+    SignedMaurerDistanceMapImageFilterType::New();
+  distanceMapImageFilter->SetInput(votingBinaryHoleFillingImageFilter->GetOutput());
+
+  distanceMapImageFilter->SetInsideIsPositive(false);
+  distanceMapImageFilter->Update();
+
+
+  using DistanceMapFileWriterType = itk::ImageFileWriter<FloatImageType>;
+  DistanceMapFileWriterType::Pointer distanceMapWriter = DistanceMapFileWriterType::New();
+  distanceMapWriter->SetFileName(argv[3]);
+  distanceMapWriter->SetInput(distanceMapImageFilter->GetOutput());
+  distanceMapWriter->Update();
+
+
+  // Apply the watershed segmentation
+  using WatershedFilterType = itk::WatershedImageFilter<FloatImageType>;
+  WatershedFilterType::Pointer watershed = WatershedFilterType::New();
+
+  const float watershedThreshold = std::stod(argv[8]);
+  const float watershedLevel = std::stod(argv[9]);
+
+  watershed->SetThreshold(watershedThreshold);
+  watershed->SetLevel(watershedLevel);
+
+  watershed->SetInput(distanceMapImageFilter->GetOutput());
+  watershed->Update();
+
+
+  using RGBFilterType = itk::ScalarToRGBColormapImageFilter<LabeledImageType, RGBImageType>;
+  RGBFilterType::Pointer colormapImageFilter = RGBFilterType::New();
+  colormapImageFilter->SetColormap(RGBFilterType::Jet);
+  colormapImageFilter->SetInput(watershed->GetOutput());
+  colormapImageFilter->Update();
+
+  using WatershedFileWriterType = itk::ImageFileWriter<RGBImageType>;
+  WatershedFileWriterType::Pointer watershedWriter = WatershedFileWriterType::New();
+  watershedWriter->SetFileName(argv[4]);
+  watershedWriter->SetInput(colormapImageFilter->GetOutput());
+  watershedWriter->Update();
+
+
+  // Clean the segmentation image: remove small objects by performing an
+  // opening morphological operation
+  using StructuringElementType =
+    itk::BinaryBallStructuringElement<LabeledImageType::PixelType, LabeledImageType::ImageDimension>;
+  StructuringElementType structuringElement;
+
+  const unsigned int cleaningStructuringElementRadius = std::stoi(argv[10]);
+  structuringElement.SetRadius(cleaningStructuringElementRadius);
+  structuringElement.CreateStructuringElement();
+
+  using BinaryMorphologicalOpeningImageFilterType =
+    itk::BinaryMorphologicalOpeningImageFilter<LabeledImageType, LabeledImageType, StructuringElementType>;
+  BinaryMorphologicalOpeningImageFilterType::Pointer openingFilter = BinaryMorphologicalOpeningImageFilterType::New();
+  openingFilter->SetInput(watershed->GetOutput());
+  openingFilter->SetKernel(structuringElement);
+  openingFilter->Update();
+
+
+  using SegmentationFileWriterType = itk::ImageFileWriter<RGBImageType>;
+  SegmentationFileWriterType::Pointer segmentationWriter = SegmentationFileWriterType::New();
+  segmentationWriter->SetFileName(argv[5]);
+  segmentationWriter->SetInput(colormapImageFilter->GetOutput());
+  segmentationWriter->Update();
+
+
+  return EXIT_SUCCESS;
+}

src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/Code.ipynb

@@ -0,0 +1,109 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "id": "9e8cfdc6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sys\n",
+    "\n",
+    "import itk\n",
+    "import numpy as np\n",
+    "\n",
+    "from itkwidgets import view, cm\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "id": "057c9f2f",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "962ea74edef74bc4b0b026a83f3dde3a",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "Viewer(cmap=['bone_Matlab'], geometries=[], gradient_opacity=0.22, point_sets=[], rendered_image=<itk.itkImage…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "stack_image = itk.imread('./PlateauBorder.tif')\n",
+    "view(stack_image, cmap=[cm.bone])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "id": "ab956229",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "c0fc5e027fe149d4b4bc866b054c6e38",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "Viewer(cmap=['bone_Matlab'], geometries=[], gradient_opacity=0.22, point_sets=[], rendered_image=<itk.itkImage…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "binarizing_radius = 2\n",
+    "majority_threshold = 2\n",
+    "\n",
+    "bubble_image = itk.voting_binary_hole_filling_image_filter(stack_image,\n",
+    "                                                           radius=binarizing_radius,\n",
+    "                                                           majority_threshold=majority_threshold,\n",
+    "                                                           background_value=0,\n",
+    "                                                           foreground_value=255)\n",
+    "bubble_image = itk.invert_intensity_image_filter(bubble_image)\n",
+    "view(bubble_image, cmap=[cm.bone,])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bfc02551",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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.8.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}