From 660ce50945f6804647261db7422fe834c37781ac Mon Sep 17 00:00:00 2001 From: Matt McCormick Date: Fri, 14 May 2021 16:49:45 -0400 Subject: [PATCH 1/2] STYLE: GeodesicActiveContour CMake style --- .../CMakeLists.txt | 48 +++++++++---------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/CMakeLists.txt b/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/CMakeLists.txt index b3e8a596c..a8a679548 100644 --- a/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/CMakeLists.txt +++ b/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/CMakeLists.txt @@ -1,19 +1,19 @@ cmake_minimum_required(VERSION 3.10.2) -project( SegmentWithGeodesicActiveContourLevelSet ) +project(SegmentWithGeodesicActiveContourLevelSet) -find_package( ITK REQUIRED ) -include( ${ITK_USE_FILE} ) +find_package(ITK REQUIRED) +include(${ITK_USE_FILE}) -add_executable( SegmentWithGeodesicActiveContourLevelSet Code.cxx ) -target_link_libraries( SegmentWithGeodesicActiveContourLevelSet ${ITK_LIBRARIES} ) +add_executable(SegmentWithGeodesicActiveContourLevelSet Code.cxx) +target_link_libraries(SegmentWithGeodesicActiveContourLevelSet ${ITK_LIBRARIES}) -install( TARGETS SegmentWithGeodesicActiveContourLevelSet +install(TARGETS SegmentWithGeodesicActiveContourLevelSet DESTINATION bin/ITKSphinxExamples/Segmentation/LevelSets COMPONENT Runtime ) -install( FILES Code.cxx CMakeLists.txt +install(FILES Code.cxx CMakeLists.txt DESTINATION share/ITKSphinxExamples/Code/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet COMPONENT Code ) @@ -21,34 +21,34 @@ install( FILES Code.cxx CMakeLists.txt enable_testing() add_test(NAME GeodesicActiveContourLeftVentricleTest - COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/SegmentWithGeodesicActiveContourLevelSet - ${CMAKE_CURRENT_BINARY_DIR}/BrainProtonDensitySlice.png - Output1.png - 81 114 5.0 1.0 -0.5 3.0 2.0 800 + COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/SegmentWithGeodesicActiveContourLevelSet + ${CMAKE_CURRENT_BINARY_DIR}/BrainProtonDensitySlice.png + Output1.png + 81 114 5.0 1.0 -0.5 3.0 2.0 800 ) add_test(NAME GeodesicActiveContourRightVentricleTest - COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/SegmentWithGeodesicActiveContourLevelSet - ${CMAKE_CURRENT_BINARY_DIR}/BrainProtonDensitySlice.png - Output2.png - 99 114 5.0 1.0 -0.5 3.0 2.0 800 + COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/SegmentWithGeodesicActiveContourLevelSet + ${CMAKE_CURRENT_BINARY_DIR}/BrainProtonDensitySlice.png + Output2.png + 99 114 5.0 1.0 -0.5 3.0 2.0 800 ) add_test(NAME GeodesicActiveContourWhiteMatterTest - COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/SegmentWithGeodesicActiveContourLevelSet - ${CMAKE_CURRENT_BINARY_DIR}/BrainProtonDensitySlice.png - Output3.png - 56 92 5.0 1.0 -0.3 2.0 10.0 800 + COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/SegmentWithGeodesicActiveContourLevelSet + ${CMAKE_CURRENT_BINARY_DIR}/BrainProtonDensitySlice.png + Output3.png + 56 92 5.0 1.0 -0.3 2.0 10.0 800 ) add_test(NAME GeodesicActiveContourGrayMatterTest - COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/SegmentWithGeodesicActiveContourLevelSet - ${CMAKE_CURRENT_BINARY_DIR}/BrainProtonDensitySlice.png - Output4.png - 40 90 5.0 .5 -0.3 2.0 10.0 800 + COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/SegmentWithGeodesicActiveContourLevelSet + ${CMAKE_CURRENT_BINARY_DIR}/BrainProtonDensitySlice.png + Output4.png + 40 90 5.0 .5 -0.3 2.0 10.0 800 ) -if( ITK_WRAP_PYTHON ) +if(ITK_WRAP_PYTHON) add_test(NAME GeodesicActiveContourLeftVentricleTestPython COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/Code.py ${CMAKE_CURRENT_BINARY_DIR}/BrainProtonDensitySlice.png From 594ae039423e8b66db69efdab81ccf3558c886c6 Mon Sep 17 00:00:00 2001 From: Matt McCormick Date: Fri, 14 May 2021 17:44:19 -0400 Subject: [PATCH 2/2] DOC: Add notebook example for SegmentWithGeodesicActiveContourLevelSet Also update the Python code style to be more Pythonic --- .../CMakeLists.txt | 2 +- .../Code.ipynb | 195 ++++++++++++++++ .../Code.py | 214 ++++++------------ .../Documentation.rst | 5 + 4 files changed, 271 insertions(+), 145 deletions(-) create mode 100644 src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Code.ipynb diff --git a/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/CMakeLists.txt b/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/CMakeLists.txt index a8a679548..81a0adf9f 100644 --- a/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/CMakeLists.txt +++ b/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/CMakeLists.txt @@ -13,7 +13,7 @@ install(TARGETS SegmentWithGeodesicActiveContourLevelSet COMPONENT Runtime ) -install(FILES Code.cxx CMakeLists.txt +install(FILES Code.cxx CMakeLists.txt Code.py DESTINATION share/ITKSphinxExamples/Code/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet COMPONENT Code ) diff --git a/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Code.ipynb b/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Code.ipynb new file mode 100644 index 000000000..bee961a19 --- /dev/null +++ b/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Code.ipynb @@ -0,0 +1,195 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "7c8b935b", + "metadata": {}, + "source": [ + "# Segment with geodesic active contour level set" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "196c6071", + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "import os\n", + "from urllib.request import urlretrieve\n", + "\n", + "import itk\n", + "\n", + "from itkwidgets import view" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d6f35e01", + "metadata": {}, + "outputs": [], + "source": [ + "input_filename = 'BrainProtonDensitySlice.png'\n", + "if not os.path.exists(input_filename):\n", + " url = 'https://data.kitware.com/api/v1/file/57b5d8028d777f10f2694bbf/download'\n", + " urlretrieve(url, input_filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "52104b31", + "metadata": {}, + "outputs": [], + "source": [ + "InputPixelType = itk.ctype('float')\n", + "\n", + "input_image = itk.imread(input_filename, InputPixelType)\n", + "\n", + "view(input_image)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4297e2ac", + "metadata": {}, + "outputs": [], + "source": [ + "smoothed = itk.curvature_anisotropic_diffusion_image_filter(input_image,\n", + " time_step=0.125,\n", + " number_of_iterations=5,\n", + " conductance_parameter=9.0)\n", + "view(smoothed)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "99186583", + "metadata": {}, + "outputs": [], + "source": [ + "sigma = 1.0\n", + "\n", + "gradient_magnitude = itk.gradient_magnitude_recursive_gaussian_image_filter(smoothed,\n", + " sigma=sigma)\n", + "view(gradient_magnitude)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8a219f34", + "metadata": {}, + "outputs": [], + "source": [ + "alpha = -0.5\n", + "beta = 3.0\n", + "\n", + "sigmoid = itk.sigmoid_image_filter(gradient_magnitude,\n", + " output_minimum=0.0,\n", + " output_maximum=1.0,\n", + " alpha=alpha,\n", + " beta=beta)\n", + "\n", + "view(sigmoid)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "499443a8", + "metadata": {}, + "outputs": [], + "source": [ + "Dimension = input_image.GetImageDimension()\n", + "seeds = itk.VectorContainer[itk.UI, itk.LevelSetNode[InputPixelType, Dimension]].New()\n", + "seeds.Initialize()\n", + "\n", + "seed_position = itk.Index[Dimension]()\n", + "seed_position[0] = 81\n", + "seed_position[1] = 114\n", + "node = itk.LevelSetNode[InputPixelType, Dimension]()\n", + "node.SetValue(-5.0)\n", + "node.SetIndex(seed_position)\n", + "seeds.InsertElement(0, node)\n", + "\n", + "fast_marching = itk.fast_marching_image_filter(trial_points=seeds,\n", + " speed_constant=1.0,\n", + " output_size=input_image.GetBufferedRegion().GetSize())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2c5e205", + "metadata": {}, + "outputs": [], + "source": [ + "propagation_scaling = 2.0\n", + "number_of_iterations = 800\n", + "\n", + "geodesic_active_contour = \\\n", + " itk.geodesic_active_contour_level_set_image_filter(fast_marching,\n", + " propagation_scaling=propagation_scaling,\n", + " curvature_scaling=1.0,\n", + " advection_scaling=1.0,\n", + " maximum_r_m_s_error=0.02,\n", + " number_of_iterations=number_of_iterations,\n", + " feature_image=sigmoid)\n", + "\n", + "view(geodesic_active_contour)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1269f1b4", + "metadata": {}, + "outputs": [], + "source": [ + "OutputPixelType = itk.ctype('unsigned char')\n", + "thresholded = itk.binary_threshold_image_filter(geodesic_active_contour,\n", + " lower_threshold=-1000.0,\n", + " upper_threshold=0.0,\n", + " outside_value=itk.NumericTraits[OutputPixelType].min(),\n", + " inside_value=itk.NumericTraits[OutputPixelType].max(),\n", + " ttype=[type(geodesic_active_contour), itk.Image[OutputPixelType,Dimension]])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "58bdf582", + "metadata": {}, + "outputs": [], + "source": [ + "view(thresholded)" + ] + } + ], + "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/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Code.py b/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Code.py index 6b936bf08..79cdc0ae2 100755 --- a/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Code.py +++ b/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Code.py @@ -25,169 +25,95 @@ ) sys.exit(1) -inputFileName = sys.argv[1] +input_filename = sys.argv[1] +output_filename = sys.argv[2] +seed_pos_x = int(sys.argv[3]) +seed_pos_y = int(sys.argv[4]) -outputFileName = sys.argv[2] -seedPosX = int(sys.argv[3]) -seedPosY = int(sys.argv[4]) - -initialDistance = float(sys.argv[5]) +initial_distance = float(sys.argv[5]) sigma = float(sys.argv[6]) alpha = float(sys.argv[7]) beta = float(sys.argv[8]) -propagationScaling = float(sys.argv[9]) -numberOfIterations = int(sys.argv[10]) -seedValue = -initialDistance +propagation_scaling = float(sys.argv[9]) +number_of_iterations = int(sys.argv[10]) +seed_value = -initial_distance -Dimension = 2 +InputPixelType = itk.ctype("float") +input_image = itk.imread(input_filename, InputPixelType) -InputPixelType = itk.F -OutputPixelType = itk.UC +Dimension = input_image.GetImageDimension() -InputImageType = itk.Image[InputPixelType, Dimension] -OutputImageType = itk.Image[OutputPixelType, Dimension] +smoothed = itk.curvature_anisotropic_diffusion_image_filter( + input_image, time_step=0.125, number_of_iterations=5, conductance_parameter=9.0 +) -ReaderType = itk.ImageFileReader[InputImageType] -WriterType = itk.ImageFileWriter[OutputImageType] - -reader = ReaderType.New() -reader.SetFileName(inputFileName) - -SmoothingFilterType = itk.CurvatureAnisotropicDiffusionImageFilter[ - InputImageType, InputImageType -] -smoothing = SmoothingFilterType.New() -smoothing.SetTimeStep(0.125) -smoothing.SetNumberOfIterations(5) -smoothing.SetConductanceParameter(9.0) -smoothing.SetInput(reader.GetOutput()) - -GradientFilterType = itk.GradientMagnitudeRecursiveGaussianImageFilter[ - InputImageType, InputImageType -] -gradientMagnitude = GradientFilterType.New() -gradientMagnitude.SetSigma(sigma) -gradientMagnitude.SetInput(smoothing.GetOutput()) - -SigmoidFilterType = itk.SigmoidImageFilter[InputImageType, InputImageType] -sigmoid = SigmoidFilterType.New() -sigmoid.SetOutputMinimum(0.0) -sigmoid.SetOutputMaximum(1.0) -sigmoid.SetAlpha(alpha) -sigmoid.SetBeta(beta) -sigmoid.SetInput(gradientMagnitude.GetOutput()) - -FastMarchingFilterType = itk.FastMarchingImageFilter[InputImageType, InputImageType] -fastMarching = FastMarchingFilterType.New() - -GeoActiveContourFilterType = itk.GeodesicActiveContourLevelSetImageFilter[ - InputImageType, InputImageType, InputPixelType -] -geodesicActiveContour = GeoActiveContourFilterType.New() -geodesicActiveContour.SetPropagationScaling(propagationScaling) -geodesicActiveContour.SetCurvatureScaling(1.0) -geodesicActiveContour.SetAdvectionScaling(1.0) -geodesicActiveContour.SetMaximumRMSError(0.02) -geodesicActiveContour.SetNumberOfIterations(numberOfIterations) -geodesicActiveContour.SetInput(fastMarching.GetOutput()) -geodesicActiveContour.SetFeatureImage(sigmoid.GetOutput()) - -ThresholdingFilterType = itk.BinaryThresholdImageFilter[InputImageType, OutputImageType] -thresholder = ThresholdingFilterType.New() -thresholder.SetLowerThreshold(-1000.0) -thresholder.SetUpperThreshold(0.0) -thresholder.SetOutsideValue(itk.NumericTraits[OutputPixelType].min()) -thresholder.SetInsideValue(itk.NumericTraits[OutputPixelType].max()) -thresholder.SetInput(geodesicActiveContour.GetOutput()) - -seedPosition = itk.Index[Dimension]() -seedPosition[0] = seedPosX -seedPosition[1] = seedPosY +gradient_magnitude = itk.gradient_magnitude_recursive_gaussian_image_filter( + smoothed, sigma=sigma +) -node = itk.LevelSetNode[InputPixelType, Dimension]() -node.SetValue(seedValue) -node.SetIndex(seedPosition) +sigmoid = itk.sigmoid_image_filter( + gradient_magnitude, output_minimum=0.0, output_maximum=1.0, alpha=alpha, beta=beta +) +Dimension = input_image.GetImageDimension() seeds = itk.VectorContainer[itk.UI, itk.LevelSetNode[InputPixelType, Dimension]].New() seeds.Initialize() + +seed_position = itk.Index[Dimension]() +seed_position[0] = seed_pos_x +seed_position[1] = seed_pos_y +node = itk.LevelSetNode[InputPixelType, Dimension]() +node.SetValue(-5.0) +node.SetIndex(seed_position) seeds.InsertElement(0, node) -fastMarching.SetTrialPoints(seeds) -fastMarching.SetSpeedConstant(1.0) - -CastFilterType = itk.RescaleIntensityImageFilter[InputImageType, OutputImageType] - -caster1 = CastFilterType.New() -caster2 = CastFilterType.New() -caster3 = CastFilterType.New() -caster4 = CastFilterType.New() - -writer1 = WriterType.New() -writer2 = WriterType.New() -writer3 = WriterType.New() -writer4 = WriterType.New() - -caster1.SetInput(smoothing.GetOutput()) -writer1.SetInput(caster1.GetOutput()) -writer1.SetFileName("GeodesicActiveContourImageFilterOutput1.png") -caster1.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min()) -caster1.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max()) -writer1.Update() - -caster2.SetInput(gradientMagnitude.GetOutput()) -writer2.SetInput(caster2.GetOutput()) -writer2.SetFileName("GeodesicActiveContourImageFilterOutput2.png") -caster2.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min()) -caster2.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max()) -writer2.Update() - -caster3.SetInput(sigmoid.GetOutput()) -writer3.SetInput(caster3.GetOutput()) -writer3.SetFileName("GeodesicActiveContourImageFilterOutput3.png") -caster3.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min()) -caster3.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max()) -writer3.Update() - -caster4.SetInput(fastMarching.GetOutput()) -writer4.SetInput(caster4.GetOutput()) -writer4.SetFileName("GeodesicActiveContourImageFilterOutput4.png") -caster4.SetOutputMinimum(itk.NumericTraits[OutputPixelType].min()) -caster4.SetOutputMaximum(itk.NumericTraits[OutputPixelType].max()) - -fastMarching.SetOutputSize(reader.GetOutput().GetBufferedRegion().GetSize()) - -writer = WriterType.New() -writer.SetFileName(outputFileName) -writer.SetInput(thresholder.GetOutput()) -writer.Update() - -print( - "Max. no. iterations: " + str(geodesicActiveContour.GetNumberOfIterations()) + "\n" +fast_marching = itk.fast_marching_image_filter( + trial_points=seeds, + speed_constant=1.0, + output_size=input_image.GetBufferedRegion().GetSize(), ) -print("Max. RMS error: " + str(geodesicActiveContour.GetMaximumRMSError()) + "\n") -print( - "No. elpased iterations: " - + str(geodesicActiveContour.GetElapsedIterations()) - + "\n" + +geodesic_active_contour = itk.geodesic_active_contour_level_set_image_filter( + fast_marching, + propagation_scaling=propagation_scaling, + curvature_scaling=1.0, + advection_scaling=1.0, + maximum_r_m_s_error=0.02, + number_of_iterations=number_of_iterations, + feature_image=sigmoid, ) -print("RMS change: " + str(geodesicActiveContour.GetRMSChange()) + "\n") -writer4.Update() -InternalWriterType = itk.ImageFileWriter[InputImageType] +OutputPixelType = itk.ctype("unsigned char") +OutputImageType = itk.Image[OutputPixelType, Dimension] +thresholded = itk.binary_threshold_image_filter( + geodesic_active_contour, + lower_threshold=-1000.0, + upper_threshold=0.0, + outside_value=itk.NumericTraits[OutputPixelType].min(), + inside_value=itk.NumericTraits[OutputPixelType].max(), + ttype=[type(geodesic_active_contour), OutputImageType], +) + +itk.imwrite(thresholded, output_filename) + +smoothed_cast = itk.rescale_intensity_image_filter( + smoothed, ttype=(type(smoothed), OutputImageType) +) +itk.imwrite(smoothed_cast, "GeodesicActiveContourImageFilterOutput1.png") -mapWriter = InternalWriterType.New() -mapWriter.SetInput(fastMarching.GetOutput()) -mapWriter.SetFileName("GeodesicActiveContourImageFilterOutput4.mha") -mapWriter.Update() +gradient_magnitude_cast = itk.rescale_intensity_image_filter( + smoothed, ttype=(type(gradient_magnitude), OutputImageType) +) +itk.imwrite(gradient_magnitude_cast, "GeodesicActiveContourImageFilterOutput2.png") -speedWriter = InternalWriterType.New() -speedWriter.SetInput(sigmoid.GetOutput()) -speedWriter.SetFileName("GeodesicActiveContourImageFilterOutput3.mha") -speedWriter.Update() +sigmoid_cast = itk.rescale_intensity_image_filter( + sigmoid, ttype=(type(sigmoid), OutputImageType) +) +itk.imwrite(sigmoid_cast, "GeodesicActiveContourImageFilterOutput3.png") -gradientWriter = InternalWriterType.New() -gradientWriter.SetInput(gradientMagnitude.GetOutput()) -gradientWriter.SetFileName("GeodesicActiveContourImageFilterOutput2.mha") -gradientWriter.Update() +fast_marching_cast = itk.rescale_intensity_image_filter( + fast_marching, ttype=(type(fast_marching), OutputImageType) +) +itk.imwrite(fast_marching_cast, "GeodesicActiveContourImageFilterOutput4.png") diff --git a/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Documentation.rst b/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Documentation.rst index 3abd90d39..bbe4c7a70 100644 --- a/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Documentation.rst +++ b/src/Segmentation/LevelSets/SegmentWithGeodesicActiveContourLevelSet/Documentation.rst @@ -43,6 +43,11 @@ Used parameters: * White matter: 56 92 5.0 1.0 -0.3 2.0 10.0 * Gray matter: 40 90 5.0 .5 -0.3 2.0 10.0 +Jupyter Notebook +---------------- + +.. image:: https://mybinder.org/badge_logo.svg + :target: https://mybinder.org/v2/gh/InsightSoftwareConsortium/ITKSphinxExamples/master?filepath=src%2FSegmentation%2FLevelSets%2FSegmentWithGeodesicActiveContourLevelSet%2FCode.ipynb Code ----