sitkImageRegistrationMethod.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  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.
 *
 *=========================================================================*/
#ifndef sitkImageRegistrationMethod_h
#define sitkImageRegistrationMethod_h
 
#include "sitkRegistration.h"
 
#include "sitkDetail.h"
#include "sitkImage.h"
#include "sitkPixelIDTokens.h"
#include "sitkMemberFunctionFactory.h"
#include "sitkProcessObject.h"
 
#include "sitkRandomSeed.h"
#include "sitkInterpolator.h"
#include "sitkTransform.h"
 
namespace itk
{
 
#ifndef SWIG
template <typename TInternalComputationValueType>
class ObjectToObjectOptimizerBaseTemplate;
template <typename TFixedImage,
          typename TMovingImage,
          typename TVirtualImage,
          typename TInternalComputationValueType,
          typename TMetricTraits>
class ImageToImageMetricv4;
 
template <typename TFixedImageType, typename TMovingImageType, typename TVirtualImageType, typename TCoordRep>
class DefaultImageToImageMetricTraitsv4;
 
template <typename TMetric>
class RegistrationParameterScalesEstimator;
 
template <unsigned int VDimension>
class SpatialObject;
 
class Command;
class EventObject;
#endif
 
 
namespace simple
{
class BSplineTransform;

class SITKRegistration_EXPORT ImageRegistrationMethod : public ProcessObject
{
public:
  using Self = ImageRegistrationMethod;
  using Superclass = ProcessObject;
 
  ~ImageRegistrationMethod() override;
 
  ImageRegistrationMethod();
 
  std::string
  GetName() const override
  {
    return std::string("ImageRegistrationMethod");
  }

  std::string
  ToString() const override;

  InterpolatorEnum
  GetInterpolator()
  {
    return this->m_Interpolator;
  }
  void
  SetInterpolator(InterpolatorEnum Interpolator)
  {
    this->m_Interpolator = Interpolator;
  }


#if !defined(SWIG) || defined(JAVASWIG) || defined(CSHARPSWIG)
  // Only wrap this method if the wrapping language is strongly typed
  void
  SetInitialTransform(const Transform & transform);
#endif
  void
  SetInitialTransform(Transform & transform, bool inPlace = true);
  Transform
  GetInitialTransform()
  {
    return this->m_InitialTransform;
  }
  bool
  GetInitialTransformInPlace() const
  {
    return this->m_InitialTransformInPlace;
  }


  void
  SetInitialTransformAsBSpline(BSplineTransform &                transform,
                               bool                              inPlace = true,
                               const std::vector<unsigned int> & scaleFactors = std::vector<unsigned int>());

  void
  SetMovingInitialTransform(const Transform & transform)
  {
    this->m_MovingInitialTransform = transform;
  }
  Transform
  GetMovingInitialTransform() const
  {
    return this->m_MovingInitialTransform;
  }


  void
  SetFixedInitialTransform(const Transform & transform)
  {
    this->m_FixedInitialTransform = transform;
  }
  Transform
  GetFixedInitialTransform() const
  {
    return this->m_FixedInitialTransform;
  }

 

  void
  SetVirtualDomain(const std::vector<uint32_t> & virtualSize,
                   const std::vector<double> &   virtualOrigin,
                   const std::vector<double> &   virtualSpacing,
                   const std::vector<double> &   virtualDirection);
  void
  SetVirtualDomainFromImage(const Image & virtualImage);

  void
  SetMetricAsANTSNeighborhoodCorrelation(unsigned int radius);

  void
  SetMetricAsCorrelation();

  void
  SetMetricAsDemons(double intensityDifferenceThreshold = 0.001);

  void
  SetMetricAsJointHistogramMutualInformation(unsigned int numberOfHistogramBins = 20,
                                             double       varianceForJointPDFSmoothing = 1.5);

  void
  SetMetricAsMeanSquares();

  void
  SetMetricAsMattesMutualInformation(unsigned int numberOfHistogramBins = 50);
 
 
  enum EstimateLearningRateType
  {
    Never,        
    Once,         
    EachIteration 
  };

