itkAdvancedImageToImageMetric.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 itkAdvancedImageToImageMetric_h
#define itkAdvancedImageToImageMetric_h
 
#include "itkImageToImageMetric.h"
 
#include "itkImageSamplerBase.h"
#include "itkGradientImageFilter.h"
#include "itkBSplineInterpolateImageFunction.h"
#include "itkReducedDimensionBSplineInterpolateImageFunction.h"
#include "itkAdvancedLinearInterpolateImageFunction.h"
#include "itkLimiterFunctionBase.h"
#include "itkFixedArray.h"
#include "itkAdvancedTransform.h"
 
#include "itkImageMaskSpatialObject.h"
 
// Needed for checking for B-spline for faster implementation
#include "itkAdvancedBSplineDeformableTransform.h"
#include "itkAdvancedCombinationTransform.h"
 
#include "elxDefaultConstruct.h"
 
#include <cassert>
#include <memory> // For unique_ptr.
#include <typeinfo>
 
namespace itk
{

 
template <typename TFixedImage, typename TMovingImage>
class ITK_TEMPLATE_EXPORT AdvancedImageToImageMetric : public ImageToImageMetric<TFixedImage, TMovingImage>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(AdvancedImageToImageMetric);

  using Self = AdvancedImageToImageMetric;
  using Superclass = ImageToImageMetric<TFixedImage, TMovingImage>;
  using Pointer = SmartPointer<Self>;
  using ConstPointer = SmartPointer<const Self>;

  itkOverrideGetNameOfClassMacro(AdvancedImageToImageMetric);

  itkStaticConstMacro(MovingImageDimension, unsigned int, TMovingImage::ImageDimension);
  itkStaticConstMacro(FixedImageDimension, unsigned int, TFixedImage::ImageDimension);

  using typename Superclass::CoordinateRepresentationType;
  using typename Superclass::MovingImageType;
  using typename Superclass::MovingImagePixelType;
  using MovingImagePointer = typename MovingImageType::Pointer;
  using typename Superclass::MovingImageConstPointer;
  using typename Superclass::FixedImageType;
  using FixedImagePointer = typename FixedImageType::Pointer;
  using typename Superclass::FixedImageConstPointer;
  using typename Superclass::FixedImageRegionType;
  using typename Superclass::TransformType;
  using typename Superclass::TransformPointer;
  using typename Superclass::InputPointType;
  using typename Superclass::OutputPointType;
  using typename Superclass::TransformParametersType;
  using typename Superclass::TransformJacobianType;
  using typename Superclass::InterpolatorType;
  using typename Superclass::InterpolatorPointer;
  using typename Superclass::RealType;
  using typename Superclass::GradientPixelType;
  using typename Superclass::GradientImageType;
  using typename Superclass::GradientImagePointer;
 
  // Overrule the mask type from its base class, ITK ImageToImageMetric.
  using FixedImageMaskType = ImageMaskSpatialObject<Self::FixedImageDimension>;
  using FixedImageMaskPointer = SmartPointer<FixedImageMaskType>;
  using FixedImageMaskConstPointer = SmartPointer<const FixedImageMaskType>;
  using MovingImageMaskType = ImageMaskSpatialObject<Self::MovingImageDimension>;
  using MovingImageMaskPointer = SmartPointer<MovingImageMaskType>;
  using MovingImageMaskConstPointer = SmartPointer<const MovingImageMaskType>;
 
  using typename Superclass::MeasureType;
  using typename Superclass::DerivativeType;
  using DerivativeValueType = typename DerivativeType::ValueType;
  using typename Superclass::ParametersType;

  using FixedImagePixelType = typename FixedImageType::PixelType;
  using MovingImageRegionType = typename MovingImageType::RegionType;
  using MovingImageDerivativeScalesType = FixedArray<double, Self::MovingImageDimension>;

  using ImageSamplerType = ImageSamplerBase<FixedImageType>;
  using ImageSamplerPointer = typename ImageSamplerType::Pointer;
  using ImageSampleContainerType = typename ImageSamplerType::OutputVectorContainerType;
  using ImageSampleContainerPointer = typename ImageSamplerType::OutputVectorContainerPointer;

