itkParzenWindowHistogramImageToImageMetric.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 itkParzenWindowHistogramImageToImageMetric_h
#define itkParzenWindowHistogramImageToImageMetric_h
 
#include "itkAdvancedImageToImageMetric.h"
#include "itkKernelFunctionBase2.h"
#include <vector>
 
 
namespace itk
{
 
template <class TFixedImage, class TMovingImage>
class ITK_TEMPLATE_EXPORT ParzenWindowHistogramImageToImageMetric
  : public AdvancedImageToImageMetric<TFixedImage, TMovingImage>
{
public:
  ITK_DISALLOW_COPY_AND_MOVE(ParzenWindowHistogramImageToImageMetric);

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

  itkTypeMacro(ParzenWindowHistogramImageToImageMetric, AdvancedImageToImageMetric);

  using typename Superclass::CoordinateRepresentationType;
  using typename Superclass::MovingImageType;
  using typename Superclass::MovingImagePixelType;
  using typename Superclass::MovingImageConstPointer;
  using typename Superclass::FixedImageType;
  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::FixedImageMaskType;
  using typename Superclass::FixedImageMaskPointer;
  using typename Superclass::MovingImageMaskType;
  using typename Superclass::MovingImageMaskPointer;
  using typename Superclass::MeasureType;
  using typename Superclass::DerivativeType;
  using typename Superclass::DerivativeValueType;
  using typename Superclass::ParametersType;
  using typename Superclass::FixedImagePixelType;
  using typename Superclass::MovingImageRegionType;
  using typename Superclass::ImageSamplerType;
  using typename Superclass::ImageSamplerPointer;
  using typename Superclass::ImageSampleContainerType;
  using typename Superclass::ImageSampleContainerPointer;
  using typename Superclass::FixedImageLimiterType;
  using typename Superclass::MovingImageLimiterType;
  using typename Superclass::FixedImageLimiterOutputType;
  using typename Superclass::MovingImageLimiterOutputType;
  using typename Superclass::MovingImageDerivativeScalesType;
  using typename Superclass::ThreadInfoType;

  itkStaticConstMacro(FixedImageDimension, unsigned int, FixedImageType::ImageDimension);

  itkStaticConstMacro(MovingImageDimension, unsigned int, MovingImageType::ImageDimension);

  void
  Initialize() override;

  void
  GetDerivative(const ParametersType & parameters, DerivativeType & Derivative) const override;

  void
  GetValueAndDerivative(const ParametersType & parameters,
                        MeasureType &          value,
                        DerivativeType &       derivative) const override;

  itkSetClampMacro(NumberOfFixedHistogramBins, unsigned long, 4, NumericTraits<unsigned long>::max());
  itkGetConstMacro(NumberOfFixedHistogramBins, unsigned long);

  itkSetClampMacro(NumberOfMovingHistogramBins, unsigned long, 4, NumericTraits<unsigned long>::max());
  itkGetConstMacro(NumberOfMovingHistogramBins, unsigned long);

  itkSetClampMacro(FixedKernelBSplineOrder, unsigned int, 0, 3);
  itkGetConstMacro(FixedKernelBSplineOrder, unsigned int);

  itkSetClampMacro(MovingKernelBSplineOrder, unsigned int, 0, 3);
  itkGetConstMacro(MovingKernelBSplineOrder, unsigned int);

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

  itkSetMacro(UseDerivative, bool);
  itkGetConstMacro(UseDerivative, bool);

  itkSetMacro(UseFiniteDifferenceDerivative, bool);
  itkGetConstMacro(UseFiniteDifferenceDerivative, bool);

  itkSetMacro(FiniteDifferencePerturbation, double);
  itkGetConstMacro(FiniteDifferencePerturbation, double);
 
protected:
  ParzenWindowHistogramImageToImageMetric();

  ~ParzenWindowHistogramImageToImageMetric() override = default;

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


  using typename Superclass::FixedImageIndexType;
  using typename Superclass::FixedImageIndexValueType;
  using OffsetValueType = typename FixedImageType::OffsetValueType;
  using typename Superclass::MovingImageIndexType;
  using typename Superclass::FixedImagePointType;
  using typename Superclass::MovingImagePointType;
  using typename Superclass::MovingImageContinuousIndexType;
  using typename Superclass::BSplineInterpolatorType;
  using typename Superclass::MovingImageDerivativeType;
  using typename Superclass::NonZeroJacobianIndicesType;

  using PDFValueType = double;
  using PDFDerivativeValueType = float;
  using MarginalPDFType = Array<PDFValueType>;
  using JointPDFType = Image<PDFValueType, 2>;
  using JointPDFPointer = typename JointPDFType::Pointer;
  using JointPDFDerivativesType = Image<PDFDerivativeValueType, 3>;
  using JointPDFDerivativesPointer = typename JointPDFDerivativesType::Pointer;
  using IncrementalMarginalPDFType = Image<PDFValueType, 2>;
  using IncrementalMarginalPDFPointer = typename IncrementalMarginalPDFType::Pointer;
  using JointPDFIndexType = JointPDFType::IndexType;
  using JointPDFRegionType = JointPDFType::RegionType;
  using JointPDFSizeType = JointPDFType::SizeType;
  using JointPDFDerivativesIndexType = JointPDFDerivativesType::IndexType;
  using JointPDFDerivativesRegionType = JointPDFDerivativesType::RegionType;
  using JointPDFDerivativesSizeType = JointPDFDerivativesType::SizeType;
  using IncrementalMarginalPDFIndexType = IncrementalMarginalPDFType::IndexType;
  using IncrementalMarginalPDFRegionType = IncrementalMarginalPDFType::RegionType;
  using IncrementalMarginalPDFSizeType = IncrementalMarginalPDFType::SizeType;
  using ParzenValueContainerType = Array<PDFValueType>;

  using KernelFunctionType = KernelFunctionBase2<PDFValueType>;
  using KernelFunctionPointer = typename KernelFunctionType::Pointer;


  mutable double         m_Alpha{ 0.0 };
  mutable DerivativeType m_PerturbedAlphaRight{};
  mutable DerivativeType m_PerturbedAlphaLeft{};

  mutable MarginalPDFType       m_FixedImageMarginalPDF{};
  mutable MarginalPDFType       m_MovingImageMarginalPDF{};
  JointPDFPointer               m_JointPDF{ nullptr };
  JointPDFDerivativesPointer    m_JointPDFDerivatives{ nullptr };
  JointPDFDerivativesPointer    m_IncrementalJointPDFRight{};
  JointPDFDerivativesPointer    m_IncrementalJointPDFLeft{};
  IncrementalMarginalPDFPointer m_FixedIncrementalMarginalPDFRight{ nullptr };
  IncrementalMarginalPDFPointer m_MovingIncrementalMarginalPDFRight{ nullptr };
  IncrementalMarginalPDFPointer m_FixedIncrementalMarginalPDFLeft{ nullptr };
  IncrementalMarginalPDFPointer m_MovingIncrementalMarginalPDFLeft{ nullptr };
  mutable JointPDFRegionType    m_JointPDFWindow{}; // no need for mutable anymore?
  double                        m_MovingImageNormalizedMin{ 0.0 };
  double                        m_FixedImageNormalizedMin{ 0.0 };
  double                        m_FixedImageBinSize{ 0.0 };
  double                        m_MovingImageBinSize{ 0.0 };
  double                        m_FixedParzenTermToIndexOffset{ 0.5 };
  double                        m_MovingParzenTermToIndexOffset{ -1.0 };

  KernelFunctionPointer m_FixedKernel{ nullptr };
  KernelFunctionPointer m_MovingKernel{ nullptr };
  KernelFunctionPointer m_DerivativeMovingKernel{ nullptr };

  void
  InitializeThreadingParameters() const override;

  void
  ThreadedComputePDFs(ThreadIdType threadId);

  void
  AfterThreadedComputePDFs() const;

  static ITK_THREAD_RETURN_FUNCTION_CALL_CONVENTION
  ComputePDFsThreaderCallback(void * arg);

  void
  LaunchComputePDFsThreaderCallback() const;

  static void
  EvaluateParzenValues(double                     parzenWindowTerm,
                       OffsetValueType            parzenWindowIndex,
                       const KernelFunctionType & kernel,
                       PDFValueType *             parzenValues);

  virtual void
  UpdateJointPDFAndDerivatives(const RealType                     fixedImageValue,
                               const RealType                     movingImageValue,
                               const DerivativeType *             imageJacobian,
                               const NonZeroJacobianIndicesType * nzji,
                               JointPDFType *                     jointPDF) const;

  virtual void
  UpdateJointPDFAndIncrementalPDFs(RealType                           fixedImageValue,
                                   RealType                           movingImageValue,
                                   RealType                           movingMaskValue,
                                   const DerivativeType &             movingImageValuesRight,
                                   const DerivativeType &             movingImageValuesLeft,
                                   const DerivativeType &             movingMaskValuesRight,
                                   const DerivativeType &             movingMaskValuesLeft,
                                   const NonZeroJacobianIndicesType & nzji) const;

  void
  UpdateJointPDFDerivatives(const JointPDFIndexType &          pdfIndex,
                            double                             factor,
                            const DerivativeType &             imageJacobian,
                            const NonZeroJacobianIndicesType & nzji) const;

  void
  NormalizeJointPDF(JointPDFType * pdf, const double factor) const;

  void
  NormalizeJointPDFDerivatives(JointPDFDerivativesType * pdf, const double factor) const;

  void
  ComputeMarginalPDF(const JointPDFType * jointPDF, MarginalPDFType & marginalPDF, const unsigned int direction) const;

  virtual void
  ComputeIncrementalMarginalPDFs(const JointPDFDerivativesType * incrementalPDF,
                                 IncrementalMarginalPDFType *    fixedIncrementalMarginalPDF,
                                 IncrementalMarginalPDFType *    movingIncrementalMarginalPDF) const;

  virtual void
  ComputePDFsAndPDFDerivatives(const ParametersType & parameters) const;

  virtual void
  ComputePDFsAndIncrementalPDFs(const ParametersType & parameters) const;

  virtual void
  ComputePDFsSingleThreaded(const ParametersType & parameters) const;
 
  virtual void
  ComputePDFs(const ParametersType & parameters) const;

  virtual void
  InitializeHistograms();
 
  virtual void
  InitializeKernels();

  virtual void
  GetValueAndAnalyticDerivative(const ParametersType & itkNotUsed(parameters),
                                MeasureType &          itkNotUsed(value),
                                DerivativeType &       itkNotUsed(derivative)) const
  {}

  virtual void
  GetValueAndFiniteDifferenceDerivative(const ParametersType & itkNotUsed(parameters),
                                        MeasureType &          itkNotUsed(value),
                                        DerivativeType &       itkNotUsed(derivative)) const
  {}
 
private:
  mutable std::vector<JointPDFPointer> m_ThreaderJointPDFs{};

  struct ParzenWindowHistogramMultiThreaderParameterType // can't we use the one from AdvancedImageToImageMetric ?
  {
    Self * m_Metric;
  };
  ParzenWindowHistogramMultiThreaderParameterType m_ParzenWindowHistogramThreaderParameters{};
 
  struct ParzenWindowHistogramGetValueAndDerivativePerThreadStruct
  {
    SizeValueType   st_NumberOfPixelsCounted;
    JointPDFPointer st_JointPDF;
  };
  itkPadStruct(ITK_CACHE_LINE_ALIGNMENT,
               ParzenWindowHistogramGetValueAndDerivativePerThreadStruct,
               PaddedParzenWindowHistogramGetValueAndDerivativePerThreadStruct);
  itkAlignedTypedef(ITK_CACHE_LINE_ALIGNMENT,
                    PaddedParzenWindowHistogramGetValueAndDerivativePerThreadStruct,
                    AlignedParzenWindowHistogramGetValueAndDerivativePerThreadStruct);
  mutable std::vector<AlignedParzenWindowHistogramGetValueAndDerivativePerThreadStruct>
    m_ParzenWindowHistogramGetValueAndDerivativePerThreadVariables;

  unsigned long m_NumberOfFixedHistogramBins{ 32 };
  unsigned long m_NumberOfMovingHistogramBins{ 32 };
  unsigned int  m_FixedKernelBSplineOrder{ 0 };
  unsigned int  m_MovingKernelBSplineOrder{ 3 };
  bool          m_UseDerivative{ false };
  bool          m_UseExplicitPDFDerivatives{ true };
  bool          m_UseFiniteDifferenceDerivative{ false };
  double        m_FiniteDifferencePerturbation{ 1.0 };
};
 
} // end namespace itk
 
#ifndef ITK_MANUAL_INSTANTIATION
#  include "itkParzenWindowHistogramImageToImageMetric.hxx"
#endif
 
#endif // end #ifndef itkParzenWindowHistogramImageToImageMetric_h