  void
  SetOptimizerAsConjugateGradientLineSearch(double                   learningRate,
                                            unsigned int             numberOfIterations,
                                            double                   convergenceMinimumValue = 1e-6,
                                            unsigned int             convergenceWindowSize = 10,
                                            double                   lineSearchLowerLimit = 0,
                                            double                   lineSearchUpperLimit = 5.0,
                                            double                   lineSearchEpsilon = 0.01,
                                            unsigned int             lineSearchMaximumIterations = 20,
                                            EstimateLearningRateType estimateLearningRate = Once,
                                            double                   maximumStepSizeInPhysicalUnits = 0.0);

  void
  SetOptimizerAsRegularStepGradientDescent(double                   learningRate,
                                           double                   minStep,
                                           unsigned int             numberOfIterations,
                                           double                   relaxationFactor = 0.5,
                                           double                   gradientMagnitudeTolerance = 1e-4,
                                           EstimateLearningRateType estimateLearningRate = Never,
                                           double                   maximumStepSizeInPhysicalUnits = 0.0);

  void
  SetOptimizerAsGradientDescent(double                   learningRate,
                                unsigned int             numberOfIterations,
                                double                   convergenceMinimumValue = 1e-6,
                                unsigned int             convergenceWindowSize = 10,
                                EstimateLearningRateType estimateLearningRate = Once,
                                double                   maximumStepSizeInPhysicalUnits = 0.0);

  void
  SetOptimizerAsGradientDescentLineSearch(double                   learningRate,
                                          unsigned int             numberOfIterations,
                                          double                   convergenceMinimumValue = 1e-6,
                                          unsigned int             convergenceWindowSize = 10,
                                          double                   lineSearchLowerLimit = 0,
                                          double                   lineSearchUpperLimit = 5.0,
                                          double                   lineSearchEpsilon = 0.01,
                                          unsigned int             lineSearchMaximumIterations = 20,
                                          EstimateLearningRateType estimateLearningRate = Once,
                                          double                   maximumStepSizeInPhysicalUnits = 0.0);

  void
  SetOptimizerAsLBFGSB(double       gradientConvergenceTolerance = 1e-5,
                       unsigned int numberOfIterations = 500,
                       unsigned int maximumNumberOfCorrections = 5,
                       unsigned int maximumNumberOfFunctionEvaluations = 2000,
                       double       costFunctionConvergenceFactor = 1e+7,
                       double       lowerBound = std::numeric_limits<double>::min(),
                       double       upperBound = std::numeric_limits<double>::max(),
                       bool         trace = false);

  void
  SetOptimizerAsLBFGS2(double       solutionAccuracy = 1e-5,
                       unsigned int numberOfIterations = 0,
                       unsigned int hessianApproximateAccuracy = 6,
                       unsigned int deltaConvergenceDistance = 0,
                       double       deltaConvergenceTolerance = 1e-5,
                       unsigned int lineSearchMaximumEvaluations = 40,
                       double       lineSearchMinimumStep = 1e-20,
                       double       lineSearchMaximumStep = 1e20,
                       double       lineSearchAccuracy = 1e-4);
 

  void
  SetOptimizerAsExhaustive(const std::vector<unsigned int> & numberOfSteps, double stepLength = 1.0);

  void
  SetOptimizerAsAmoeba(double       simplexDelta,
                       unsigned int numberOfIterations,
                       double       parametersConvergenceTolerance = 1e-8,
                       double       functionConvergenceTolerance = 1e-4,
                       bool         withRestarts = false);

  void
  SetOptimizerWeights(const std::vector<double> & weights);
  std::vector<double>
  GetOptimizerWeights() const;

  void
  SetOptimizerAsPowell(unsigned int numberOfIterations = 100,
                       unsigned int maximumLineIterations = 100,
                       double       stepLength = 1,
                       double       stepTolerance = 1e-6,
                       double       valueTolerance = 1e-6);

