itkOptImageToImageMetric.h

Go to the documentation of this file.

/*=========================================================================

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkOptImageToImageMetric.h,v $
  Language:  C++
  Date:      $Date: 2010-03-05 16:15:37 $
  Version:   $Revision: 1.34 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notices for more information.

=========================================================================*/
#ifndef __itkOptImageToImageMetric_h
#define __itkOptImageToImageMetric_h

#include "itkSingleValuedCostFunction.h"
#include "itkImageBase.h"
#include "itkTransform.h"
#include "itkInterpolateImageFunction.h"
#include "itkExceptionObject.h"
#include "itkGradientRecursiveGaussianImageFilter.h"
#include "itkSpatialObject.h"
#include "itkBSplineDeformableTransform.h"
#include "itkCentralDifferenceImageFunction.h"
#include "itkCovariantVector.h"

#include "itkMultiThreader.h"

#include "itkOptBSplineInterpolateImageFunction.h"

namespace itk
{
  
template <class TFixedImage,  class TMovingImage> 
class ITK_EXPORT ImageToImageMetric 
: public SingleValuedCostFunction 
{
public:
  typedef ImageToImageMetric           Self;
  typedef SingleValuedCostFunction     Superclass;
  typedef SmartPointer<Self>           Pointer;
  typedef SmartPointer<const Self>     ConstPointer;

  typedef typename Superclass::ParametersValueType CoordinateRepresentationType;

  itkTypeMacro(ImageToImageMetric, SingleValuedCostFunction);

  typedef TMovingImage                               MovingImageType;
  typedef typename TMovingImage::PixelType           MovingImagePixelType;
  typedef typename MovingImageType::ConstPointer     MovingImageConstPointer;

  typedef TFixedImage                                FixedImageType;
  typedef typename TFixedImage::PixelType            FixedImagePixelType;
  typedef typename FixedImageType::ConstPointer      FixedImageConstPointer;
  typedef typename FixedImageType::RegionType        FixedImageRegionType;

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

  typedef Transform<CoordinateRepresentationType, 
                    itkGetStaticConstMacro(MovingImageDimension),
                    itkGetStaticConstMacro(FixedImageDimension)> 
                                                     TransformType;

  typedef typename TransformType::Pointer            TransformPointer;
  typedef typename TransformType::InputPointType     InputPointType;
  typedef typename TransformType::OutputPointType    OutputPointType;
  typedef typename TransformType::ParametersType     TransformParametersType;
  typedef typename TransformType::JacobianType       TransformJacobianType;

  typedef typename FixedImageType::IndexType           FixedImageIndexType;
  typedef typename FixedImageIndexType::IndexValueType FixedImageIndexValueType;
  typedef typename MovingImageType::IndexType          MovingImageIndexType;
  typedef typename TransformType::InputPointType       FixedImagePointType;
  typedef typename TransformType::OutputPointType      MovingImagePointType;

  typedef std::vector<FixedImageIndexType>             FixedImageIndexContainer;

  typedef InterpolateImageFunction< MovingImageType, 
                                    CoordinateRepresentationType >
                                                     InterpolatorType;

  typedef typename NumericTraits<MovingImagePixelType>::RealType 
                                                     RealType;
  typedef CovariantVector<RealType,
                          itkGetStaticConstMacro(MovingImageDimension)>
                                                     GradientPixelType;
  typedef Image<GradientPixelType,
                itkGetStaticConstMacro(MovingImageDimension)> 
                                                     GradientImageType;
  typedef SmartPointer<GradientImageType>            GradientImagePointer;
  typedef GradientRecursiveGaussianImageFilter< MovingImageType,
                                                GradientImageType >
                                                     GradientImageFilterType;  
  typedef typename GradientImageFilterType::Pointer  GradientImageFilterPointer;


  typedef typename InterpolatorType::Pointer         InterpolatorPointer;


  typedef SpatialObject< itkGetStaticConstMacro(FixedImageDimension) >
                                                     FixedImageMaskType;
  typedef typename FixedImageMaskType::Pointer       FixedImageMaskPointer;
  typedef typename FixedImageMaskType::ConstPointer  FixedImageMaskConstPointer;

  typedef SpatialObject< itkGetStaticConstMacro(MovingImageDimension) >
                                                     MovingImageMaskType;
  typedef typename MovingImageMaskType::Pointer      MovingImageMaskPointer;
  typedef typename MovingImageMaskType::ConstPointer MovingImageMaskConstPointer;


  typedef typename Superclass::MeasureType                    MeasureType;

  typedef typename Superclass::DerivativeType                 DerivativeType;

  typedef typename Superclass::ParametersType                 ParametersType;

  itkSetConstObjectMacro( FixedImage, FixedImageType );

  itkGetConstObjectMacro( FixedImage, FixedImageType );

  itkSetConstObjectMacro( MovingImage, MovingImageType );

  itkGetConstObjectMacro( MovingImage, MovingImageType );

  itkSetObjectMacro( Transform, TransformType );

  itkGetConstObjectMacro( Transform, TransformType );

  itkSetObjectMacro( Interpolator, InterpolatorType );

  itkGetConstObjectMacro( Interpolator, InterpolatorType );

  unsigned long GetNumberOfMovingImageSamples( void )
    {
    return this->GetNumberOfPixelsCounted();
    }
  itkGetConstReferenceMacro( NumberOfPixelsCounted, unsigned long );

  void SetFixedImageRegion( const FixedImageRegionType reg );

  itkGetConstReferenceMacro( FixedImageRegion, FixedImageRegionType );

  itkSetObjectMacro( MovingImageMask, MovingImageMaskType );
  itkSetConstObjectMacro( MovingImageMask, MovingImageMaskType );
  itkGetConstObjectMacro( MovingImageMask, MovingImageMaskType );

  itkSetObjectMacro( FixedImageMask, FixedImageMaskType );
  itkSetConstObjectMacro( FixedImageMask, FixedImageMaskType );
  itkGetConstObjectMacro( FixedImageMask, FixedImageMaskType );

  void SetFixedImageIndexes( const FixedImageIndexContainer & indexes );
  void SetUseFixedImageIndexes( bool useIndex );
  itkGetConstReferenceMacro( UseFixedImageIndexes, bool );

  void SetNumberOfThreads( unsigned int numberOfThreads );
  itkGetConstReferenceMacro( NumberOfThreads, unsigned int );

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

  virtual void ComputeGradient( void );

  itkGetConstObjectMacro( GradientImage, GradientImageType );

  void SetTransformParameters( const ParametersType & parameters ) const;

  unsigned int GetNumberOfParameters( void ) const 
    {
    return m_Transform->GetNumberOfParameters(); 
    }

  virtual void Initialize( void ) throw ( ExceptionObject );

  virtual void MultiThreadingInitialize( void ) throw ( ExceptionObject );

  virtual void SetNumberOfFixedImageSamples( unsigned long numSamples );
  itkGetConstReferenceMacro( NumberOfFixedImageSamples, unsigned long ); 

  void SetNumberOfSpatialSamples( unsigned long num )
    {
    this->SetNumberOfFixedImageSamples( num ); 
    }
  unsigned long GetNumberOfSpatialSamples( void )
    {
    return this->GetNumberOfFixedImageSamples(); 
    }

  void SetFixedImageSamplesIntensityThreshold( const FixedImagePixelType & thresh );
  itkGetConstReferenceMacro( FixedImageSamplesIntensityThreshold, FixedImagePixelType );

  void SetUseFixedImageSamplesIntensityThreshold( bool useThresh );
  itkGetConstReferenceMacro( UseFixedImageSamplesIntensityThreshold, bool );

  void SetUseAllPixels( bool useAllPixels );
  void UseAllPixelsOn( void )
    {
    this->SetUseAllPixels( true );
    }
  void UseAllPixelsOff( void )
    {
    this->SetUseAllPixels( false );
    }
  itkGetConstReferenceMacro( UseAllPixels, bool );

  void SetUseSequentialSampling( bool sequentialSampling );
  itkGetConstReferenceMacro( UseSequentialSampling, bool );

  void ReinitializeSeed();
  void ReinitializeSeed( int seed );  

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

  typedef MultiThreader               MultiThreaderType;
  itkGetConstObjectMacro( Threader, MultiThreaderType );
  const TransformPointer* GetThreaderTransform()
    {
    return m_ThreaderTransform;
    }

protected:
  ImageToImageMetric();
  virtual ~ImageToImageMetric();

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


  class FixedImageSamplePoint
    {
    public:
    FixedImageSamplePoint()
      { 
      point.Fill(0.0);
      value = 0;
      valueIndex = 0;
      }
    ~FixedImageSamplePoint() {};

    public:
      FixedImagePointType           point;
      double                        value;
      unsigned int                  valueIndex;
    };
  
  bool                      m_UseFixedImageIndexes;
  FixedImageIndexContainer  m_FixedImageIndexes;

  bool                      m_UseFixedImageSamplesIntensityThreshold;
  FixedImagePixelType       m_FixedImageSamplesIntensityThreshold;

  typedef std::vector<FixedImageSamplePoint> FixedImageSampleContainer;

  virtual void SampleFixedImageRegion( FixedImageSampleContainer & samples) const;

  virtual void SampleFixedImageIndexes( FixedImageSampleContainer & 
                                           samples) const;

  virtual void SampleFullFixedImageRegion( FixedImageSampleContainer & 
                                           samples) const;

  FixedImageSampleContainer   m_FixedImageSamples;

  unsigned long               m_NumberOfParameters;
  mutable ParametersType      m_Parameters;

  unsigned long               m_NumberOfFixedImageSamples;
  //m_NumberOfPixelsCounted must be mutable because the const
  //thread consolidation functions merge each threads valus
  //onto this accumulator variable.
  mutable unsigned long       m_NumberOfPixelsCounted;

  FixedImageConstPointer      m_FixedImage;
  MovingImageConstPointer     m_MovingImage;

  TransformPointer            m_Transform;

  TransformPointer          * m_ThreaderTransform; 

  InterpolatorPointer         m_Interpolator;

  bool                        m_ComputeGradient;
  GradientImagePointer        m_GradientImage;

  FixedImageMaskConstPointer  m_FixedImageMask;
  MovingImageMaskConstPointer m_MovingImageMask;

  unsigned int                m_NumberOfThreads;

  bool                        m_UseAllPixels;
  bool                        m_UseSequentialSampling;

  bool                        m_ReseedIterator;

  int                         m_RandomSeed;

  bool                        m_TransformIsBSpline;

  unsigned long               m_NumBSplineWeights;

  itkStaticConstMacro(DeformationSplineOrder, unsigned int, 3 );

  typedef BSplineDeformableTransform< CoordinateRepresentationType,
                      ::itk::GetImageDimension<FixedImageType>::ImageDimension,
                                      itkGetStaticConstMacro(DeformationSplineOrder) >             BSplineTransformType;

  typedef typename BSplineTransformType::WeightsType       BSplineTransformWeightsType;
  typedef typename BSplineTransformWeightsType::ValueType  WeightsValueType;
  typedef          Array2D<WeightsValueType>               BSplineTransformWeightsArrayType;

  typedef typename BSplineTransformType::ParameterIndexArrayType 
                                                              BSplineTransformIndexArrayType;
  typedef typename BSplineTransformIndexArrayType::ValueType  IndexValueType;
  typedef          Array2D<IndexValueType>                    BSplineTransformIndicesArrayType;

  typedef          std::vector<MovingImagePointType>          MovingImagePointArrayType;
  typedef          std::vector<bool>                          BooleanArrayType;
  typedef          FixedArray< unsigned long, 
                               ::itk::GetImageDimension<FixedImageType>
                                ::ImageDimension >            BSplineParametersOffsetType;
  typedef BSplineInterpolateImageFunction<MovingImageType,
                                             CoordinateRepresentationType> 
                                                       BSplineInterpolatorType;

  typedef CentralDifferenceImageFunction<MovingImageType,
                                         CoordinateRepresentationType> 
                                                       DerivativeFunctionType;
  typedef          CovariantVector< double,
                           itkGetStaticConstMacro(MovingImageDimension) >
                                                          ImageDerivativesType;


  typename BSplineTransformType::Pointer   m_BSplineTransform;

  BSplineTransformWeightsArrayType         m_BSplineTransformWeightsArray;
  BSplineTransformIndicesArrayType         m_BSplineTransformIndicesArray;
  MovingImagePointArrayType                m_BSplinePreTransformPointsArray;
  BooleanArrayType                         m_WithinBSplineSupportRegionArray;
  
  BSplineParametersOffsetType              m_BSplineParametersOffset;

  // Variables needed for optionally caching values when using a BSpline transform.
  bool                                    m_UseCachingOfBSplineWeights;
  mutable BSplineTransformWeightsType     m_BSplineTransformWeights;
  mutable BSplineTransformIndexArrayType  m_BSplineTransformIndices;

  mutable BSplineTransformWeightsType     * m_ThreaderBSplineTransformWeights;
  mutable BSplineTransformIndexArrayType  * m_ThreaderBSplineTransformIndices;

  virtual void PreComputeTransformValues( void );

  virtual void TransformPoint( unsigned int sampleNumber,
                               MovingImagePointType& mappedPoint,
                               bool& sampleWithinSupportRegion,
                               double& movingImageValue,
                               unsigned int threadID ) const;

  virtual void TransformPointWithDerivatives( unsigned int sampleNumber,
                               MovingImagePointType& mappedPoint,
                               bool& sampleWithinSupportRegion,
                               double& movingImageValue,
                               ImageDerivativesType & gradient,
                               unsigned int threadID ) const;

  bool                                       m_InterpolatorIsBSpline;

  typename BSplineInterpolatorType::Pointer  m_BSplineInterpolator;

  typename DerivativeFunctionType::Pointer   m_DerivativeCalculator;

  virtual void ComputeImageDerivatives( const MovingImagePointType & mappedPoint,
                                        ImageDerivativesType & gradient,
                                        unsigned int threadID ) const;


  struct MultiThreaderParameterType
    {
    ImageToImageMetric                   * metric;
    };

  MultiThreaderType::Pointer               m_Threader;
  MultiThreaderParameterType               m_ThreaderParameter;
  mutable unsigned int                   * m_ThreaderNumberOfMovingImageSamples;
  bool                                     m_WithinThreadPreProcess;
  bool                                     m_WithinThreadPostProcess;

  void                           GetValueMultiThreadedPreProcessInitiate( void ) const;
  void                           GetValueMultiThreadedInitiate( void ) const;
  void                           GetValueMultiThreadedPostProcessInitiate( void ) const;
  static ITK_THREAD_RETURN_TYPE  GetValueMultiThreadedPreProcess( void * arg );
  static ITK_THREAD_RETURN_TYPE  GetValueMultiThreaded( void * arg );
  static ITK_THREAD_RETURN_TYPE  GetValueMultiThreadedPostProcess( void * arg );

  virtual inline void       GetValueThread( unsigned int threadID ) const;
  virtual inline void       GetValueThreadPreProcess( 
                                 unsigned int itkNotUsed(threadID),
                                 bool itkNotUsed(withinSampleThread) ) const
    { };
  virtual inline bool       GetValueThreadProcessSample( 
                                 unsigned int itkNotUsed(threadID),
                                 unsigned long itkNotUsed(fixedImageSample),
                                 const MovingImagePointType & itkNotUsed(mappedPoint),
                                 double itkNotUsed(movingImageValue)) const
    { return false; };
  virtual inline void       GetValueThreadPostProcess( 
                                 unsigned int itkNotUsed(threadID),
                                 bool itkNotUsed(withinSampleThread) ) const
    { };
      
  void                          GetValueAndDerivativeMultiThreadedPreProcessInitiate( void ) const;
  void                          GetValueAndDerivativeMultiThreadedInitiate( void ) const;
  void                          GetValueAndDerivativeMultiThreadedPostProcessInitiate( void ) const;
  static ITK_THREAD_RETURN_TYPE GetValueAndDerivativeMultiThreadedPreProcess( void * arg );
  static ITK_THREAD_RETURN_TYPE GetValueAndDerivativeMultiThreaded(void * arg);
  static ITK_THREAD_RETURN_TYPE GetValueAndDerivativeMultiThreadedPostProcess(void * arg);

  virtual inline void  GetValueAndDerivativeThread(unsigned int threadID) const;
  virtual inline void  GetValueAndDerivativeThreadPreProcess(
                                unsigned int itkNotUsed(threadID),
                                bool itkNotUsed(withinSampleThread)) const 
    { };
  virtual inline bool  GetValueAndDerivativeThreadProcessSample( 
                                unsigned int itkNotUsed(threadID),
                                unsigned long itkNotUsed(fixedImageSample),
                                const MovingImagePointType & itkNotUsed(mappedPoint),
                                double itkNotUsed(movingImageValue),
                                const ImageDerivativesType & itkNotUsed(movingImageGradientValue) ) const
    { return false; }
  virtual inline void  GetValueAndDerivativeThreadPostProcess( 
                                unsigned int itkNotUsed(threadID),
                                bool itkNotUsed(withinSampleThread) ) const 
    { };

  virtual void SynchronizeTransforms() const;

private:
  ImageToImageMetric(const Self&); //purposely not implemented
  void operator=(const Self&); //purposely not implemented
  
  FixedImageRegionType        m_FixedImageRegion;  

};

} // end namespace itk

#ifndef ITK_MANUAL_INSTANTIATION
#include "itkOptImageToImageMetric.txx"
#endif

#endif