itkMattesMutualInformationImageToImageMetric.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 itkMattesMutualInformationImageToImageMetric_h
#define itkMattesMutualInformationImageToImageMetric_h

#include "itkImageToImageMetric.h"
#include "itkPoint.h"
#include "itkIndex.h"
#include "itkBSplineDerivativeKernelFunction.h"
#include "itkArray2D.h"

#include <mutex>


namespace itk
{
template <typename TFixedImage, typename TMovingImage>
class ITK_TEMPLATE_EXPORT MattesMutualInformationImageToImageMetric:
  public ImageToImageMetric<TFixedImage, TMovingImage>
{
public:
  ITK_DISALLOW_COPY_AND_ASSIGN(MattesMutualInformationImageToImageMetric);

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

  itkNewMacro(Self);

  itkTypeMacro(MattesMutualInformationImageToImageMetric,
               ImageToImageMetric);

  using TransformType = typename Superclass::TransformType;
  using TransformPointer = typename Superclass::TransformPointer;
  using TransformJacobianType = typename Superclass::TransformJacobianType;
  using InterpolatorType = typename Superclass::InterpolatorType;
  using MeasureType = typename Superclass::MeasureType;
  using DerivativeType = typename Superclass::DerivativeType;
  using ParametersType = typename Superclass::ParametersType;
  using FixedImageType = typename Superclass::FixedImageType;
  using MovingImageType = typename Superclass::MovingImageType;
  using MovingImagePointType = typename Superclass::MovingImagePointType;
  using FixedImageConstPointer = typename Superclass::FixedImageConstPointer;
  using MovingImageConstPointer = typename Superclass::MovingImageConstPointer;
  using BSplineTransformWeightsType = typename Superclass::BSplineTransformWeightsType;
  using BSplineTransformIndexArrayType = typename Superclass::BSplineTransformIndexArrayType;

  using CoordinateRepresentationType = typename Superclass::CoordinateRepresentationType;
  using FixedImageSampleContainer = typename Superclass::FixedImageSampleContainer;
  using ImageDerivativesType = typename Superclass::ImageDerivativesType;
  using WeightsValueType = typename Superclass::WeightsValueType;
  using IndexValueType = typename Superclass::IndexValueType;

  using OffsetValueType = typename FixedImageType::OffsetValueType;

  static constexpr unsigned int MovingImageDimension = MovingImageType::ImageDimension;

  void Initialize() override;

  MeasureType GetValue(const ParametersType & parameters) const override;

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

  void GetValueAndDerivative(const ParametersType & parameters, MeasureType & Value, DerivativeType & Derivative) const override;

  itkSetClampMacro( NumberOfHistogramBins, SizeValueType,
                    5, NumericTraits<SizeValueType>::max() );
  itkGetConstReferenceMacro(NumberOfHistogramBins, SizeValueType);

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

  using PDFValueType = double; //NOTE:  floating point precision is not as stable.  Double precision proves faster and more robust in real-world testing.

  using JointPDFType = Image<PDFValueType, 2>;
  using JointPDFDerivativesType = Image<PDFValueType, 3>;

  const typename JointPDFType::Pointer GetJointPDF () const
    {
    if( this->m_MMIMetricPerThreadVariables == nullptr )
      {
      return JointPDFType::Pointer(nullptr);
      }
    return this->m_MMIMetricPerThreadVariables[0].JointPDF;
    }

  const typename JointPDFDerivativesType::Pointer GetJointPDFDerivatives () const
    {
    if( this->m_MMIMetricPerThreadVariables == nullptr )
      {
      return JointPDFDerivativesType::Pointer(nullptr);
      }
    return this->m_MMIMetricPerThreadVariables[0].JointPDFDerivatives;
    }

protected:

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

private:
  using JointPDFIndexType = JointPDFType::IndexType;
  using JointPDFValueType = JointPDFType::PixelType;
  using JointPDFRegionType = JointPDFType::RegionType;
  using JointPDFSizeType = JointPDFType::SizeType;
  using JointPDFDerivativesIndexType = JointPDFDerivativesType::IndexType;
  using JointPDFDerivativesValueType = JointPDFDerivativesType::PixelType;
  using JointPDFDerivativesRegionType = JointPDFDerivativesType::RegionType;
  using JointPDFDerivativesSizeType = JointPDFDerivativesType::SizeType;

  using CubicBSplineFunctionType = BSplineKernelFunction<3,PDFValueType>;
  using CubicBSplineDerivativeFunctionType = BSplineDerivativeKernelFunction<3,PDFValueType>;

  void ComputeFixedImageParzenWindowIndices( FixedImageSampleContainer & samples);

  void ComputePDFDerivatives(ThreadIdType threadId, unsigned int sampleNumber, int movingImageParzenWindowIndex,
                                     const ImageDerivativesType
                                     &  movingImageGradientValue,
                                     PDFValueType cubicBSplineDerivativeValue) const;

  void GetValueThreadPreProcess(ThreadIdType threadId, bool withinSampleThread) const override;
  void GetValueThreadPostProcess(ThreadIdType threadId, bool withinSampleThread) const override;
  //NOTE:  The signature in base class requires that movingImageValue is of type double
  bool GetValueThreadProcessSample(ThreadIdType threadId, SizeValueType fixedImageSample,
                                                  const MovingImagePointType & mappedPoint,
                                                  double movingImageValue) const override;

  void GetValueAndDerivativeThreadPreProcess( ThreadIdType threadId, bool withinSampleThread) const override;
  void GetValueAndDerivativeThreadPostProcess( ThreadIdType threadId, bool withinSampleThread) const override;
  //NOTE:  The signature in base class requires that movingImageValue is of type double
  bool GetValueAndDerivativeThreadProcessSample(ThreadIdType threadId, SizeValueType fixedImageSample,
                                                               const MovingImagePointType & mappedPoint,
                                                               double movingImageValue, const ImageDerivativesType &
                                                               movingImageGradientValue) const override;

  SizeValueType m_NumberOfHistogramBins{50};
  PDFValueType  m_MovingImageNormalizedMin{0.0};
  PDFValueType  m_FixedImageNormalizedMin{0.0};
  PDFValueType  m_FixedImageTrueMin{0.0};
  PDFValueType  m_FixedImageTrueMax{0.0};
  PDFValueType  m_MovingImageTrueMin{0.0};
  PDFValueType  m_MovingImageTrueMax{0.0};
  PDFValueType  m_FixedImageBinSize{0.0};
  PDFValueType  m_MovingImageBinSize{0.0};

  typename CubicBSplineFunctionType::Pointer           m_CubicBSplineKernel;
  typename CubicBSplineDerivativeFunctionType::Pointer m_CubicBSplineDerivativeKernel;

  using PRatioType = PDFValueType;
  using PRatioArrayType = Array2D<PRatioType>;

  mutable PRatioArrayType m_PRatioArray;

  using MarginalPDFType = std::vector< PDFValueType >;
  mutable MarginalPDFType             m_MovingImageMarginalPDF;

  struct MMIMetricPerThreadStruct
  {
    int JointPDFStartBin;
    int JointPDFEndBin;

    PDFValueType JointPDFSum;

    DerivativeType MetricDerivative;

    typename JointPDFType::Pointer            JointPDF;
    typename JointPDFDerivativesType::Pointer JointPDFDerivatives;

    typename TransformType::JacobianType Jacobian;

    MarginalPDFType FixedImageMarginalPDF;
  };

#if !defined(ITK_WRAPPING_PARSER)
  itkPadStruct( ITK_CACHE_LINE_ALIGNMENT, MMIMetricPerThreadStruct,
                                            PaddedMMIMetricPerThreadStruct);
  itkAlignedTypedef( ITK_CACHE_LINE_ALIGNMENT, PaddedMMIMetricPerThreadStruct,
                                               AlignedMMIMetricPerThreadStruct );
  // Due to a bug in older version of Visual Studio where std::vector resize
  // uses a value instead of a const reference, this must be a pointer.
  // See
  //   http://thetweaker.wordpress.com/2010/05/05/stdvector-of-aligned-elements/
  //   http://connect.microsoft.com/VisualStudio/feedback/details/692988
  mutable AlignedMMIMetricPerThreadStruct * m_MMIMetricPerThreadVariables;
#endif

  bool         m_UseExplicitPDFDerivatives{true};
  mutable bool m_ImplicitDerivativesSecondPass{false};
};
} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkMattesMutualInformationImageToImageMetric.hxx"
#endif

#endif