  void
  SetOptimizerAsOnePlusOneEvolutionary(unsigned int numberOfIterations = 100,
                                       double       epsilon = 1.5e-4,
                                       double       initialRadius = 1.01,
                                       double       growthFactor = -1.0,
                                       double       shrinkFactor = -1.0,
                                       unsigned int seed = sitkWallClock);
 

  void
  SetOptimizerScales(const std::vector<double> & scales);

  void
  SetOptimizerScalesFromJacobian(unsigned int centralRegionRadius = 5);

  void
  SetOptimizerScalesFromIndexShift(unsigned int centralRegionRadius = 5, double smallParameterVariation = 0.01);

  void
  SetOptimizerScalesFromPhysicalShift(unsigned int centralRegionRadius = 5, double smallParameterVariation = 0.01);
 

  void
  SetMetricFixedMask(const Image & binaryMask);

  void
  SetMetricMovingMask(const Image & binaryMask);

  void
  SetMetricSamplingPercentage(double percentage, unsigned int seed = sitkWallClock);
  void
  SetMetricSamplingPercentagePerLevel(const std::vector<double> & percentage, unsigned int seed = sitkWallClock);

  const std::vector<double> &
  GetMetricSamplingPercentagePerLevel() const;

  enum MetricSamplingStrategyType
  {
    NONE,
    REGULAR,
    RANDOM
  };

  void
  SetMetricSamplingStrategy(MetricSamplingStrategyType strategy);

  void
  SetMetricUseFixedImageGradientFilter(bool);
  void
  MetricUseFixedImageGradientFilterOn()
  {
    return this->SetMetricUseFixedImageGradientFilter(true);
  }
  void
  MetricUseFixedImageGradientFilterOff()
  {
    return this->SetMetricUseFixedImageGradientFilter(false);
  }

 

  void
  SetMetricUseMovingImageGradientFilter(bool);
  void
  MetricUseMovingImageGradientFilterOn()
  {
    return this->SetMetricUseMovingImageGradientFilter(true);
  }
  void
  MetricUseMovingImageGradientFilterOff()
  {
    return this->SetMetricUseMovingImageGradientFilter(false);
  }

 

  void
  SetShrinkFactorsPerLevel(const std::vector<unsigned int> & shrinkFactors);

  void
  SetSmoothingSigmasPerLevel(const std::vector<double> & smoothingSigmas);

  void
  SetSmoothingSigmasAreSpecifiedInPhysicalUnits(bool arg);
  void
  SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()
  {
    this->SetSmoothingSigmasAreSpecifiedInPhysicalUnits(true);
  }
  void
  SmoothingSigmasAreSpecifiedInPhysicalUnitsOff()
  {
    this->SetSmoothingSigmasAreSpecifiedInPhysicalUnits(false);
  }

 

  Transform
  Execute(const Image & fixed, const Image & moving);
 

  double
  MetricEvaluate(const Image & fixed, const Image & moving);
 

  unsigned int
  GetOptimizerIteration() const;
  std::vector<double>
  GetOptimizerPosition() const;
  double
  GetOptimizerLearningRate() const;
  double
  GetOptimizerConvergenceValue() const;
  double
  GetMetricValue() const;

  uint64_t
  GetMetricNumberOfValidPoints() const;
 
 
  unsigned int
  GetCurrentLevel() const;

  std::vector<double>
  GetOptimizerScales() const;

  std::string
  GetOptimizerStopConditionDescription() const;
 

  bool
  StopRegistration();
 
protected:
  template <class TImage>
  Transform
  ExecuteInternal(const Image & fixed, const Image & moving);
 
  template <class TImage>
  double
  EvaluateInternal(const Image & fixed, const Image & moving);
 
 
  itk::ObjectToObjectOptimizerBaseTemplate<double> *
  CreateOptimizer(unsigned int numberOfTransformParameters);
 
  template <unsigned int VDimension>
  itk::SpatialObject<VDimension> *
  CreateSpatialObjectMask(const Image & mask);
 
