WIP: ENH: Add BlockMatching registration Python wrappings + example

examples/itkBlockMatchingMultiResolutionImageRegistrationMethodTest.py

@@ -0,0 +1,91 @@
+#!/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.
+
+import os
+from urllib.request import urlretrieve
+
+import itk
+
+# Verify wrapping succeeded
+assert 'MultiResolutionFixedBlockRadiusCalculator' in dir(itk.Ultrasound)
+
+# Get image input for registration
+FIXED_IMAGE_PATH = 'Input/rf_pre15.mha'
+MOVING_IMAGE_PATH = 'Input/rf_post15.mha'
+
+if not os.path.exists(FIXED_IMAGE_PATH):
+    url = 'https://data.kitware.com/api/v1/file/58ee3b778d777f16d095fd8a/download'
+    urlretrieve(url, FIXED_IMAGE_PATH)
+
+if not os.path.exists(MOVING_IMAGE_PATH):
+    url = 'https://data.kitware.com/api/v1/file/58ee3b758d777f16d095fd87/download'
+    urlretrieve(url, MOVING_IMAGE_PATH)
+
+fixed_image = itk.imread(FIXED_IMAGE_PATH, itk.F)
+moving_image = itk.imread(MOVING_IMAGE_PATH, itk.F)
+
+dimension = fixed_image.GetImageDimension()
+displacement_image_type = itk.Image[itk.Vector[itk.F,dimension],dimension]
+
+block_radius_calculator = itk.Ultrasound.MultiResolutionFixedBlockRadiusCalculator[type(fixed_image)].New()
+block_radius_calculator.SetRadius([12,4])
+
+# Create schedule for iterative registration
+search_region_source = itk.Ultrasound.MultiResolutionFixedSearchRegionImageSource[type(fixed_image), 
+                                                                                  type(fixed_image),
+                                                                                  displacement_image_type].New()
+pyramid_schedule = itk.Array2D[itk.UI]()
+pyramid_schedule.SetSize(3,2)
+pyramid_schedule.SetElement(0,0,3)
+pyramid_schedule.SetElement(0,1,2)
+pyramid_schedule.SetElement(1,0,2)
+pyramid_schedule.SetElement(1,1,1)
+pyramid_schedule.SetElement(2,0,1)
+pyramid_schedule.SetElement(2,1,1)
+
+search_region_source.SetPyramidSchedule(pyramid_schedule)
+search_region_source.SetSearchRegionRadiusSchedule([50,6])
+search_region_source.SetOverlapSchedule(1.0)
+
+metric_image_filter = itk.Ultrasound.NormalizedCrossCorrelationNeighborhoodIteratorMetricImageFilter[type(fixed_image),
+                                                                                                     type(fixed_image),
+                                                                                                     type(fixed_image)].New()
+
+level_registration_method = itk.Ultrasound.BlockMatchingImageRegistrationMethod[type(fixed_image),
+                                                                   type(fixed_image),
+                                                                   type(fixed_image),
+                                                                   displacement_image_type,
+                                                                   itk.F].New()
+level_registration_method.SetMetricImageFilter(metric_image_filter)
+
+# Set up the multi-resolution registration object
+multi_res_registration_method = itk.Ultrasound.BlockMatchingMultiResolutionImageRegistrationMethod[type(fixed_image),
+                                                                   type(fixed_image),
+                                                                   type(fixed_image),
+                                                                   displacement_image_type,
+                                                                   itk.F].New()
+multi_res_registration_method.SetFixedImage(fixed_image)
+multi_res_registration_method.SetMovingImage(moving_image)
+multi_res_registration_method.SetBlockRadiusCalculator(block_radius_calculator)
+multi_res_registration_method.SetSearchRegionImageSource(search_region_source)
+multi_res_registration_method.SetSchedules(pyramid_schedule, pyramid_schedule)
+multi_res_registration_method.SetImageRegistrationMethod(level_registration_method)
+
+# Run the actual registration
+multi_res_registration_method.Update()
+
+# Write out results
+itk.imwrite(multi_res_registration_method.GetOutput(), 'Output/rf_post15_registered.mha')

include/itkBlockMatchingBayesianRegularizationDisplacementCalculator.hxx