  using FixedImageLimiterType = LimiterFunctionBase<RealType, FixedImageDimension>;
  using FixedImageLimiterPointer = typename FixedImageLimiterType::Pointer;
  using FixedImageLimiterOutputType = typename FixedImageLimiterType::OutputType;
  using MovingImageLimiterType = LimiterFunctionBase<RealType, MovingImageDimension>;
  using MovingImageLimiterPointer = typename MovingImageLimiterType::Pointer;
  using MovingImageLimiterOutputType = typename MovingImageLimiterType::OutputType;

  using ScalarType = typename TransformType::ScalarType;
  using AdvancedTransformType = AdvancedTransform<ScalarType, FixedImageDimension, MovingImageDimension>;
  using NumberOfParametersType = typename AdvancedTransformType::NumberOfParametersType;

  using CombinationTransformType = AdvancedCombinationTransform<ScalarType, FixedImageDimension>;
  using BSplineOrder1TransformType = AdvancedBSplineDeformableTransform<ScalarType, FixedImageDimension, 1>;
  using BSplineOrder2TransformType = AdvancedBSplineDeformableTransform<ScalarType, FixedImageDimension, 2>;
  using BSplineOrder3TransformType = AdvancedBSplineDeformableTransform<ScalarType, FixedImageDimension, 3>;
  using BSplineOrder1TransformPointer = typename BSplineOrder1TransformType::Pointer;
  using BSplineOrder2TransformPointer = typename BSplineOrder2TransformType::Pointer;
  using BSplineOrder3TransformPointer = typename BSplineOrder3TransformType::Pointer;

  using ThreadInfoType = MultiThreaderBase::WorkUnitInfo;

 
  virtual void
  SetFixedImageMask(const FixedImageMaskType * const arg)
  {
    assert(arg == nullptr || typeid(*arg) == typeid(FixedImageMaskType));
    Superclass::SetFixedImageMask(arg);
  }
 
  virtual void
  SetMovingImageMask(const MovingImageMaskType * const arg)
  {
    assert(arg == nullptr || typeid(*arg) == typeid(MovingImageMaskType));
    Superclass::SetMovingImageMask(arg);
  }
 
  const FixedImageMaskType *
  GetFixedImageMask() const override
  {
    const auto * const mask = Superclass::GetFixedImageMask();
    assert(mask == nullptr || typeid(*mask) == typeid(FixedImageMaskType));
    return static_cast<const FixedImageMaskType *>(mask);
  }
 
  const MovingImageMaskType *
  GetMovingImageMask() const override
  {
    const auto * const mask = Superclass::GetMovingImageMask();
    assert(mask == nullptr || typeid(*mask) == typeid(MovingImageMaskType));
    return static_cast<const MovingImageMaskType *>(mask);
  }

  virtual void
  SetTransform(AdvancedTransformType * arg)
  {
    this->Superclass::SetTransform(arg);
    if (m_AdvancedTransform != arg)
    {
      m_AdvancedTransform = arg;
      this->Modified();
    }
  }
 

  AdvancedTransformType *
  GetTransform() override
  {
    return m_AdvancedTransform.GetPointer();
  }

  const AdvancedTransformType *
  GetTransform() const override
  {
    return m_AdvancedTransform.GetPointer();
  }
 

  itkSetObjectMacro(ImageSampler, ImageSamplerType);
  ImageSamplerType *
  GetImageSampler() const
  {
    return m_ImageSampler.GetPointer();
  }
 

  itkGetConstMacro(UseImageSampler, bool);

  itkSetMacro(RequiredRatioOfValidSamples, double);
  itkGetConstMacro(RequiredRatioOfValidSamples, double);

  itkSetObjectMacro(MovingImageLimiter, MovingImageLimiterType);
  itkGetConstObjectMacro(MovingImageLimiter, MovingImageLimiterType);
  itkSetObjectMacro(FixedImageLimiter, FixedImageLimiterType);
  itkGetConstObjectMacro(FixedImageLimiter, FixedImageLimiterType);

