diff --git a/src/Segmentation/Watersheds/CMakeLists.txt b/src/Segmentation/Watersheds/CMakeLists.txt index 688e8b868..4931d0abd 100644 --- a/src/Segmentation/Watersheds/CMakeLists.txt +++ b/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 ) diff --git a/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/BubbleVolume.csv.sha512 b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/BubbleVolume.csv.sha512 new file mode 100644 index 000000000..e1832505b --- /dev/null +++ b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/BubbleVolume.csv.sha512 @@ -0,0 +1 @@ +0d5403fda7b3c94d3229bf5e42b523a4a5bb2cd84a45973011520daf5b3001192e74ff3d290b43c0feb0ebd84c17a328f4646a7887408c1f63e52927676f2360 diff --git a/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/CMakeLists.txt b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/CMakeLists.txt new file mode 100644 index 000000000..305596274 --- /dev/null +++ b/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() diff --git a/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/Code.cxx b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/Code.cxx new file mode 100644 index 000000000..e058485ca --- /dev/null +++ b/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 +// +// +// 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; + using FloatImageType = itk::Image; + using RGBPixelType = itk::RGBPixel; + using RGBImageType = itk::Image; + using LabeledImageType = itk::Image; + + + using FileReaderType = itk::ImageFileReader; + 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; + 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; + 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; + SignedMaurerDistanceMapImageFilterType::Pointer distanceMapImageFilter = + SignedMaurerDistanceMapImageFilterType::New(); + distanceMapImageFilter->SetInput(votingBinaryHoleFillingImageFilter->GetOutput()); + + distanceMapImageFilter->SetInsideIsPositive(false); + distanceMapImageFilter->Update(); + + + using DistanceMapFileWriterType = itk::ImageFileWriter; + DistanceMapFileWriterType::Pointer distanceMapWriter = DistanceMapFileWriterType::New(); + distanceMapWriter->SetFileName(argv[3]); + distanceMapWriter->SetInput(distanceMapImageFilter->GetOutput()); + distanceMapWriter->Update(); + + + // Apply the watershed segmentation + using WatershedFilterType = itk::WatershedImageFilter; + 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; + RGBFilterType::Pointer colormapImageFilter = RGBFilterType::New(); + colormapImageFilter->SetColormap(RGBFilterType::RGBColormapFilterEnum::Jet); + colormapImageFilter->SetInput(watershed->GetOutput()); + colormapImageFilter->Update(); + + using WatershedFileWriterType = itk::ImageFileWriter; + 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; + StructuringElementType structuringElement; + + const unsigned int cleaningStructuringElementRadius = std::stoi(argv[10]); + structuringElement.SetRadius(cleaningStructuringElementRadius); + structuringElement.CreateStructuringElement(); + + using BinaryMorphologicalOpeningImageFilterType = + itk::BinaryMorphologicalOpeningImageFilter; + BinaryMorphologicalOpeningImageFilterType::Pointer openingFilter = BinaryMorphologicalOpeningImageFilterType::New(); + openingFilter->SetInput(watershed->GetOutput()); + openingFilter->SetKernel(structuringElement); + openingFilter->Update(); + + + using SegmentationFileWriterType = itk::ImageFileWriter; + SegmentationFileWriterType::Pointer segmentationWriter = SegmentationFileWriterType::New(); + segmentationWriter->SetFileName(argv[5]); + segmentationWriter->SetInput(colormapImageFilter->GetOutput()); + segmentationWriter->Update(); + + + return EXIT_SUCCESS; +} diff --git a/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/Code.ipynb b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/Code.ipynb new file mode 100644 index 000000000..7dcf139c7 --- /dev/null +++ b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/Code.ipynb @@ -0,0 +1,596 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "source": [ + "# Watershed with Distance Map\n", + "This example illustrates how to segment an image using the watershed method\n", + "and the signed Maurer distance map." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "import itk\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import pandas as pd\n" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "Define the types to be used" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "dimension = 3\n", + "\n", + "uchar_pixel_type = itk.UC\n", + "uchar_image_type = itk.Image[uchar_pixel_type, dimension]\n", + "\n", + "float_pixel_type = itk.F\n", + "float_image_type = itk.Image[float_pixel_type, dimension]\n", + "\n", + "rgb_pixel_type = itk.RGBPixel[uchar_pixel_type]\n", + "RGBImageType = itk.Image[rgb_pixel_type, dimension]" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "Display the input image and print its shape." + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%% md\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "stack_image = itk.imread(\"PlateauBorder.tif\")\n", + "print(stack_image.shape, stack_image.dtype)" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "# Display the input image projections\n", + "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 4))\n", + "for i, (cax, clabel) in enumerate(zip([ax1, ax2, ax3], [\"xy\", \"zy\", \"zx\"])):\n", + " cax.imshow(np.sum(stack_image, i).squeeze(), cmap=\"bone_r\")\n", + " cax.set_title(\"{} projection\".format(clabel))\n", + " cax.set_xlabel(clabel[0])\n", + " cax.set_ylabel(clabel[1])" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "## Create a bubble image\n", + "The bubble image is the reverse of the plateau border image: there cannot be\n", + "air where there is water." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "index_radius = itk.Size[dimension]()\n", + "index_radius.Fill(args.binarizing_radius)\n", + "\n", + "bubble_image = itk.voting_binary_iterative_hole_filling_image_filter(\n", + " stack_image,\n", + " radius=index_radius,\n", + " background_value=0,\n", + " foreground_value=255,\n", + " majority_threshold=args.majority_threshold,\n", + ")\n", + "\n", + "plt.imshow(bubble_image[5], cmap=\"bone\")\n", + "\n", + "# Write bubble image\n", + "itk.imwrite(bubble_image, \"ReversedInputImageTest01.tif\")" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "## Watershed on bubbles\n", + "Use the ITK watershed operation on the distance map of the input image." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "%%time\n", + "\n", + "bbl_cast_image = itk.cast_image_filter(\n", + " bubble_image,\n", + " ttype=(uchar_image_type, float_image_type),\n", + ")\n", + "\n", + "# Normalize the image to the [0, 255] range\n", + "bubble_image_preclamp = itk.multiply_image_filter(\n", + " bbl_cast_image,\n", + " constant=255.0,\n", + ")\n", + "bubble_image_clamp = itk.clamp_image_filter(\n", + " bubble_image_preclamp,\n", + " bounds=(0, 255),\n", + ")\n", + "\n", + "# Get the distance map of the input image\n", + "distance_map_image = itk.signed_maurer_distance_map_image_filter(\n", + " bubble_image_clamp,\n", + " inside_is_positive=False,\n", + ")\n", + "itk.imwrite(distance_map_image, args.distance_map_output_filename)\n", + "\n", + "# Apply the watershed segmentation\n", + "watershed_image = itk.watershed_image_filter(\n", + " distance_map_image,\n", + " threshold=args.watershed_threshold,\n", + " level=args.level,\n", + ")\n", + "\n", + "# Cast to unsigned char so that it can be written as a TIFF image\n", + "# WatershedImageFilter produces itk.ULL, but CastImageFilter does not wrap\n", + "# itk.ULL: see ITK issue 2551\n", + "# ws_cast_image = itk.cast_image_filter(watershed_image, ttype=(watershed_image.__class__, uchar_image_type))\n", + "# Workaround\n", + "rgb_ws_image = itk.scalar_to_rgb_colormap_image_filter(\n", + " watershed_image,\n", + " colormap=itk.ScalarToRGBColormapImageFilterEnums.RGBColormapFilter_Jet,\n", + ")\n", + "\n", + "# itk.imwrite(ws_cast_image, args.watershed_output_filename)\n", + "itk.imwrite(rgb_ws_image, args.watershed_output_filename)" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "## Cut-through\n", + "Show the values at a slice in the middle as a way to get a feeling for what the\n", + "watershed and distance map did." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "mid_slice = watershed_image.shape[0] // 2\n", + "ws_vol_arr = itk.array_view_from_image(watershed_image)\n", + "\n", + "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 7))\n", + "ax1.imshow(bubble_image[mid_slice], cmap=\"bone\")\n", + "ax1.set_title(\"Bubble image\")\n", + "m_val = np.abs(dmap_vol[mid_slice]).std()\n", + "ax2.imshow(dmap_vol[mid_slice], cmap=\"RdBu\", vmin=-m_val, vmax=m_val)\n", + "ax2.set_title(\"Distance image\\nmin: {:2.2f}; max: {:2.2f}; mean: {:2.2f}\".format(\n", + " distance_map_image[mid_slice].min(),\n", + " distance_map_image[mid_slice].max(),\n", + " distance_map_image[mid_slice].mean())\n", + ")\n", + "ax3.imshow(ws_vol[mid_slice], cmap=\"nipy_spectral\")\n", + "ax3.set_title(\"Watershed\\nLabels found: {}\".format(\n", + " len(np.unique(ws_vol_arr[ws_vol_arr > 0])))\n", + ")\n", + "\n", + "fig.savefig(\"Pipeline_images.png\")" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "\n", + "# Show segmentation result projections\n", + "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 4))\n", + "for i, (cax, clabel) in enumerate(zip([ax1, ax2, ax3], [\"xy\", \"zy\", \"zx\"])):\n", + " cax.imshow(np.max(ws_vol_arr, i).squeeze(), cmap=\"nipy_spectral\")\n", + " cax.set_title(\"{} projection\".format(clabel))\n", + " cax.set_xlabel(clabel[0])\n", + " cax.set_ylabel(clabel[1])\n", + "\n", + "fig.savefig(\"Segmentation_projections.png\")" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "## Clean the segmentation and relabel\n", + "Clean the segmentation image: remove small objects by performing an\n", + "opening morphological operation and relabel in order." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "bubble_label_image = np.zeros(watershed_image.shape).astype(np.uint16)\n", + "new_idx = 1\n", + "bubble_ids = [\n", + " (idx, np.sum(watershed_image[watershed_image==idx] > 0))\n", + " for idx in np.unique(watershed_image[watershed_image > 0])\n", + "]\n", + "\n", + "dimension = len(np.shape(bubble_image))\n", + "structuring_element_type = itk.FlatStructuringElement[dimension]\n", + "# Bubbles are round\n", + "structuring_element = structuring_element_type.Ball(\n", + " args.cleaning_structuring_element_radius\n", + ")\n", + "\n", + "from tqdm import tqdm\n", + "\n", + "# Count the kept bubbles in bubble label image\n", + "for old_idx, vol in tqdm(sorted(bubble_ids, key = lambda x: x[1])):\n", + " if 40000 < vol < 400000:\n", + " old_img = watershed_image==old_idx\n", + " cleaned_img = itk.binary_morphological_opening_image_filter(\n", + " old_img,\n", + " kernel=structuring_element,\n", + " )\n", + " bubble_label_image[old_img] = new_idx\n", + " new_idx += 1\n", + "\n", + "print(\"Total bubbles kept: {}/{}\".format(new_idx, len(bubble_ids)))\n", + "\n", + "# itk.imwrite(segmented_clean_image, \"Clean_segmentation.tif\")" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 4))\n", + "for i, (cax, clabel) in enumerate(zip([ax1, ax2, ax3], [\"xy\", \"zy\", \"zx\"])):\n", + " cax.imshow(\n", + " np.max(bubble_label_image, i),\n", + " cmap=\"jet\",\n", + " vmin=0,\n", + " vmax=new_idx,\n", + " )\n", + " cax.set_title(\"{} projection\".format(clabel))\n", + " cax.set_xlabel(clabel[0])\n", + " cax.set_ylabel(clabel[1])\n", + "\n", + "fig.savefig(\"Clean_segmentation_projections.png\")" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "## Show 3D rendering" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "from mpl_toolkits.mplot3d.art3d import Poly3DCollection\n", + "from skimage import measure\n", + "from tqdm import tqdm\n", + "\n", + "# Show 3D rendering\n", + "def show_3d_mesh(image, thresholds):\n", + " p = image[::-1].swapaxes(1, 2)\n", + " cmap = plt.cm.get_cmap(\"nipy_spectral_r\")\n", + " _fig = plt.figure(figsize=(10, 10))\n", + " ax = _fig.add_subplot(111, projection=\"3d\")\n", + " for _i, c_threshold in tqdm(list(enumerate(thresholds))):\n", + " verts, faces, _, _ = measure.marching_cubes(p==c_threshold, level=0)\n", + " mesh = Poly3DCollection(\n", + " verts[faces],\n", + " alpha=0.25,\n", + " edgecolor=None,\n", + " linewidth=0.1,\n", + " )\n", + " mesh.set_facecolor(cmap(_i / len(thresholds))[:3])\n", + " mesh.set_edgecolor([1, 0, 0])\n", + " ax.add_collection3d(mesh)\n", + "\n", + " ax.set_xlim(0, p.shape[0])\n", + " ax.set_ylim(0, p.shape[1])\n", + " ax.set_zlim(0, p.shape[2])\n", + "\n", + " ax.view_init(45, 45)\n", + " return _fig" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "fig = show_3d_mesh(bubble_label_image, range(1, np.max(bubble_label_image), 10))\n", + "\n", + "# Write 3D rendering of segmented image\n", + "fig.savefig(\"Volume_rendering.png\")" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "## Calculate bubble centers" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "def meshgrid3d_like(in_img):\n", + " return np.meshgrid(\n", + " range(in_img.shape[1]),range(in_img.shape[0]), range(in_img.shape[2])\n", + " )\n", + "\n", + "zz, xx, yy = meshgrid3d_like(bubble_label_image)\n", + "\n", + "out_results = []\n", + "for c_label in np.unique(bubble_label_image): # one bubble at a time\n", + " if c_label > 0: # ignore background\n", + " cur_roi = bubble_label_image == c_label\n", + " out_results += [\n", + " {\n", + " \"x\": xx[cur_roi].mean(),\n", + " \"y\": yy[cur_roi].mean(),\n", + " \"z\": zz[cur_roi].mean(),\n", + " \"volume\": np.sum(cur_roi),\n", + " }\n", + " ]\n", + "\n", + "# Write the bubble volume stats\n", + "import pandas as pd\n", + "out_table = pd.DataFrame(out_results)\n", + "out_table.to_csv(\"bubble_volume_out.csv\")\n", + "\n", + "# Write the bubble volume stats sample table\n", + "volume_sample = out_table.sample(5)\n", + "volume_sample.save(\"Bubble_volume_stats_sample.ong\")" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "# Write the bubble volume density plot\n", + "bubble_volume_density = out_table[\"volume\"].plot.density()\n", + "bubble_volume_density.save(\"Bubble_volume_density_stats.png\")" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "# Write the bubble center plot\n", + "bubble_centers_plot = out_table.plot.hexbin(\"x\", \"y\", gridsize=(5, 5))\n", + "bubble_centers_plot.save(\"Bubble_centers.png\")" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "markdown", + "source": [ + "## Compare with the training values" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "train_values = pd.read_csv(\"bubble_volume.csv\")" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))\n", + "_, n_bins, _ = ax1.hist(\n", + " np.log10(train_values[\"volume\"]), bins=20, label=\"Training volumes\"\n", + ")\n", + "ax1.hist(np.log10(\n", + " out_table[\"volume\"]),\n", + " n_bins,\n", + " alpha=0.5,\n", + " label=\"Watershed volumes\",\n", + ")\n", + "ax1.legend()\n", + "ax1.set_title(\"Volume comparison\\n(Log10)\")\n", + "ax2.plot(\n", + " out_table[\"x\"],\n", + " out_table[\"y\"],\n", + " \"r.\",\n", + " train_values[\"x\"],\n", + " train_values[\"y\"],\n", + " \"b.\",\n", + ")\n", + "ax2.legend([\"Watershed bubbles\", \"Training bubbles\"])\n", + "\n", + "fig.savefig(\"Bubble_stats.png\")" + ], + "metadata": { + "collapsed": false, + "pycharm": { + "name": "#%%\n" + } + } + } + ], + "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 +} diff --git a/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/Code.py b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/Code.py new file mode 100755 index 000000000..133ec8ab2 --- /dev/null +++ b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/Code.py @@ -0,0 +1,129 @@ +#!/usr/bin/env python + +# 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. + +# Run with: +# ./Code.py +# +# +# +# e.g. +# ./Code.py 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.) +# +# threshold: absolute minimum height value used during processing. +# Raising this threshold percentage effectively decreases the number of local minima in the input, +# resulting in an initial segmentation with fewer regions. +# The assumption is that the shallow regions that thresholding removes are of of less interest. +# level: controls the depth of metaphorical flooding of the image. +# That is, it sets the maximum saliency value of interest in the result. +# Raising and lowering the Level influences the number of segments +# in the basic segmentation that are merged to produce the final output. +# A level of 1.0 is analogous to flooding the image up to a +# depth that is 100 percent of the maximum value in the image. +# A level of 0.0 produces the basic segmentation, which will typically be very oversegmented. +# Level values of interest are typically low (i.e. less than about 0.40 or 40%), +# since higher values quickly start to undersegment the image. + + +import argparse + +import itk +import numpy as np + +parser = argparse.ArgumentParser( + description="Segment an image using the watershed method and the signed Maurer distance map." +) +parser.add_argument("input_image_filename") +parser.add_argument("reversed_input_output_filename") +parser.add_argument("distance_map_output_filename") +parser.add_argument("watershed_output_filename") +parser.add_argument("segmentation_result_output_filename") +parser.add_argument("binarizing_radius", type=int) +parser.add_argument("majority_threshold", type=int) +parser.add_argument("watershed_threshold", type=float) +parser.add_argument("level", type=float) +parser.add_argument("cleaning_structuring_element_radius", type=int) +args = parser.parse_args() + +dimension = 3 + +uchar_pixel_type = itk.UC +uchar_image_type = itk.Image[uchar_pixel_type, dimension] + +float_pixel_type = itk.F +float_image_type = itk.Image[float_pixel_type, dimension] + +rgb_pixel_type = itk.RGBPixel[uchar_pixel_type] +RGBImageType = itk.Image[rgb_pixel_type, dimension] + +stack_image = itk.imread(args.input_image_filename) + +# Create bubble image: get a binarized version of the input image +index_radius = itk.Size[dimension]() +index_radius.Fill(args.binarizing_radius) + +bubble_image = itk.voting_binary_iterative_hole_filling_image_filter( + stack_image, + radius=index_radius, + background_value=0, + foreground_value=255, + majority_threshold=args.majority_threshold, +) + +# Write bubble image +itk.imwrite(bubble_image, args.reversed_input_output_filename) + +# Watershed on bubbles +bbl_cast_image = itk.cast_image_filter(bubble_image, ttype=(uchar_image_type, float_image_type)) + +# Normalize the image to the [0, 255] range +bubble_image_preclamp = itk.multiply_image_filter(bbl_cast_image, constant=255.0) + +bubble_image_clamp = itk.clamp_image_filter(bubble_image_preclamp, bounds=(0, 255)) + +# Get the distance map of the input image +distance_map_image = itk.signed_maurer_distance_map_image_filter( + bubble_image_clamp, inside_is_positive=False +) + +itk.imwrite(distance_map_image, args.distance_map_output_filename) + +# Apply the watershed segmentation +watershed_image = itk.watershed_image_filter( + distance_map_image, threshold=args.watershed_threshold, level=args.level, +) + +# Reduce the maximum label number by using sequential labels +relabeled_image = itk.relabel_component_image_filter(watershed_image) + +# Debugging output +itk.imwrite(relabeled_image, args.watershed_output_filename + '.nrrd') + +# Cast to unsigned short so that it can be written as a TIFF image +label_image_type = itk.Image[itk.US, dimension] +ws_cast_image = itk.cast_image_filter(relabeled_image, ttype=(relabeled_image.__class__, label_image_type)) +itk.imwrite(ws_cast_image, args.watershed_output_filename) + +# Clean the segmentation image: remove small objects by performing an +# opening morphological operation +structuring_element_type = itk.FlatStructuringElement[dimension] +structuring_element = structuring_element_type.Ball(args.cleaning_structuring_element_radius) + +segmented_clean_image = itk.binary_morphological_opening_image_filter(ws_cast_image, kernel=structuring_element) + +itk.imwrite(segmented_clean_image, args.segmentation_result_output_filename) diff --git a/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/Documentation.rst b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/Documentation.rst new file mode 100644 index 000000000..651389107 --- /dev/null +++ b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/Documentation.rst @@ -0,0 +1,126 @@ +Watershed with Distance Map +=========================== + +.. index:: + single: WatershedWithDistanceMap + +Synopsis +-------- + +This example illustrates how to segment an image using the watershed method +and the signed Maurer distance map. + +Results +------- + +.. figure:: BubbleImage.png + :scale: 50% + :alt: Input image + + Input image + +.. figure:: SegmentWithWatershedAndDistanceMapTest01Baseline.png + :scale: 50% + :alt: Segmented image + + Segmented image + +Input image +^^^^^^^^^^^ + +.. figure:: PlateauBorderXYProjection.png + :scale: 50% + :alt: Plateau border XY projection image + + Plateau border XY projection image + +Pipeline images +^^^^^^^^^^^^^^^ + +.. figure:: PipelineImages.png + :scale: 50% + :alt: Input image, distance image and watershed segmentation image + + Input image, distance image and watershed segmentation image + +.. figure:: BubbleImageSegmentedProjectionsImage.png + :scale: 50% + :alt: Segmented bubble projection images nipy spectral colormap + + Segmented bubble projection images nipy spectral colormap + +Clean and relabel +^^^^^^^^^^^^^^^^^ + +.. figure:: CleanAndRelabelProjectionImages.png + :scale: 50% + :alt: Segmented bubble projection images after cleaning and consecutive labelling + + Segmented bubble projection images after cleaning and consecutive labelling + +3D rendering +^^^^^^^^^^^^ + +.. figure:: BubbleImageSegmentedVolume3D.png + :scale: 50% + :alt: Segmented bubble 3D image + + Segmented bubble 3D image + +Bubble statistics +^^^^^^^^^^^^^^^^^ + +.. figure:: BubbleVolumeStatsSampleTable.png + :scale: 50% + :alt: Bubble volume stats sample table + + Bubble volume stats sample table + +.. figure:: BubbleVolumeDensityStats.png + :scale: 50% + :alt: Bubble volume densities + + Bubble volume densities + +.. figure:: BubbleCenters.png + :scale: 50% + :alt: Bubble centers + + Bubble centers + +Comparison +^^^^^^^^^^ + +.. figure:: BubbleStatsComparison.png + :scale: 50% + :alt: Bubble volume histogram and center location comparison + + Bubble volume histogram and center location comparison + +Jupyter Notebook +---------------- + +.. image:: https://mybinder.org/badge_logo.svg + :target: https://mybinder.org/v2/gh/InsightSoftwareConsortium/ITKSphinxExamples/master?filepath=src%2FSegmentation%2FWatersheds%2FSegmentWithWatershedAndDistanceMap%2FCode.ipynb + +Code +---- + +Python +...... + +.. literalinclude:: Code.py + :language: python + :lines: 1, 18- + +C++ +... + +.. literalinclude:: Code.cxx + :lines: 18- + +Classes demonstrated +-------------------- + +.. breathelink:: itk::SignedMaurerDistanceMapImageFilter +.. breathelink:: itk::WatershedImageFilter diff --git a/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/PlateauBorder.tif.sha512 b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/PlateauBorder.tif.sha512 new file mode 100644 index 000000000..9eb69a053 --- /dev/null +++ b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/PlateauBorder.tif.sha512 @@ -0,0 +1 @@ +1424a5e7a304d8cbb830321182f5d77742acdfe4bc3d72cf58e079cf994e4d51bce1ec88530865a98292c163ef8487dec720d84b128a47c277e29141facdc006 diff --git a/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/segmentation_result_analysis.py b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/segmentation_result_analysis.py new file mode 100755 index 000000000..02b5f0500 --- /dev/null +++ b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/segmentation_result_analysis.py @@ -0,0 +1,103 @@ +#!/usr/bin/env python + +# 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. + +# Run with: +# ./segmentation_result_analysis.py <> +# e.g. +# ./segmentation_result_analysis.py P +# +# This is a companion code to the SegmentWithWatershedAndDistanceMap example +# to generate the example segmentation results analysis visualizations. +# + + +import argparse + +import itk +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + + +parser = argparse.ArgumentParser(description="Analyze the segmentation result.") +parser.add_argument("clean_segmentation_filename") +parser.add_argument("training_bubble_volume_stats_filename") +parser.add_argument("bubble_volume_stats_output_filename") +parser.add_argument("bubble_volume_stats_sample_table_output_filename") +parser.add_argument("bubble_volume_density_stats_plot_output_filename") +parser.add_argument("bubble_centers_plot_output_filename") +parser.add_argument("bubble_stats_comparison_plot_output_filename") +args = parser.parse_args() + +# Read segmentation image +bubble_label_image = itk.imread(args.clean_segmentation_filename) + + +# Calculate bubble centers +def meshgrid3d_like(in_img): + return np.meshgrid( + range(in_img.shape[1]), range(in_img.shape[0]), range(in_img.shape[2]) + ) + + +zz, xx, yy = meshgrid3d_like(bubble_label_image) +out_results = [] +for c_label in np.unique(bubble_label_image): # one bubble at a time + if c_label > 0: # ignore background + cur_roi = bubble_label_image == c_label + out_results += [ + { + "x": xx[cur_roi].mean(), + "y": yy[cur_roi].mean(), + "z": zz[cur_roi].mean(), + "volume": np.sum(cur_roi), + } + ] + +# Write the bubble volume stats +out_table = pd.DataFrame(out_results) +out_table.to_csv(args.bubble_volume_stats_output_filename) # "bubble_volume.csv" + +# Write the bubble volume stats sample table +volume_sample = out_table.sample(5) +volume_sample.save(args.bubble_volume_stats_sample_table_output_filename) + +# Write the bubble volume density plot +bubble_volume_density = out_table["volume"].plot.density() +bubble_volume_density.save(args.bubble_volume_density_stats_plot_output_filename) + +# Write the bubble center plot +bubble_centers_plot = out_table.plot.hexbin("x", "y", gridsize=(5, 5)) +bubble_centers_plot.save(args.bubble_centers_plot_output_filename) + +# Compare to the training values +train_values = pd.read_csv( + args.training_bubble_volume_stats_filename +) # "../input/bubble_volume.csv" + +fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4)) +_, n_bins, _ = ax1.hist( + np.log10(train_values["volume"]), bins=20, label="Training volumes" +) +ax1.hist(np.log10(out_table["volume"]), n_bins, alpha=0.5, label="Watershed volumes") +ax1.legend() +ax1.set_title("Volume comparison\n(Log10)") +ax2.plot( + out_table["x"], out_table["y"], "r.", train_values["x"], train_values["y"], "b." +) +ax2.legend(["Watershed bubbles", "Training bubbles"]) + +fig.savefig(args.bubble_stats_comparison_plot_output_filename) diff --git a/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/segmentation_result_visualization.py b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/segmentation_result_visualization.py new file mode 100755 index 000000000..302509147 --- /dev/null +++ b/src/Segmentation/Watersheds/SegmentWithWatershedAndDistanceMap/segmentation_result_visualization.py @@ -0,0 +1,143 @@ +#!/usr/bin/env python + +# 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. + +# Run with: +# ./segmentation_result_visualization.py <> +# e.g. +# ./segmentation_result_visualization.py P +# +# This is a companion code to the SegmentWithWatershedAndDistanceMap example +# to generate the example segmentation results visualizations. +# + + +import argparse + +import itk +import matplotlib.pyplot as plt +import numpy as np +from mpl_toolkits.mplot3d.art3d import Poly3DCollection +from skimage import measure +from tqdm import tqdm + + +parser = argparse.ArgumentParser( + description="Visualize the segmentation result." +) +parser.add_argument("input_image_filename") +parser.add_argument("bubble_filename") +parser.add_argument("distance_map_filename") +parser.add_argument("watershed_filename") +parser.add_argument("clean_segmentation_filename") +parser.add_argument("input_projections_output_filename") +parser.add_argument("pipeline_images_output_filename") +parser.add_argument("segmentation_projections_output_filename") +parser.add_argument("clean_segmentation_projections_output_filename") +parser.add_argument("volume_rendering_output_filename") +args = parser.parse_args() + + +stack_image = itk.imread(args.input_image_filename) + +# Display the input image projections +fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 4)) +for i, (cax, clabel) in enumerate(zip([ax1, ax2, ax3], ["xy", "zy", "zx"])): + cax.imshow(np.sum(stack_image, i).squeeze(), cmap="bone_r") + cax.set_title("{} projection".format(clabel)) + cax.set_xlabel(clabel[0]) + cax.set_ylabel(clabel[1]) + +fig.savefig(args.input_projections_output_filename) + +bubble_image = itk.imread(args.bubble_filename) + +plt.imshow(bubble_image[5], cmap="bone") +fig.savefig(args.bubble_filename) + +# Cut-through: display the segmentation pipeline image results +# Show the values at a slice in the middle as a way to get a feeling for what +# the distance map and the watershed did + +ws_vol = itk.imread(args.watershed_filename) +ws_vol_arr = itk.array_view_from_image(ws_vol) + +dmap_vol = itk.imread(args.distance_map_filename) + +mid_slice = ws_vol.shape[0] // 2 +fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 7)) +ax1.imshow(bubble_image[mid_slice], cmap="bone") +ax1.set_title("Bubble image") +m_val = np.abs(dmap_vol[mid_slice]).std() +ax2.imshow(dmap_vol[mid_slice], cmap="RdBu", vmin=-m_val, vmax=m_val) +ax2.set_title("Distance image\nmin: {:2.2f}; max: {:2.2f}; mean: {:2.2f}".format( + dmap_vol[mid_slice].min(), dmap_vol[mid_slice].max(), dmap_vol[mid_slice].mean())) +ax3.imshow(ws_vol[mid_slice], cmap="nipy_spectral") +ax3.set_title("Watershed\nLabels found: {}".format(len(np.unique(ws_vol_arr[ws_vol_arr > 0])))) + +fig.savefig(args.pipeline_images_output_filename) + +# Show segmentation result projections +fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 4)) +for i, (cax, clabel) in enumerate(zip([ax1, ax2, ax3], ["xy", "zy", "zx"])): + cax.imshow(np.max(ws_vol, i).squeeze(), cmap="nipy_spectral") + cax.set_title("{} projection".format(clabel)) + cax.set_xlabel(clabel[0]) + cax.set_ylabel(clabel[1]) + +fig.savefig(args.segmentation_projections_output_filename) + +# Clean segmentation result projections +bubble_label_image = itk.imread(args.clean_segmentation_filename) + +# ToDo +# Count the kept bubbles in bubble_label_image +new_idx = 1234 + +fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 4)) +for i, (cax, clabel) in enumerate(zip([ax1, ax2, ax3], ["xy", "zy", "zx"])): + cax.imshow(np.max(bubble_label_image, i), cmap="jet", vmin=0, vmax=new_idx) + cax.set_title("{} projection".format(clabel)) + cax.set_xlabel(clabel[0]) + cax.set_ylabel(clabel[1]) + +fig.savefig(args.clean_segmentation_projections_output_filename) + + +# Show 3D rendering +def show_3d_mesh(image, thresholds): + p = image[::-1].swapaxes(1, 2) + cmap = plt.cm.get_cmap("nipy_spectral_r") + _fig = plt.figure(figsize=(10, 10)) + ax = _fig.add_subplot(111, projection="3d") + for _i, c_threshold in tqdm(list(enumerate(thresholds))): + verts, faces, _, _ = measure.marching_cubes(p == c_threshold, level=0) + mesh = Poly3DCollection(verts[faces], alpha=0.25, edgecolor=None, linewidth=0.1) + mesh.set_facecolor(cmap(_i / len(thresholds))[:3]) + mesh.set_edgecolor([1, 0, 0]) + ax.add_collection3d(mesh) + + ax.set_xlim(0, p.shape[0]) + ax.set_ylim(0, p.shape[1]) + ax.set_zlim(0, p.shape[2]) + + ax.view_init(45, 45) + return _fig + + +fig = show_3d_mesh(bubble_label_image, range(1, np.max(bubble_label_image), 10)) + +# Write 3D rendering of segmented image +fig.savefig(args.volume_rendering_output_filename) diff --git a/src/Segmentation/Watersheds/index.rst b/src/Segmentation/Watersheds/index.rst index 33105cfc3..1ccb83400 100644 --- a/src/Segmentation/Watersheds/index.rst +++ b/src/Segmentation/Watersheds/index.rst @@ -5,4 +5,5 @@ Watersheds :maxdepth: 1 MorphologicalWatershedSegmentation/Documentation.rst + SegmentWithWatershedAndDistanceMap/Documentation.rst SegmentWithWatershedImageFilter/Documentation.rst