sitkImageRegistrationMethod.h

Go to the documentation of this file.

/*=========================================================================
*
*  Copyright Insight Software Consortium
*
*  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:

    typedef ImageRegistrationMethod Self;
    typedef ProcessObject Superclass;

    virtual ~ImageRegistrationMethod();

    ImageRegistrationMethod();

    std::string GetName() const { return std::string("ImageRegistrationMethod"); }

    std::string ToString() const;

    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;

  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);

    virtual void PreUpdate( itk::ProcessObject *p );
    virtual void OnActiveProcessDelete( ) SITK_NOEXCEPT;
    virtual unsigned long AddITKObserver(const itk::EventObject &, itk::Command *);
    virtual void RemoveITKObserver( EventCommand &e );

  private:

    nsstd::function<unsigned int()> m_pfGetOptimizerIteration;
    nsstd::function<std::vector<double>()> m_pfGetOptimizerPosition;
    nsstd::function<double()> m_pfGetOptimizerLearningRate;
    nsstd::function<double()> m_pfGetOptimizerConvergenceValue;
    nsstd::function<double()> m_pfGetMetricValue;
    nsstd::function<uint64_t()> m_pfGetMetricNumberOfValidPoints;
    nsstd::function<std::vector<double>()> m_pfGetOptimizerScales;
    nsstd::function<std::string()> m_pfGetOptimizerStopConditionDescription;


    nsstd::function<unsigned int()> m_pfGetCurrentLevel;

    nsstd::function<void (itk::TransformBase *outTransform)> m_pfUpdateWithBestValue;

    template < class TMemberFunctionPointer >
      struct EvaluateMemberFunctionAddressor
    {
      typedef typename ::detail::FunctionTraits<TMemberFunctionPointer>::ClassType ObjectType;

      template< typename TImageType >
      TMemberFunctionPointer operator() ( void ) 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>;
    nsstd::auto_ptr<detail::MemberFunctionFactory<MemberFunctionType> > m_MemberFactory;
    nsstd::auto_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