  itkSetMacro(MovingLimitRangeRatio, double);
  itkGetConstMacro(MovingLimitRangeRatio, double);
  itkSetMacro(FixedLimitRangeRatio, double);
  itkGetConstMacro(FixedLimitRangeRatio, double);

  itkGetConstMacro(UseFixedImageLimiter, bool);
  itkGetConstMacro(UseMovingImageLimiter, bool);

  itkSetMacro(UseMovingImageDerivativeScales, bool);
  itkGetConstMacro(UseMovingImageDerivativeScales, bool);
 
  itkSetMacro(ScaleGradientWithRespectToMovingImageOrientation, bool);
  itkGetConstMacro(ScaleGradientWithRespectToMovingImageOrientation, bool);
 
  itkSetMacro(MovingImageDerivativeScales, MovingImageDerivativeScalesType);
  itkGetConstReferenceMacro(MovingImageDerivativeScales, MovingImageDerivativeScalesType);

  void
  Initialize() override;

  itkSetMacro(UseMetricSingleThreaded, bool);
  itkGetConstReferenceMacro(UseMetricSingleThreaded, bool);
  itkBooleanMacro(UseMetricSingleThreaded);

  // \todo: maybe these can be united, check base class.
  itkSetMacro(UseMultiThread, bool);
  itkGetConstReferenceMacro(UseMultiThread, bool);
  itkBooleanMacro(UseMultiThread);

  virtual void
  BeforeThreadedGetValueAndDerivative(const TransformParametersType & parameters) const;
 
  void
  SetRandomVariateGenerator(Statistics::MersenneTwisterRandomVariateGenerator & randomVariateGenerator)
  {
    m_RandomVariateGenerator = &randomVariateGenerator;
  }
 
protected:
  AdvancedImageToImageMetric();

  ~AdvancedImageToImageMetric() override = default;

  void
  PrintSelf(std::ostream & os, Indent indent) const override;


  using FixedImageIndexType = typename FixedImageType::IndexType;
  using FixedImageIndexValueType = typename FixedImageIndexType::IndexValueType;
  using MovingImageIndexType = typename MovingImageType::IndexType;
  using FixedImagePointType = typename TransformType::InputPointType;
  using MovingImagePointType = typename TransformType::OutputPointType;
  using MovingImageContinuousIndexType = typename InterpolatorType::ContinuousIndexType;

  using BSplineInterpolatorType =
    BSplineInterpolateImageFunction<MovingImageType, CoordinateRepresentationType, double>;
  using BSplineInterpolatorPointer = typename BSplineInterpolatorType::Pointer;
  using BSplineInterpolatorFloatType =
    BSplineInterpolateImageFunction<MovingImageType, CoordinateRepresentationType, float>;
  using BSplineInterpolatorFloatPointer = typename BSplineInterpolatorFloatType::Pointer;
  using ReducedBSplineInterpolatorType =
    ReducedDimensionBSplineInterpolateImageFunction<MovingImageType, CoordinateRepresentationType, double>;
  using ReducedBSplineInterpolatorPointer = typename ReducedBSplineInterpolatorType::Pointer;
  using LinearInterpolatorType = AdvancedLinearInterpolateImageFunction<MovingImageType, CoordinateRepresentationType>;
  using LinearInterpolatorPointer = typename LinearInterpolatorType::Pointer;
  using MovingImageDerivativeType = typename BSplineInterpolatorType::CovariantVectorType;

  using NonZeroJacobianIndicesType = typename AdvancedTransformType::NonZeroJacobianIndicesType;


  mutable ImageSamplerPointer m_ImageSampler{ nullptr };

  typename AdvancedTransformType::Pointer m_AdvancedTransform{ nullptr };

  mutable bool m_TransformIsBSpline{ false };

  FixedImagePixelType          m_FixedImageTrueMin{ 0 };
  FixedImagePixelType          m_FixedImageTrueMax{ 1 };
  MovingImagePixelType         m_MovingImageTrueMin{ 0 };
  MovingImagePixelType         m_MovingImageTrueMax{ 1 };
  FixedImageLimiterOutputType  m_FixedImageMinLimit{ 0 };
  FixedImageLimiterOutputType  m_FixedImageMaxLimit{ 1 };
  MovingImageLimiterOutputType m_MovingImageMinLimit{ 0 };
  MovingImageLimiterOutputType m_MovingImageMaxLimit{ 1 };