  template <class TImageType>
  itk::ImageToImageMetricv4<TImageType,
                            TImageType,
                            TImageType,
                            double,
                            itk::DefaultImageToImageMetricTraitsv4<TImageType, TImageType, TImageType, double>> *
  CreateMetric();
 
  template <class TImageType>
  void
  SetupMetric(
    itk::ImageToImageMetricv4<TImageType,
                              TImageType,
                              TImageType,
                              double,
                              itk::DefaultImageToImageMetricTraitsv4<TImageType, TImageType, TImageType, double>> *,
    const TImageType *,
    const TImageType *);
 
  template <typename TMetric>
  itk::RegistrationParameterScalesEstimator<TMetric> *
  CreateScalesEstimator();
 
  template <typename TTransformAdaptorPointer, typename TRegistrationMethod>
  std::vector<TTransformAdaptorPointer>
  CreateTransformParametersAdaptor(TRegistrationMethod * method);
 
  void
  PreUpdate(itk::ProcessObject * p) override;
  void
  OnActiveProcessDelete() noexcept override;
  unsigned long
  AddITKObserver(const itk::EventObject &, itk::Command *) override;
  void
  RemoveITKObserver(EventCommand & e) override;
 
private:
  std::function<unsigned int()>        m_pfGetOptimizerIteration;
  std::function<std::vector<double>()> m_pfGetOptimizerPosition;
  std::function<double()>              m_pfGetOptimizerLearningRate;
  std::function<double()>              m_pfGetOptimizerConvergenceValue;
  std::function<double()>              m_pfGetMetricValue;
  std::function<uint64_t()>            m_pfGetMetricNumberOfValidPoints;
  std::function<std::vector<double>()> m_pfGetOptimizerScales;
  std::function<std::string()>         m_pfGetOptimizerStopConditionDescription;
  std::function<void()>                m_pfOptimizerStopRegistration;
 
 
  std::function<unsigned int()> m_pfGetCurrentLevel;
 
  std::function<void(itk::TransformBase * outTransform)> m_pfUpdateWithBestValue;
 
  template <class TMemberFunctionPointer>
  struct EvaluateMemberFunctionAddressor
  {
    using ObjectType = typename ::detail::FunctionTraits<TMemberFunctionPointer>::ClassType;
 
    template <typename TImageType>
    constexpr TMemberFunctionPointer
    operator()() const
    {
      return &ObjectType::template EvaluateInternal<TImageType>;
    }
  };
 
  typedef Transform (ImageRegistrationMethod::*MemberFunctionType)(const Image & fixed, const Image & moving);
  typedef double (ImageRegistrationMethod::*EvaluateMemberFunctionType)(const Image & fixed, const Image & moving);
  friend struct detail::MemberFunctionAddressor<MemberFunctionType>;
 
  InterpolatorEnum m_Interpolator;
  Transform        m_InitialTransform;
  bool             m_InitialTransformInPlace;
  Transform        m_MovingInitialTransform;
  Transform        m_FixedInitialTransform;
 
  std::vector<uint32_t> m_VirtualDomainSize;
  std::vector<double>   m_VirtualDomainOrigin;
  std::vector<double>   m_VirtualDomainSpacing;
  std::vector<double>   m_VirtualDomainDirection;
 
