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; }
    SITK_RETURN_SELF_TYPE_HEADER SetInterpolator ( InterpolatorEnum Interpolator )
      { this->m_Interpolator = Interpolator; return *this; }
#if !defined(SWIG) || defined(JAVASWIG) || defined(CSHARPSWIG)
    // Only wrap this method if the wrapping language is strongly typed
    SITK_RETURN_SELF_TYPE_HEADER SetInitialTransform ( const Transform &transform );
#endif
    SITK_RETURN_SELF_TYPE_HEADER 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>() );
 
    SITK_RETURN_SELF_TYPE_HEADER SetMovingInitialTransform( const Transform &transform )
    { this->m_MovingInitialTransform = transform; return *this; }
    Transform GetMovingInitialTransform( ) const
    { return this->m_MovingInitialTransform; }
    SITK_RETURN_SELF_TYPE_HEADER SetFixedInitialTransform( const Transform &transform )
    { this->m_FixedInitialTransform = transform; return *this; }
    Transform GetFixedInitialTransform( ) const
    { return this->m_FixedInitialTransform; }
    SITK_RETURN_SELF_TYPE_HEADER SetVirtualDomain( const std::vector<uint32_t> &virtualSize,
                            const std::vector<double> &virtualOrigin,
                            const std::vector<double> &virtualSpacing,
                            const std::vector<double> &virtualDirection );
    SITK_RETURN_SELF_TYPE_HEADER SetVirtualDomainFromImage( const Image &virtualImage );
    SITK_RETURN_SELF_TYPE_HEADER SetMetricAsANTSNeighborhoodCorrelation( unsigned int radius );
 
    SITK_RETURN_SELF_TYPE_HEADER SetMetricAsCorrelation( );
 
    SITK_RETURN_SELF_TYPE_HEADER SetMetricAsDemons( double intensityDifferenceThreshold = 0.001 );
 
    SITK_RETURN_SELF_TYPE_HEADER SetMetricAsJointHistogramMutualInformation( unsigned int numberOfHistogramBins = 20,
                                                      double varianceForJointPDFSmoothing = 1.5);
 
    SITK_RETURN_SELF_TYPE_HEADER SetMetricAsMeanSquares( );
 
    SITK_RETURN_SELF_TYPE_HEADER SetMetricAsMattesMutualInformation( unsigned int numberOfHistogramBins = 50 );
 
 
    enum EstimateLearningRateType { Never, 
                                    Once,  
                                    EachIteration 
    };
 
    SITK_RETURN_SELF_TYPE_HEADER 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 );
 
    SITK_RETURN_SELF_TYPE_HEADER SetOptimizerAsRegularStepGradientDescent( double learningRate,
                                                    double minStep,
                                                    unsigned int numberOfIterations,
                                                    double relaxationFactor=0.5,
                                                    double gradientMagnitudeTolerance = 1e-4,
                                                    EstimateLearningRateType estimateLearningRate = Never,
                                                    double maximumStepSizeInPhysicalUnits = 0.0);
 
    SITK_RETURN_SELF_TYPE_HEADER SetOptimizerAsGradientDescent( double learningRate,
                                         unsigned int numberOfIterations,
                                         double convergenceMinimumValue = 1e-6,
                                         unsigned int convergenceWindowSize = 10,
                                         EstimateLearningRateType estimateLearningRate = Once,
                                         double maximumStepSizeInPhysicalUnits = 0.0);
 
    SITK_RETURN_SELF_TYPE_HEADER 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 );
 
    SITK_RETURN_SELF_TYPE_HEADER 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 );
 
    SITK_RETURN_SELF_TYPE_HEADER 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);
 
 
    SITK_RETURN_SELF_TYPE_HEADER SetOptimizerAsExhaustive( const std::vector<unsigned int> &numberOfSteps,
                                    double stepLength = 1.0 );
 
    SITK_RETURN_SELF_TYPE_HEADER SetOptimizerAsAmoeba(double simplexDelta,
                               unsigned int numberOfIterations,
                               double parametersConvergenceTolerance=1e-8,
                               double functionConvergenceTolerance=1e-4,
                               bool withRestarts = false );
 
    SITK_RETURN_SELF_TYPE_HEADER SetOptimizerWeights(const std::vector<double> &weights);
    std::vector<double> GetOptimizerWeights( ) const;
    SITK_RETURN_SELF_TYPE_HEADER SetOptimizerAsPowell(unsigned int numberOfIterations = 100,
                               unsigned int maximumLineIterations = 100,
                               double stepLength = 1,
                               double stepTolerance = 1e-6,
                               double valueTolerance = 1e-6 );
 
    SITK_RETURN_SELF_TYPE_HEADER 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);
 
 
 
    SITK_RETURN_SELF_TYPE_HEADER SetOptimizerScales( const std::vector<double> &scales );
 
    SITK_RETURN_SELF_TYPE_HEADER SetOptimizerScalesFromJacobian( unsigned int centralRegionRadius = 5 );
 
    SITK_RETURN_SELF_TYPE_HEADER SetOptimizerScalesFromIndexShift( unsigned int centralRegionRadius = 5,
                                           double smallParameterVariation =  0.01 );
 
    SITK_RETURN_SELF_TYPE_HEADER SetOptimizerScalesFromPhysicalShift( unsigned int centralRegionRadius = 5,
                                              double smallParameterVariation =  0.01 );
 
 
    SITK_RETURN_SELF_TYPE_HEADER SetMetricFixedMask( const Image &binaryMask );
 
    SITK_RETURN_SELF_TYPE_HEADER SetMetricMovingMask( const Image &binaryMask );
 
    SITK_RETURN_SELF_TYPE_HEADER SetMetricSamplingPercentage(double percentage, unsigned int seed = sitkWallClock);
    SITK_RETURN_SELF_TYPE_HEADER SetMetricSamplingPercentagePerLevel(const std::vector<double> &percentage, unsigned int seed = sitkWallClock);
    const std::vector<double> &GetMetricSamplingPercentagePerLevel() const;
 
    enum MetricSamplingStrategyType {
      NONE,
      REGULAR,
      RANDOM
    };
 
    SITK_RETURN_SELF_TYPE_HEADER SetMetricSamplingStrategy( MetricSamplingStrategyType strategy);
 
    SITK_RETURN_SELF_TYPE_HEADER SetMetricUseFixedImageGradientFilter(bool);
    SITK_RETURN_SELF_TYPE_HEADER MetricUseFixedImageGradientFilterOn() {return this->SetMetricUseFixedImageGradientFilter(true);}
    SITK_RETURN_SELF_TYPE_HEADER MetricUseFixedImageGradientFilterOff() {return this->SetMetricUseFixedImageGradientFilter(false);}
    SITK_RETURN_SELF_TYPE_HEADER SetMetricUseMovingImageGradientFilter(bool);
    SITK_RETURN_SELF_TYPE_HEADER MetricUseMovingImageGradientFilterOn() {return this->SetMetricUseMovingImageGradientFilter(true);}
    SITK_RETURN_SELF_TYPE_HEADER MetricUseMovingImageGradientFilterOff() {return this->SetMetricUseMovingImageGradientFilter(false);}
    SITK_RETURN_SELF_TYPE_HEADER SetShrinkFactorsPerLevel( const std::vector<unsigned int> &shrinkFactors );
 
    SITK_RETURN_SELF_TYPE_HEADER SetSmoothingSigmasPerLevel( const std::vector<double> &smoothingSigmas );
 
    SITK_RETURN_SELF_TYPE_HEADER SetSmoothingSigmasAreSpecifiedInPhysicalUnits(bool arg);
    SITK_RETURN_SELF_TYPE_HEADER SmoothingSigmasAreSpecifiedInPhysicalUnitsOn() { this->SetSmoothingSigmasAreSpecifiedInPhysicalUnits(true); return *this;}
    SITK_RETURN_SELF_TYPE_HEADER SmoothingSigmasAreSpecifiedInPhysicalUnitsOff()  { this->SetSmoothingSigmasAreSpecifiedInPhysicalUnits(false); return *this;}
    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 >
      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>;
    std::unique_ptr<detail::MemberFunctionFactory<MemberFunctionType> > m_MemberFactory;
    std::unique_ptr<detail::MemberFunctionFactory<EvaluateMemberFunctionType> > m_EvaluateMemberFactory;
 
    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;
  };
 
}
}
 
#endif // sitkImageRegistrationMethod_h