  virtual void
  ThreadedGetValue(ThreadIdType) const
  {}

  virtual void
  AfterThreadedGetValue(MeasureType &) const
  {}

  static ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION
  GetValueThreaderCallback(void * arg);

  void
  LaunchGetValueThreaderCallback() const;

  virtual void
  ThreadedGetValueAndDerivative(ThreadIdType) const
  {}

  virtual void
  AfterThreadedGetValueAndDerivative(MeasureType &, DerivativeType &) const
  {}

  static ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION
  GetValueAndDerivativeThreaderCallback(void * arg);

  void
  LaunchGetValueAndDerivativeThreaderCallback() const;

  static ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION
  AccumulateDerivativesThreaderCallback(void * arg);

  bool m_UseMetricSingleThreaded{ true };
  bool m_UseMultiThread{ false };

  struct MultiThreaderParameterType
  {
    // To give the threads access to all members.
    AdvancedImageToImageMetric * st_Metric;
    // Used for accumulating derivatives
    DerivativeValueType * st_DerivativePointer;
    DerivativeValueType   st_NormalizationFactor;
  };
  mutable MultiThreaderParameterType m_ThreaderMetricParameters{};

 
  // test per thread struct with padding and alignment
  struct GetValueAndDerivativePerThreadStruct
  {
    SizeValueType  st_NumberOfPixelsCounted;
    MeasureType    st_Value;
    DerivativeType st_Derivative;
  };
  itkPadStruct(ITK_CACHE_LINE_ALIGNMENT,
               GetValueAndDerivativePerThreadStruct,
               PaddedGetValueAndDerivativePerThreadStruct);
  itkAlignedTypedef(ITK_CACHE_LINE_ALIGNMENT,
                    PaddedGetValueAndDerivativePerThreadStruct,
                    AlignedGetValueAndDerivativePerThreadStruct);
  mutable std::unique_ptr<AlignedGetValueAndDerivativePerThreadStruct[]> m_GetValueAndDerivativePerThreadVariables{
    nullptr
  };
  mutable ThreadIdType m_GetValueAndDerivativePerThreadVariablesSize{ 0 };

  virtual void
  InitializeThreadingParameters() const;



  virtual void
  InitializeImageSampler();

  itkSetMacro(UseImageSampler, bool);

  void
  CheckNumberOfSamples() const;


  void
  CheckForBSplineInterpolator();

  virtual bool
  EvaluateMovingImageValueAndDerivative(const MovingImagePointType & mappedPoint,
                                        RealType &                   movingImageValue,
                                        MovingImageDerivativeType *  gradient) const
  {
    return EvaluateMovingImageValueAndDerivativeWithOptionalThreadId(mappedPoint, movingImageValue, gradient);
  }
 
  /* A faster version of `EvaluateMovingImageValueAndDerivative`: Non-virtual, using multithreading, and doing less
   * dynamic memory allocation/decallocation operations, internally. */
  bool
  FastEvaluateMovingImageValueAndDerivative(const MovingImagePointType & mappedPoint,
                                            RealType &                   movingImageValue,
                                            MovingImageDerivativeType *  gradient,
                                            const ThreadIdType           threadId) const
  {
    return EvaluateMovingImageValueAndDerivativeWithOptionalThreadId(mappedPoint, movingImageValue, gradient, threadId);
  }

  virtual void
  EvaluateTransformJacobianInnerProduct(const TransformJacobianType &     jacobian,
                                        const MovingImageDerivativeType & movingImageDerivative,
                                        DerivativeType &                  imageJacobian) const;


  void
  CheckForAdvancedTransform();

  void
  CheckForBSplineTransform() const;

  MovingImagePointType
  TransformPoint(const FixedImagePointType & fixedImagePoint) const;

  bool
  EvaluateTransformJacobian(const FixedImagePointType &  fixedImagePoint,
                            TransformJacobianType &      jacobian,
                            NonZeroJacobianIndicesType & nzji) const;