  // optimizer
  enum OptimizerType
  {
    ConjugateGradientLineSearch,
    RegularStepGradientDescent,
    GradientDescent,
    GradientDescentLineSearch,
    LBFGSB,
    Exhaustive,
    Amoeba,
    Powell,
    OnePlusOneEvolutionary,
    LBFGS2
  };
  OptimizerType             m_OptimizerType;
  double                    m_OptimizerLearningRate;
  double                    m_OptimizerMinimumStepLength;
  unsigned int              m_OptimizerNumberOfIterations;
  double                    m_OptimizerLineSearchLowerLimit;
  double                    m_OptimizerLineSearchUpperLimit;
  double                    m_OptimizerLineSearchEpsilon;
  unsigned int              m_OptimizerLineSearchMaximumIterations;
  EstimateLearningRateType  m_OptimizerEstimateLearningRate;
  double                    m_OptimizerMaximumStepSizeInPhysicalUnits;
  double                    m_OptimizerRelaxationFactor;
  double                    m_OptimizerGradientMagnitudeTolerance;
  double                    m_OptimizerConvergenceMinimumValue;
  unsigned int              m_OptimizerConvergenceWindowSize;
  double                    m_OptimizerGradientConvergenceTolerance;
  unsigned int              m_OptimizerMaximumNumberOfCorrections;
  unsigned int              m_OptimizerMaximumNumberOfFunctionEvaluations;
  double                    m_OptimizerCostFunctionConvergenceFactor;
  double                    m_OptimizerLowerBound;
  double                    m_OptimizerUpperBound;
  bool                      m_OptimizerTrace;
  std::vector<unsigned int> m_OptimizerNumberOfSteps;
  double                    m_OptimizerStepLength;
  double                    m_OptimizerSimplexDelta;
  double                    m_OptimizerParametersConvergenceTolerance;
  double                    m_OptimizerFunctionConvergenceTolerance;
  bool                      m_OptimizerWithRestarts;
  unsigned int              m_OptimizerMaximumLineIterations;
  double                    m_OptimizerStepTolerance;
  double                    m_OptimizerValueTolerance;
  double                    m_OptimizerEpsilon;
  double                    m_OptimizerInitialRadius;
  double                    m_OptimizerGrowthFactor;
  double                    m_OptimizerShrinkFactor;
  unsigned int              m_OptimizerSeed;
  double                    m_OptimizerSolutionAccuracy;
  unsigned int              m_OptimizerHessianApproximationAccuracy;
  unsigned int              m_OptimizerDeltaConvergenceDistance;
  double                    m_OptimizerDeltaConvergenceTolerance;
  unsigned int              m_OptimizerLineSearchMaximumEvaluations;
  double                    m_OptimizerLineSearchMinimumStep;
  double                    m_OptimizerLineSearchMaximumStep;
  double                    m_OptimizerLineSearchAccuracy;
 
 
  std::vector<unsigned int> m_TransformBSplineScaleFactors;
 
  std::vector<double> m_OptimizerWeights;
 
  enum OptimizerScalesType
  {
    Manual,
    Jacobian,
    IndexShift,
    PhysicalShift
  };
  OptimizerScalesType m_OptimizerScalesType;
  std::vector<double> m_OptimizerScales;
  unsigned int        m_OptimizerScalesCentralRegionRadius;
  double              m_OptimizerScalesSmallParameterVariation;
 
  // metric
  enum MetricType
  {
    ANTSNeighborhoodCorrelation,
    Correlation,
    Demons,
    JointHistogramMutualInformation,
    MeanSquares,
    MattesMutualInformation
  };
  MetricType   m_MetricType;
  unsigned int m_MetricRadius;
  double       m_MetricIntensityDifferenceThreshold;
  unsigned int m_MetricNumberOfHistogramBins;
  double       m_MetricVarianceForJointPDFSmoothing;
 
  Image m_MetricFixedMaskImage;
  Image m_MetricMovingMaskImage;
 
  std::vector<double>        m_MetricSamplingPercentage;
  MetricSamplingStrategyType m_MetricSamplingStrategy;
  unsigned int               m_MetricSamplingSeed;
 
  bool m_MetricUseFixedImageGradientFilter;
  bool m_MetricUseMovingImageGradientFilter;
 
  std::vector<unsigned int> m_ShrinkFactorsPerLevel;
  std::vector<double>       m_SmoothingSigmasPerLevel;
  bool                      m_SmoothingSigmasAreSpecifiedInPhysicalUnits;
 
  std::string  m_StopConditionDescription;
  double       m_MetricValue;
  unsigned int m_Iteration;
  uint64_t     m_NumberOfValidPoints;
 
  itk::ObjectToObjectOptimizerBaseTemplate<double> * m_ActiveOptimizer;
};
 
} // namespace simple
} // namespace itk
 
#endif // sitkImageRegistrationMethod_h