Merge pull request #198 from tbirdso/mutual-information-affine-example

ENH: Fix MutualInformationAffineExample and add Binder notebooks

Author: thewtex

requirements.txt

@@ -0,0 +1,67 @@
+async-generator==1.10
+atomicwrites==1.4.0
+backcall==0.2.0
+brotlipy==0.7.0
+certifi==2020.12.5
+colorama==0.4.4
+colorcet==2.0.2
+cycler==0.10.0
+decorator==4.4.2
+defusedxml==0.6.0
+entrypoints==0.3
+imageio==2.9.0
+iniconfig==1.1.1
+ipydatawidgets==4.2.0
+ipykernel==5.4.3
+ipympl==0.6.2
+ipython==7.19.0
+ipython-genutils==0.2.0
+ipywidgets==7.6.3
+itk==5.2rc2
+itk-boneenhancement==0.4.0
+itk-bonemorphometry==1.2.0
+itk-core==5.2rc2
+itk-filtering==5.2rc2
+itk-hasi==0.1.1
+itk-io==5.2rc2
+itk-meshtopolydata==0.6.3
+itk-numerics==5.2rc2
+itk-registration==5.2rc2
+itk-segmentation==5.2rc2
+itkwidgets==0.32.0
+jedi==0.18.0
+Jinja2==2.11.2
+json5==0.9.5
+jupyter-client==6.1.11
+jupyter-core==4.7.0
+jupyterlab==2.2.6
+jupyterlab-widgets==1.0.0
+kiwisolver==1.3.1
+MarkupSafe==1.1.1
+matplotlib==3.3.3
+mistune==0.8.4
+numpy==1.19.5
+param==1.10.1
+parso==0.8.1
+pickleshare==0.7.5
+Pillow==8.1.0
+pluggy==0.13.1
+prompt-toolkit==3.0.10
+py==1.10.0
+pyct==0.4.8
+Pygments==2.7.4
+pyparsing==2.4.7
+pytest==6.2.1
+python-dateutil==2.8.1
+pyzmq==20.0.0
+six==1.15.0
+testpath==0.4.4
+toml==0.10.2
+tornado==6.1
+traitlets==5.0.5
+traittypes==0.2.1
+wcwidth==0.2.5
+webencodings==0.5.1
+widgetsnbextension==3.5.1
+wincertstore==0.2
+zstandard==0.15.1

src/Core/Transform/CMakeLists.txt

@@ -29,15 +29,13 @@ compare_to_baseline(EXAMPLE_NAME GlobalRegistrationTwoImagesAffine
     output
   )
 
