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 <vnl/vnl_sparse_matrix.h>
 
#include "itkImageMaskSpatialObject.h"
 
// Needed for checking for B-spline for faster implementation
#include "itkAdvancedBSplineDeformableTransform.h"
#include "itkAdvancedCombinationTransform.h"
 
#include "itkPlatformMultiThreader.h"
 
#include <memory> // For unique_ptr.
 
namespace itk
{
 
template <class TFixedImage, class 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>;
 
  itkTypeMacro(AdvancedImageToImageMetric, ImageToImageMetric);
 
  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;
  using typename Superclass::GradientImageFilterType;
  using typename Superclass::GradientImageFilterPointer;
  using typename Superclass::FixedImageMaskType;
  using typename Superclass::FixedImageMaskPointer;
  using typename Superclass::MovingImageMaskType;
  using typename Superclass::MovingImageMaskPointer;
  using typename Superclass::MeasureType;
  using typename Superclass::DerivativeType;
  using DerivativeValueType = typename DerivativeType::ValueType;
  using typename Superclass::ParametersType;
 
  using FixedImageMaskSpatialObject2Type = ImageMaskSpatialObject<Self::FixedImageDimension>;
  using MovingImageMaskSpatialObject2Type = ImageMaskSpatialObject<Self::MovingImageDimension>;
 
  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 HessianValueType = typename DerivativeType::ValueType;
  using HessianType = vnl_sparse_matrix<HessianValueType>;
 
  using ThreaderType = itk::PlatformMultiThreader;
  using ThreadInfoType = typename ThreaderType::WorkUnitInfo;
 
  virtual void
  SetTransform(AdvancedTransformType * arg)
  {
    this->Superclass::SetTransform(arg);
    if (this->m_AdvancedTransform != arg)
    {
      this->m_AdvancedTransform = arg;
      this->Modified();
    }
  }
 
 
  const AdvancedTransformType *
  GetTransform() const override
  {
    return this->m_AdvancedTransform.GetPointer();
  }
 
 
  itkSetObjectMacro(ImageSampler, ImageSamplerType);
  virtual ImageSamplerType *
  GetImageSampler() const
  {
    return this->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;
 
  virtual void
  GetSelfHessian(const TransformParametersType & parameters, HessianType & H) const;
 
  virtual void
  SetNumberOfWorkUnits(ThreadIdType numberOfThreads);
 
  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;
 
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 CentralDifferenceGradientFilterType = GradientImageFilter<MovingImageType, RealType, RealType>;
  using CentralDifferenceGradientFilterPointer = typename CentralDifferenceGradientFilterType::Pointer;
 
  using NonZeroJacobianIndicesType = typename AdvancedTransformType::NonZeroJacobianIndicesType;
 
  mutable ImageSamplerPointer m_ImageSampler{ nullptr };
 
  bool                              m_InterpolatorIsLinear{ false };
  bool                              m_InterpolatorIsBSpline{ false };
  bool                              m_InterpolatorIsBSplineFloat{ false };
  bool                              m_InterpolatorIsReducedBSpline{ false };
  LinearInterpolatorPointer         m_LinearInterpolator{ nullptr };
  BSplineInterpolatorPointer        m_BSplineInterpolator{ nullptr };
  BSplineInterpolatorFloatPointer   m_BSplineInterpolatorFloat{ nullptr };
  ReducedBSplineInterpolatorPointer m_ReducedBSplineInterpolator{ nullptr };
 
  CentralDifferenceGradientFilterPointer m_CentralDifferenceGradientFilter{ nullptr };
 
  bool                                    m_TransformIsAdvanced{ false };
  typename AdvancedTransformType::Pointer m_AdvancedTransform{ nullptr };
  mutable bool                            m_TransformIsBSpline{ false };
 
  FixedImageLimiterPointer     m_FixedImageLimiter{ nullptr };
  MovingImageLimiterPointer    m_MovingImageLimiter{ nullptr };
  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 threadID)
  {}
 
  virtual void
  AfterThreadedGetValue(MeasureType & value) const
  {}
 
  static ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION
  GetValueThreaderCallback(void * arg);
 
  void
  LaunchGetValueThreaderCallback() const;
 
  virtual void
  ThreadedGetValueAndDerivative(ThreadIdType threadID)
  {}
 
  virtual void
  AfterThreadedGetValueAndDerivative(MeasureType & value, DerivativeType & derivative) 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 };
  bool m_UseOpenMP;
 
  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 GetValuePerThreadStruct
  {
    SizeValueType st_NumberOfPixelsCounted;
    MeasureType   st_Value;
  };
  itkPadStruct(ITK_CACHE_LINE_ALIGNMENT, GetValuePerThreadStruct, PaddedGetValuePerThreadStruct);
  itkAlignedTypedef(ITK_CACHE_LINE_ALIGNMENT, PaddedGetValuePerThreadStruct, AlignedGetValuePerThreadStruct);
  mutable std::unique_ptr<AlignedGetValuePerThreadStruct[]> m_GetValuePerThreadVariables{ nullptr };
  mutable ThreadIdType                                      m_GetValuePerThreadVariablesSize{ 0 };
 
  // 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);
 
  virtual void
  CheckNumberOfSamples(unsigned long wanted, unsigned long found) const;
 
  virtual 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;
 
  virtual void
  CheckForAdvancedTransform();
 
  virtual void
  CheckForBSplineTransform() const;
 
  MovingImagePointType
  TransformPoint(const FixedImagePointType & fixedImagePoint) const;
 
  virtual bool
  EvaluateTransformJacobian(const FixedImagePointType &  fixedImagePoint,
                            TransformJacobianType &      jacobian,
                            NonZeroJacobianIndicesType & nzji) const;
 
  virtual bool
  IsInsideMovingMask(const MovingImagePointType & point) const;
 
  virtual void
  InitializeLimiters();
 
  itkSetMacro(UseFixedImageLimiter, bool);
  itkSetMacro(UseMovingImageLimiter, bool);
 
  double m_FixedLimitRangeRatio{ 0.01 };
  double m_MovingLimitRangeRatio{ 0.01 };
 
private:
  template <typename... TOptionalThreadId>
  bool
  EvaluateMovingImageValueAndDerivativeWithOptionalThreadId(const MovingImagePointType & mappedPoint,
                                                            RealType &                   movingImageValue,
                                                            MovingImageDerivativeType *  gradient,
                                                            const TOptionalThreadId... optionalThreadId) const;
 
  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) };
};
 
} // end namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkAdvancedImageToImageMetric.hxx"
#endif
 
#endif // end #ifndef itkAdvancedImageToImageMetric_h