elxAdaptiveStochasticLBFGS.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright UMC Utrecht and contributors
 *
 *  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 __elxAdaptiveStochasticLBFGS_h
#define __elxAdaptiveStochasticLBFGS_h
 
#include "elxIncludes.h" // include first to avoid MSVS warning
#include "itkAdaptiveStochasticLBFGSOptimizer.h"
 
#include "elxProgressCommand.h"
#include "itkAdvancedTransform.h"
#include "itkMersenneTwisterRandomVariateGenerator.h"
#include "itkComputeJacobianTerms.h"
#include "itkComputeDisplacementDistribution.h"
#include "itkPlatformMultiThreader.h"
#include "itkImageRandomSampler.h"
#include "itkLineSearchOptimizer.h"
#include "itkMoreThuenteLineSearchOptimizer.h"
 
 
namespace elastix
{
template <class TElastix>
class AdaptiveStochasticLBFGS :
  public itk::AdaptiveStochasticLBFGSOptimizer,
  public OptimizerBase<TElastix>
{
public:
 
  typedef AdaptiveStochasticLBFGS           Self;
  typedef AdaptiveStochasticLBFGSOptimizer  Superclass1;
  typedef OptimizerBase<TElastix>           Superclass2;
  typedef itk::SmartPointer<Self>           Pointer;
  typedef itk::SmartPointer<const Self>     ConstPointer;
 
  itkNewMacro( Self );
 
  itkTypeMacro( AdaptiveStochasticLBFGS,
    AdaptiveStochasticLBFGSOptimizer );
 
  elxClassNameMacro( "AdaptiveStochasticLBFGS" );
 
  typedef Superclass1::CostFunctionType               CostFunctionType;
  typedef Superclass1::CostFunctionPointer            CostFunctionPointer;
  typedef Superclass1::StopConditionType              StopConditionType;
 
  typedef typename Superclass2::ElastixType           ElastixType;
  typedef typename Superclass2::ElastixPointer        ElastixPointer;
  typedef typename Superclass2::ConfigurationType     ConfigurationType;
  typedef typename Superclass2::ConfigurationPointer  ConfigurationPointer;
  typedef typename Superclass2::RegistrationType      RegistrationType;
  typedef typename Superclass2::RegistrationPointer   RegistrationPointer;
  typedef typename Superclass2::ITKBaseType           ITKBaseType;
  typedef itk::SizeValueType                          SizeValueType;
 
//  typedef LineSearchOptimizer           LineSearchOptimizerType;
 
//  typedef LineSearchOptimizerType::Pointer LineSearchOptimizerPointer;
 
  typedef typename Superclass1::ParametersType        ParametersType;
 
  typedef itk::LineSearchOptimizer                    LineSearchOptimizerType;
 
  typedef LineSearchOptimizerType::Pointer            LineSearchOptimizerPointer;
  typedef itk::MoreThuenteLineSearchOptimizer         LineOptimizerType;
  typedef LineOptimizerType::Pointer                  LineOptimizerPointer;
 
  void BeforeRegistration( void ) override;
  void BeforeEachResolution( void ) override;
  void AfterEachResolution( void ) override;
  void AfterEachIteration( void ) override;
  void AfterRegistration( void ) override;
 
  void StartOptimization( void ) override;
 
  virtual void LBFGSUpdate( void );
 
  void AdvanceOneStep( void ) override;
 
  void ResumeOptimization( void ) override;
 
  void MetricErrorResponse( itk::ExceptionObject & err ) override;
 
  void StopOptimization( void ) override;
 
  itkSetMacro( AutomaticParameterEstimation, bool );
  itkGetConstMacro( AutomaticParameterEstimation, bool );
 
  itkSetMacro( AutomaticLBFGSStepsizeEstimation, bool );
  itkGetConstMacro( AutomaticLBFGSStepsizeEstimation, bool );
 
  itkSetMacro( MaximumStepLength, double );
  itkGetConstMacro( MaximumStepLength, double );
 
  itkSetMacro( MaximumNumberOfSamplingAttempts, SizeValueType );
 
  itkGetConstReferenceMacro( MaximumNumberOfSamplingAttempts, SizeValueType );
 
  typedef unsigned int  ThreadIdType;
 
  void SetNumberOfWorkUnits( ThreadIdType numberOfThreads )
  {
    this->m_Threader->SetNumberOfWorkUnits( numberOfThreads );
  }
 
protected:
 
  typedef typename RegistrationType::FixedImageType   FixedImageType;
  typedef typename RegistrationType::MovingImageType  MovingImageType;
 
  typedef typename FixedImageType::RegionType         FixedImageRegionType;
  typedef typename FixedImageType::IndexType          FixedImageIndexType;
  typedef typename FixedImageType::PointType          FixedImagePointType;
  typedef typename RegistrationType::ITKBaseType      itkRegistrationType;
  typedef typename itkRegistrationType::TransformType TransformType;
  typedef typename TransformType::JacobianType        JacobianType;
  typedef itk::ComputeJacobianTerms<
    FixedImageType,TransformType >                    ComputeJacobianTermsType;
  typedef typename JacobianType::ValueType            JacobianValueType;
  struct SettingsType { double a, A, alpha, fmax, fmin, omega; };
  typedef typename std::vector<SettingsType>          SettingsVectorType;
 
  typedef itk::ComputeDisplacementDistribution<
    FixedImageType,TransformType >                    ComputeDisplacementDistributionType;
 
  typedef itk::ImageSamplerBase<FixedImageType>       ImageSamplerBaseType;
  typedef typename ImageSamplerBaseType::Pointer      ImageSamplerBasePointer;
  typedef itk::ImageRandomSamplerBase<FixedImageType> ImageRandomSamplerBaseType;
  typedef typename
    ImageRandomSamplerBaseType::Pointer               ImageRandomSamplerBasePointer;
  typedef
    itk::ImageRandomCoordinateSampler<FixedImageType> ImageRandomCoordinateSamplerType;
  typedef typename
    ImageRandomCoordinateSamplerType::Pointer         ImageRandomCoordinateSamplerPointer;
  typedef typename ImageSamplerBaseType::ImageSampleType       ImageSampleType;
 
  typedef itk::ImageRandomSampler< FixedImageType >   ImageRandomSamplerType;
  typedef typename ImageRandomSamplerType::Pointer    ImageRandomSamplerPointer;
  typedef typename
    ImageRandomSamplerType::ImageSampleContainerType  ImageRadomSampleContainerType;
  typedef typename
    ImageRadomSampleContainerType::Pointer            ImageRadomSampleContainerPointer;
 
  typedef itk::ImageGridSampler< FixedImageType >     ImageGridSamplerType;
  typedef typename ImageGridSamplerType::Pointer      ImageGridSamplerPointer;
  typedef typename
    ImageGridSamplerType::ImageSampleContainerType    ImageSampleContainerType;
  typedef typename ImageSampleContainerType::Pointer  ImageSampleContainerPointer;
 
  typedef itk::Statistics::MersenneTwisterRandomVariateGenerator RandomGeneratorType;
  typedef ProgressCommand                             ProgressCommandType;
  typedef typename ProgressCommand::Pointer           ProgressCommandPointer;
 
  typedef JacobianType                                TransformJacobianType;
  itkStaticConstMacro( FixedImageDimension, unsigned int, FixedImageType::ImageDimension );
  itkStaticConstMacro( MovingImageDimension, unsigned int, MovingImageType::ImageDimension );
  typedef typename TransformType::ScalarType          CoordinateRepresentationType;
  typedef itk::AdvancedTransform<
    CoordinateRepresentationType,
    itkGetStaticConstMacro(FixedImageDimension),
    itkGetStaticConstMacro(MovingImageDimension) >    AdvancedTransformType;
  typedef typename
    AdvancedTransformType::NonZeroJacobianIndicesType NonZeroJacobianIndicesType;
 
  typedef itk::Array< double >               RhoType;
  typedef std::vector< ParametersType >      SType;
  typedef std::vector< DerivativeType >      YType;
  typedef itk::Array< double >               DiagonalMatrixType;
 
  AdaptiveStochasticLBFGS();
  ~AdaptiveStochasticLBFGS() override {};
 
  SettingsVectorType m_SettingsVector;
 
  SizeValueType m_NumberOfGradientMeasurements;
  SizeValueType m_NumberOfJacobianMeasurements;
  SizeValueType m_NumberOfSamplesForExactGradient;
 
  typename AdvancedTransformType::Pointer           m_AdvancedTransform;
 
  typename RandomGeneratorType::Pointer             m_RandomGenerator;
 
  double m_SigmoidScaleFactor;
 
  virtual void CheckForAdvancedTransform( void );
 
  virtual void PrintSettingsVector( const SettingsVectorType & settings ) const;
 
  virtual void AutomaticParameterEstimation( void );
 
  virtual void AutomaticParameterEstimationOriginal( void );
 
  virtual void AutomaticParameterEstimationUsingDisplacementDistribution( void );
 
  virtual void AutomaticLBFGSStepsizeEstimation( void );
 
  virtual void SampleGradients( const ParametersType & mu0,
    double perturbationSigma, double & gg, double & ee );
 
  virtual void GetScaledDerivativeWithExceptionHandling(
    const ParametersType & parameters, DerivativeType & derivative );
 
  virtual void AddRandomPerturbation( ParametersType & parameters, double sigma );
 
  virtual void StoreCurrentPoint(
    const ParametersType & step,
    const DerivativeType & grad_dif );
  virtual void ComputeDiagonalMatrix( DiagonalMatrixType & diag_H0 );
 
  virtual void ComputeSearchDirection(
    const DerivativeType & gradient,
    DerivativeType & searchDir );
 
  DerivativeType                m_PreviousCurvatureGradient;
 
  double                        m_NoiseFactor;
  unsigned int                  m_CurrentT;
  unsigned int                  m_PreviousT;
  unsigned int                  m_Bound;
 
  RhoType m_Rho;
  SType   m_S;
  YType   m_Y;
  RhoType m_HessianFillValue;
  double  m_WindowScale;
 
private:
 
  AdaptiveStochasticLBFGS( const Self& );  // purposely not implemented
  void operator=( const Self& );           // purposely not implemented
 
  // multi-threaded AdvanceOneStep:
  struct MultiThreaderParameterType
  {
    ParametersType *  t_NewPosition;
    Self *            t_Optimizer;
  };
 
  static itk::ITK_THREAD_RETURN_TYPE AdvanceOneStepThreaderCallback( void * arg );
 
  inline void ThreadedAdvanceOneStep( ThreadIdType threadId, ParametersType & newPosition );
 
  bool    m_AutomaticParameterEstimation;
  bool    m_AutomaticLBFGSStepsizeEstimation;
  double  m_MaximumStepLength;
 
  SizeValueType m_MaximumNumberOfSamplingAttempts;
  SizeValueType m_CurrentNumberOfSamplingAttempts;
  SizeValueType m_PreviousErrorAtIteration;
  bool          m_AutomaticParameterEstimationDone;
 
  SizeValueType m_OutsideIterations;
 
  SizeValueType m_MaxBandCovSize;
  SizeValueType m_NumberOfBandStructureSamples;
  SizeValueType m_NumberOfInnerLoopSamples;
  SizeValueType m_NumberOfSpatialSamples;
 
  bool m_UseNoiseCompensation;
  bool m_OriginalButSigmoidToDefault;
  bool m_UseNoiseFactor;
 
  bool                        m_UseAdaptiveLBFGSStepSizes;
 
}; // end class AdaptiveStochasticLBFGS
 
 
} // end namespace elastix
 
#ifndef ITK_MANUAL_INSTANTIATION
#include "elxAdaptiveStochasticLBFGS.hxx"
#endif
 
#endif // end #ifndef __elxAdaptiveStochasticLBFGS_h