-# TODO revisit input/baseline images and example paramters
-#   before reenabling
-#add_example(MutualInformationAffine)
-#compare_to_baseline(EXAMPLE_NAME MutualInformationAffine
-#  BASELINE_PREFIX
-#    fixed
-#    moving
-#    output
-#  )
+add_example(MutualInformationAffine)
+compare_to_baseline(EXAMPLE_NAME MutualInformationAffine
+  BASELINE_PREFIX
+    fixed
+    moving
+    OutputBaseline
+  )
 
 add_example(GlobalRegistrationTwoImagesBSpline)
 compare_to_baseline(EXAMPLE_NAME GlobalRegistrationTwoImagesBSpline

src/Core/Transform/MutualInformationAffine/CMakeLists.txt

@@ -22,5 +22,8 @@ install(FILES Code.cxx CMakeLists.txt
 
 enable_testing()
 add_test(NAME MutualInformationAffineTest
-  COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/MutualInformationAffine)
+  COMMAND ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/MutualInformationAffine
+  fixed.png
+  moving.png
+  output.png)
 

src/Core/Transform/MutualInformationAffine/Code.cxx

@@ -27,38 +27,32 @@
 #include "itkEllipseSpatialObject.h"
 #include "itkSpatialObjectToImageFilter.h"
 #include "itkImageFileWriter.h"
+#include "itkImageFileReader.h"
 
 constexpr unsigned int Dimension = 2;
 using PixelType = unsigned char;
 
 using ImageType = itk::Image<PixelType, Dimension>;
 
-static void
-CreateEllipseImage(ImageType::Pointer image);
-static void
-CreateCircleImage(ImageType::Pointer image);
-
 int
-main(int itkNotUsed(argc), char * itkNotUsed(argv)[])
+main(int argc, char * argv[])
 {
-  // Generate synthetic fixed and moving images
-  ImageType::Pointer fixedImage = ImageType::New();
-  CreateCircleImage(fixedImage);
-  ImageType::Pointer movingImage = ImageType::New();
-  CreateEllipseImage(movingImage);
-
-  // Write the two synthetic inputs
-  using WriterType = itk::ImageFileWriter<ImageType>;
+  using ReaderType = itk::ImageFileReader<ImageType>;
 
-  WriterType::Pointer fixedWriter = WriterType::New();
-  fixedWriter->SetFileName("fixed.png");
-  fixedWriter->SetInput(fixedImage);
-  fixedWriter->Update();
+  if (argc < 4)
+  {
+    std::cout << "Usage: " << argv[0] << " imageFile1 imageFile2 outputFile" << std::endl;
+    return EXIT_FAILURE;
+  }
+  ReaderType::Pointer fixedReader = ReaderType::New();
+  fixedReader->SetFileName(argv[1]);
+  fixedReader->Update();
+  ImageType::Pointer fixedImage = fixedReader->GetOutput();
 
-  WriterType::Pointer movingWriter = WriterType::New();
-  movingWriter->SetFileName("moving.png");
-  movingWriter->SetInput(movingImage);
-  movingWriter->Update();
+  ReaderType::Pointer movingReader = ReaderType::New();
+  movingReader->SetFileName(argv[2]);
+  movingReader->Update();
+  ImageType::Pointer movingImage = movingReader->GetOutput();
 
   // We use floats internally
   using InternalImageType = itk::Image<float, Dimension>;
@@ -113,8 +107,8 @@ main(int itkNotUsed(argc), char * itkNotUsed(argv)[])
   //  which have been normalized to a mean of zero and unit variance.  We
   //  will follow this empirical rule in this example.
 
-  metric->SetFixedImageStandardDeviation(0.4);
-  metric->SetMovingImageStandardDeviation(0.4);
+  metric->SetFixedImageStandardDeviation(5.0);
+  metric->SetMovingImageStandardDeviation(5.0);
 
   registration->SetFixedImage(fixedSmoother->GetOutput());
   registration->SetMovingImage(movingSmoother->GetOutput());
@@ -172,11 +166,8 @@ main(int itkNotUsed(argc), char * itkNotUsed(argv)[])
 
   metric->SetNumberOfSpatialSamples(numberOfSamples);
 
-  // optimizer->SetLearningRate( 15.0 ); //"All the sampled point mapped to outside of the moving image"
-  // optimizer->SetLearningRate( 1.0 );
-  optimizer->SetLearningRate(0.1);
-  optimizer->SetNumberOfIterations(1000);
-  optimizer->MaximizeOn(); // We want to maximize mutual information (the default of the optimizer is to minimize)
+  // For consistent results when regression testing.
+  metric->ReinitializeSeed(121212);
 
   // Note that large values of the learning rate will make the optimizer
   // unstable. Small values, on the other hand, may result in the optimizer
@@ -193,6 +184,30 @@ main(int itkNotUsed(argc), char * itkNotUsed(argv)[])
   // optimizer step length is proportional to the Metric values themselves.
   // Metrics with large values will require you to use smaller values for the
   // learning rate in order to maintain a similar optimizer behavior.
+  optimizer->SetLearningRate(1.0);
+
+  // Note that the only stop condition for the v3 GradientDescentOptimizer class
+  // is that the maximum number of iterations is reached.
+  // For the option to exit early on convergence use GradientDescentOptimizerv4
+  // with an accompanying v4 metric class.
+  optimizer->SetNumberOfIterations(200);
+  optimizer->MaximizeOn(); // We want to maximize mutual information (the default of the optimizer is to minimize)
+
+  auto scales = optimizer->GetScales();
+
+  // Let optimizer take
+  // large steps along translation parameters,
+  // moderate steps along rotational parameters,
+  // and small steps along scale parameters
+  scales.SetSize(6);
+  scales.SetElement(0, 100);
+  scales.SetElement(1, 0.5);
+  scales.SetElement(2, 0.5);
+  scales.SetElement(3, 100);
+  scales.SetElement(4, 0.0001);
+  scales.SetElement(5, 0.0001);
+
+  optimizer->SetScales(scales);
 
   try
   {
@@ -240,107 +255,12 @@ main(int itkNotUsed(argc), char * itkNotUsed(argv)[])
   resample->SetOutputDirection(fixedImage->GetDirection());
   resample->SetDefaultPixelValue(100);
 
+  using WriterType = itk::ImageFileWriter<ImageType>;
+
   WriterType::Pointer writer = WriterType::New();
-  writer->SetFileName("output.png");
+  writer->SetFileName(argv[3]);
   writer->SetInput(resample->GetOutput());
   writer->Update();
 
   return EXIT_SUCCESS;
 }
-
-
-void
-CreateEllipseImage(ImageType::Pointer image)
-{
-  using EllipseType = itk::EllipseSpatialObject<Dimension>;
-
-  using SpatialObjectToImageFilterType = itk::SpatialObjectToImageFilter<EllipseType, ImageType>;
-
-  SpatialObjectToImageFilterType::Pointer imageFilter = SpatialObjectToImageFilterType::New();
-
-  ImageType::SizeType size;
-  size[0] = 100;
-  size[1] = 100;
-
-  imageFilter->SetSize(size);
-
-  ImageType::SpacingType spacing;
-  spacing.Fill(1);
-  imageFilter->SetSpacing(spacing);
-
-  EllipseType::Pointer   ellipse = EllipseType::New();
-  EllipseType::ArrayType radiusArray;
-  radiusArray[0] = 10;
-  radiusArray[1] = 20;
-  ellipse->SetRadiusInObjectSpace(radiusArray);
-
-  using TransformType = EllipseType::TransformType;
-  TransformType::Pointer transform = TransformType::New();
-  transform->SetIdentity();
-
-  TransformType::OutputVectorType translation;
-  translation[0] = 65;
-  translation[1] = 45;
-  transform->Translate(translation, false);
-
-  ellipse->SetObjectToParentTransform(transform);
-
-  imageFilter->SetInput(ellipse);
-
-  ellipse->SetDefaultInsideValue(255);
-  ellipse->SetDefaultOutsideValue(0);
-  imageFilter->SetUseObjectValue(true);
-  imageFilter->SetOutsideValue(0);
-
-  imageFilter->Update();
-
-  image->Graft(imageFilter->GetOutput());
-}
-
-void
-CreateCircleImage(ImageType::Pointer image)
-{
-  using EllipseType = itk::EllipseSpatialObject<Dimension>;
-
-  using SpatialObjectToImageFilterType = itk::SpatialObjectToImageFilter<EllipseType, ImageType>;
-
-  SpatialObjectToImageFilterType::Pointer imageFilter = SpatialObjectToImageFilterType::New();
-
-  ImageType::SizeType size;
-  size[0] = 100;
-  size[1] = 100;
-
-  imageFilter->SetSize(size);
-
-  ImageType::SpacingType spacing;
-  spacing.Fill(1);
-  imageFilter->SetSpacing(spacing);
-
-  EllipseType::Pointer   ellipse = EllipseType::New();
-  EllipseType::ArrayType radiusArray;
-  radiusArray[0] = 10;
-  radiusArray[1] = 10;
-  ellipse->SetRadiusInObjectSpace(radiusArray);
-
-  using TransformType = EllipseType::TransformType;
-  TransformType::Pointer transform = TransformType::New();
-  transform->SetIdentity();
-
-  TransformType::OutputVectorType translation;
-  translation[0] = 50;
-  translation[1] = 50;
-  transform->Translate(translation, false);
-
-  ellipse->SetObjectToParentTransform(transform);
-
-  imageFilter->SetInput(ellipse);
-
-  ellipse->SetDefaultInsideValue(255);
-  ellipse->SetDefaultOutsideValue(0);
-  imageFilter->SetUseObjectValue(true);
-  imageFilter->SetOutsideValue(0);
-
-  imageFilter->Update();
-
-  image->Graft(imageFilter->GetOutput());
-}

src/Core/Transform/MutualInformationAffine/Documentation.rst

@@ -4,6 +4,7 @@ Mutual Information Affine
 .. index::
    single: AffineTransform
    single: MutualInformationImageToImageMetric
+   single: GradientDescentOptimizer
 
 Synopsis
 --------
@@ -25,20 +26,21 @@ Results
 
   moving.png
 
-.. figure:: output.png
+.. figure:: OutputBaseline.png
   :scale: 70%
-  :alt: output.png
+  :alt: OutputBaseline.png
 
-  output.png
+  OutputBaseline.png
 
 Output::
 
-  Optimizer stop condition: GradientDescentOptimizer: Maximum number of iterations (1000) exceeded.
-  Final Parameters: [2.03502150081667, 0.4815329593844025, -0.524998748481682, 0.5494736862921054, 0.016827072908497744, -0.018211282766070345]
+  Optimizer stop condition: GradientDescentOptimizer: Maximum number of iterations (200) exceeded.
+  Final Parameters: [1.0028069041777101, -0.009647107048100798, -0.010717116855425006, 1.0029040091646872, 13.26279335943067, 12.226979744206531]
+
   Result =
-   Iterations    = 1000
-   Metric value  = -0.000162484
-   Numb. Samples = 100
+   Iterations    = 200
+   Metric value  = 0.000971698
+   Numb. Samples = 567
 
 Code
 ----
@@ -49,8 +51,16 @@ C++
 .. literalinclude:: Code.cxx
    :lines: 18-
 
+
+Binder
+------
+.. image:: https://mybinder.org/badge_logo.svg
+ :target: https://mybinder.org/v2/gh/InsightSoftwareConsortium/ITKExamples/master?filepath=src%2FCore%2FTransform%2FMutualInformationAffine%2FMutualInformationAffine.ipynb
+
+
 Classes demonstrated
 --------------------
 
 .. breathelink:: itk::AffineTransform
 .. breathelink:: itk::MutualInformationImageToImageMetric
+.. breathelink:: itk::GradientDescentOptimizer