  virtual bool
  IsInsideMovingMask(const MovingImagePointType & point) const;

  void
  InitializeLimiters();

  itkSetMacro(UseFixedImageLimiter, bool);
  itkSetMacro(UseMovingImageLimiter, bool);
 
  double m_FixedLimitRangeRatio{ 0.01 };
  double m_MovingLimitRangeRatio{ 0.01 };
 
  // Prevent accidentally calling SetFixedImageMask or SetMovingImageMask through the ITK ImageToImageMetric interface.
  void
  SetFixedImageMask(typename Superclass::FixedImageMaskType *) final
  {
    itkExceptionMacro("Intentionally left unimplemented!");
  }
 
  void
  SetFixedImageMask(const typename Superclass::FixedImageMaskType *) final
  {
    itkExceptionMacro("Intentionally left unimplemented!");
  }
 
  void
  SetMovingImageMask(typename Superclass::MovingImageMaskType *) final
  {
    itkExceptionMacro("Intentionally left unimplemented!");
  }
 
  void
  SetMovingImageMask(const typename Superclass::MovingImageMaskType *) final
  {
    itkExceptionMacro("Intentionally left unimplemented!");
  }
 
  Statistics::MersenneTwisterRandomVariateGenerator &
  GetRandomVariateGenerator()
  {
    return *m_RandomVariateGenerator;
  }
 
  // Note: Bypasses logical const-correctness
  Statistics::MersenneTwisterRandomVariateGenerator &
  GetMutableRandomVariateGenerator() const
  {
    return *m_RandomVariateGenerator;
  }
 
  // Protected using-declaration, to avoid `-Woverloaded-virtual` warnings from GCC (GCC 11.4) or clang (macos-12).
  using Superclass::SetTransform;
 
private:
  template <typename... TOptionalThreadId>
  bool
  EvaluateMovingImageValueAndDerivativeWithOptionalThreadId(const MovingImagePointType & mappedPoint,
                                                            RealType &                   movingImageValue,
                                                            MovingImageDerivativeType *  gradient,
                                                            const TOptionalThreadId... optionalThreadId) const;

  FixedImageLimiterPointer          m_FixedImageLimiter{ nullptr };
  MovingImageLimiterPointer         m_MovingImageLimiter{ nullptr };
  LinearInterpolatorPointer         m_LinearInterpolator{ nullptr };
  BSplineInterpolatorPointer        m_BSplineInterpolator{ nullptr };
  BSplineInterpolatorFloatPointer   m_BSplineInterpolatorFloat{ nullptr };
  ReducedBSplineInterpolatorPointer m_ReducedBSplineInterpolator{ nullptr };

  bool   m_UseImageSampler{ false };
  bool   m_UseFixedImageLimiter{ false };
  bool   m_UseMovingImageLimiter{ false };
  double m_RequiredRatioOfValidSamples{ 0.25 };
  bool   m_UseMovingImageDerivativeScales{ false };
  bool   m_ScaleGradientWithRespectToMovingImageOrientation{ false };
 
  MovingImageDerivativeScalesType m_MovingImageDerivativeScales{ MovingImageDerivativeScalesType::Filled(1.0) };
 
  mutable elx::DefaultConstruct<Statistics::MersenneTwisterRandomVariateGenerator> m_DefaultRandomVariateGenerator{};
  Statistics::MersenneTwisterRandomVariateGenerator * m_RandomVariateGenerator{ &m_DefaultRandomVariateGenerator };
 
  // Private using-declarations, to avoid `-Woverloaded-virtual` warnings from GCC (GCC 11.4) or clang (macos-12).
  using Superclass::TransformPoint;
 
  struct DummyMask
  {};
 
  // Prevent accidentally accessing m_FixedImageMask or m_MovingImageMask directly from the ITK's ImageToImageMetric.
  static constexpr DummyMask m_FixedImageMask{};
  static constexpr DummyMask m_MovingImageMask{};
};
 
} // end namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkAdvancedImageToImageMetric.hxx"
#endif
 
#endif // end #ifndef itkAdvancedImageToImageMetric_h