@@ -342,8 +342,8 @@ BayesianRegularizationDisplacementCalculator<TMetricImage, TDisplacementImage>::
   typename FaceCalculatorType::FaceListType::iterator fit;
 
   MetricImageImagePointerType tempMetricImageImagePtr;
-  MetricImagePointerType      postImage;
-  MetricImagePointerType      priorImage;
+  //MetricImagePointerType      postImage;
+  //MetricImagePointerType      priorImage;
   m_CurrentIteration = 0;
   // Hack to make the SplitRequestedRegion method operate on the face list
   // regions;
@@ -425,14 +425,14 @@ BayesianRegularizationDisplacementCalculator<TMetricImage, TDisplacementImage>::
     {
       priorImage = priorImageImageIt.GetPrevious(direction);
       // If we are inside the boundary.
-      if (priorImage.GetPointer() != nullptr)
+      if (priorImage != nullptr)
       {
         shift.Fill(0.0);
         shift[direction] = -1 * spacing[direction];
         this->ImpartLikelihood(postImage, priorImage, direction, shift);
       }
       priorImage = priorImageImageIt.GetNext(direction);
-      if (priorImage.GetPointer() != nullptr)
+      if (priorImage != nullptr)
       {
         shift.Fill(0.0);
         shift[direction] = spacing[direction];

include/itkBlockMatchingDisplacementPipeline.h

@@ -129,7 +129,7 @@ class ITK_TEMPLATE_EXPORT DisplacementPipeline
 
   /** The registration method. */
   using LevelRegistrationMethodType = BlockMatching::
-    ImageRegistrationMethod<FixedImageType, MovingImageType, MetricImageType, DisplacementImageType, CoordRepType>;
+    BlockMatchingImageRegistrationMethod<FixedImageType, MovingImageType, MetricImageType, DisplacementImageType, CoordRepType>;
 
   /** Interpolation classes. */
   using ParabolicInterpolatorType =
@@ -174,7 +174,7 @@ class ITK_TEMPLATE_EXPORT DisplacementPipeline
     BlockMatching::BayesianRegularizationDisplacementCalculator<MetricImageType, DisplacementImageType>;
 
   /** Multi-resolution registration method. */
-  using RegistrationMethodType = BlockMatching::MultiResolutionImageRegistrationMethod<FixedImageType,
+  using RegistrationMethodType = BlockMatching::BlockMatchingMultiResolutionImageRegistrationMethod<FixedImageType,
                                                                                        MovingImageType,
                                                                                        MetricImageType,
                                                                                        DisplacementImageType,

include/itkBlockMatchingImageRegistrationMethod.h

@@ -28,7 +28,7 @@ namespace itk
 namespace BlockMatching
 {
 
-/** \class ImageRegistrationMethod
+/** \class BlockMatchingImageRegistrationMethod
  *
  * This class defines the interface to perform deformable image registration by
  * block matching.
@@ -43,7 +43,7 @@ namespace BlockMatching
  * area.  The information from the search region image (origin, spacing, region,
  * etc) determines the information in the output displacement image.
  *
- * \sa ImageRegistrationMethod
+ * \sa BlockMatchingImageRegistrationMethod
  *
  * \ingroup RegistrationFilters
  * \ingroup Ultrasound
@@ -53,12 +53,12 @@ template <typename TFixedImage,
           typename TMetricImage,
           typename TDisplacementImage,
           typename TCoordRep>
-class ITK_TEMPLATE_EXPORT ImageRegistrationMethod
+class ITK_TEMPLATE_EXPORT BlockMatchingImageRegistrationMethod
   : public ImageToImageFilter<itk::Image<typename TMovingImage::RegionType, TDisplacementImage::ImageDimension>,
                               TDisplacementImage>
 {
 public:
-  ITK_DISALLOW_COPY_AND_ASSIGN(ImageRegistrationMethod);
+  ITK_DISALLOW_COPY_AND_ASSIGN(BlockMatchingImageRegistrationMethod);
 
   /** ImageDimension enumeration. */
   itkStaticConstMacro(ImageDimension, unsigned int, TDisplacementImage::ImageDimension);
@@ -92,7 +92,7 @@ class ITK_TEMPLATE_EXPORT ImageRegistrationMethod
   using SearchRegionImageType = Image<typename MovingImageType::RegionType, ImageDimension>;
 
   /** Standard class type alias. */
-  using Self = ImageRegistrationMethod;
+  using Self = BlockMatchingImageRegistrationMethod;
   using Superclass = ImageToImageFilter<SearchRegionImageType, TDisplacementImage>;
   using Pointer = SmartPointer<Self>;
   using ConstPointer = SmartPointer<const Self>;
@@ -101,7 +101,7 @@ class ITK_TEMPLATE_EXPORT ImageRegistrationMethod
   itkNewMacro(Self);
 
   /** Run-time type information (and related methods). */
-  itkTypeMacro(ImageRegistrationMethod, ImageSource);
+  itkTypeMacro(BlockMatchingImageRegistrationMethod, ImageSource);
 
   /** Type of the point use for determing the location in the fixed image of a
    * block's center. */
@@ -188,8 +188,8 @@ class ITK_TEMPLATE_EXPORT ImageRegistrationMethod
 
 
 protected:
-  ImageRegistrationMethod();
-  virtual ~ImageRegistrationMethod() {}
+  BlockMatchingImageRegistrationMethod();
+  virtual ~BlockMatchingImageRegistrationMethod() {}
 
   void
   GenerateOutputInformation() override;

include/itkBlockMatchingImageRegistrationMethod.hxx

@@ -36,8 +36,8 @@ template <typename TFixedImage,
           typename TMetricImage,
           typename TDisplacementImage,
           typename TCoordRep>
-ImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::
-  ImageRegistrationMethod()
+BlockMatchingImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::
+  BlockMatchingImageRegistrationMethod()
   : m_UseStreaming(false)
 {
   m_FixedImage = nullptr;
@@ -55,7 +55,7 @@ template <typename TFixedImage,
           typename TDisplacementImage,
           typename TCoordRep>
 void
-ImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::SetFixedImage(
+BlockMatchingImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::SetFixedImage(
   FixedImageType * fixedImage)
 {
   if (this->m_FixedImage.GetPointer() != fixedImage)
@@ -72,7 +72,7 @@ template <typename TFixedImage,
           typename TDisplacementImage,
           typename TCoordRep>
 void
-ImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::SetMovingImage(
+BlockMatchingImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::SetMovingImage(
   MovingImageType * movingImage)
 {
   if (this->m_MovingImage.GetPointer() != movingImage)
@@ -89,7 +89,7 @@ template <typename TFixedImage,
           typename TDisplacementImage,
           typename TCoordRep>
 void
-ImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::
+BlockMatchingImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::
   GenerateOutputInformation()
 {
   Superclass::GenerateOutputInformation();
@@ -109,7 +109,7 @@ template <typename TFixedImage,
           typename TDisplacementImage,
           typename TCoordRep>
 void
-ImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::
+BlockMatchingImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::
   GenerateInputRequestedRegion()
 {
   Superclass::GenerateInputRequestedRegion();
@@ -132,7 +132,7 @@ template <typename TFixedImage,
           typename TDisplacementImage,
           typename TCoordRep>
 void
-ImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::
+BlockMatchingImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::
   EnlargeOutputRequestedRegion(DataObject * data)
 {
   this->m_MetricImageToDisplacementCalculator->ModifyEnlargeOutputRequestedRegion(data);
@@ -145,7 +145,7 @@ template <typename TFixedImage,
           typename TDisplacementImage,
           typename TCoordRep>
 void
-ImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::GenerateData()
+BlockMatchingImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::GenerateData()
 {
   const SearchRegionImageType * input = this->GetInput();
 
@@ -212,7 +212,7 @@ template <typename TFixedImage,
           typename TDisplacementImage,
           typename TCoordRep>
 void
-ImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::Initialize()
+BlockMatchingImageRegistrationMethod<TFixedImage, TMovingImage, TMetricImage, TDisplacementImage, TCoordRep>::Initialize()
 {
   if (!m_FixedImage)
   {

include/itkBlockMatchingMetricImageFilter.h

@@ -31,7 +31,7 @@ namespace BlockMatching
  * region.
  *
  * This class is intended to be used internally by
- * itk::BlockMatching::ImageRegistrationMethod and its ilk.
+ * itk::BlockMatching::BlockMatchingImageRegistrationMethod and its ilk.
  *
  * There are two inputs to this image, the Fixed Image, and the Moving Image.
  * A block or kernel from the FixedImage is specified with SetFixedRegion().

include/itkBlockMatchingMetricImageToDisplacementCalculator.h

@@ -32,7 +32,7 @@ namespace BlockMatching
  *
  * \brief Calculates the displacement of a block from the MetricImage(s).
  *
- * This class tightly integrates with BlockMatching::ImageRegistrationMethod.
+ * This class tightly integrates with BlockMatching::BlockMatchingImageRegistrationMethod.
  * It is the responsibility of this class to calculate the displacement given
  * MetricImage(s) in a block matching registration method.
  *
@@ -47,7 +47,7 @@ namespace BlockMatching
  *
  * Caching of the MetricImage can be enabled by SetCacheMetricImageOn();
  *
- * The behavior of the associated BlockMatching::ImageRegistrationMethod
+ * The behavior of the associated BlockMatching::BlockMatchingImageRegistrationMethod
  * GenerateInputRequestedRegion() and EnlargeOutputRequestedRegion() with
  * ModifyGenerateInputRequestedRegion() and
  * ModifyEnlargeOutputRequestedRegion().
@@ -71,7 +71,7 @@ class ITK_TEMPLATE_EXPORT MetricImageToDisplacementCalculator : public Object
 
   /** Type of the metric image (input pixels). */
   using MetricImageType = TMetricImage;
-  using MetricImagePointerType = typename MetricImageType::Pointer;
+  using MetricImagePointerType = MetricImageType*; // TODO ::Pointer
   using IndexType = typename MetricImageType::IndexType;
 
   /** Type of the displacement image (output). */
@@ -134,13 +134,13 @@ class ITK_TEMPLATE_EXPORT MetricImageToDisplacementCalculator : public Object
   virtual void
   Compute() = 0;
 
-  /** Modify the associated BlockMatching::ImageRegistrationMethod's
+  /** Modify the associated BlockMatching::BlockMatchingImageRegistrationMethod's
    * GenerateInputRequestedRegion().  */
   virtual void
   ModifyGenerateInputRequestedRegion(RegionType & region)
   {}
 
-  /** Modify the associated BlockMatching::ImageRegistrationMethod's
+  /** Modify the associated BlockMatching::BlockMatchingImageRegistrationMethod's
    * EnlargeOutputRequestedRegion().  */
   virtual void
   ModifyEnlargeOutputRequestedRegion(DataObject * data)

include/itkBlockMatchingMultiResolutionBlockRadiusCalculator.h

@@ -29,7 +29,7 @@ namespace BlockMatching
 /** \class MultiResolutionBlockRadiusCalculator
  *
  * \brief Base class calculate the fixed image matching kernel radius for each level in a
- * BlockMatching::MultiResolutionImageRegistrationMethod.
+ * BlockMatching::BlockMatchingMultiResolutionImageRegistrationMethod.
  *
  * This must be able to produce the fixed image block radius for every level of
  * the MultiResolutionPyramidImage filter.
@@ -74,7 +74,7 @@ class ITK_TEMPLATE_EXPORT MultiResolutionBlockRadiusCalculator : public Object
   itkGetConstReferenceMacro(PyramidSchedule, ScheduleType);
 
   /** Set/Get the Fixed image. The fixed image is set during the
-   * BlockMatching::MultiResolutionImageRegistrationMethod and is available for
+   * BlockMatching::BlockMatchingMultiResolutionImageRegistrationMethod and is available for
    * child classes if the choose to use it.  */
   itkSetConstObjectMacro(FixedImage, FixedImageType);
   itkGetConstObjectMacro(FixedImage, FixedImageType);

include/itkBlockMatchingMultiResolutionFixedBlockRadiusCalculator.h

@@ -30,7 +30,7 @@ namespace BlockMatching
  * \brief A fixed radius is used for every level.
  *
  * This class generates the fixed image matching kernel radius in a
- * BlockMatching::MultiResolutionImageRegistrationMethod.  A fixed block radius
+ * BlockMatching::BlockMatchingMultiResolutionImageRegistrationMethod.  A fixed block radius
  * is used, which means the size of the block in physical coordinates then
  * scales with the pyramid schedule.
  *