itkOptImageToImageMetric.h

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkOptImageToImageMetric.h,v $
  Language:  C++
  Date:      $Date: 2008-07-30 20:56:17 $
  Version:   $Revision: 1.22 $

  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 "itkSingleValuedCostFunction.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 );

  itkGetConstReferenceMacro( NumberOfMovingImageSamples, unsigned long );
  unsigned long GetNumberOfPixelsCounter( void )
    {
    return GetNumberOfMovingImageSamples(); 
    }

  itkSetMacro( FixedImageRegion, FixedImageRegionType );

  itkGetConstReferenceMacro( FixedImageRegion, FixedImageRegionType );

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

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

  void SetFixedImageIndexes( const FixedImageIndexContainer & indexes );

  itkSetMacro( NumberOfThreads, unsigned int );
  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 )
    {
    // This is not a macro because we need to also set m_UseAllPixels
    itkDebugMacro("Setting NumberOfFixedImageSamples to " << numSamples ); 
    if (this->m_NumberOfFixedImageSamples != numSamples)
      {
      if( this->m_NumberOfFixedImageSamples == 0 )
        {
        // make sure that there is at least one sample
        this->m_NumberOfFixedImageSamples = 1;
        }
      this->m_UseAllPixels = false;
      this->m_NumberOfFixedImageSamples = numSamples;
      this->Modified();
      } 
    } 

  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 );
  itkSetMacro( UseFixedImageSamplesIntensityThreshold, bool );
  itkGetConstReferenceMacro( UseFixedImageSamplesIntensityThreshold, bool );

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

  itkGetConstReferenceMacro( NumberOfPixelsCounted, unsigned long );

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

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

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

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

  mutable unsigned long       m_NumberOfPixelsCounted;


  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;
  FixedImagePixelType       m_FixedImageSamplesIntensityThreshold;
  bool                      m_UseFixedImageSamplesIntensityThreshold;

  typedef std::vector<FixedImageSamplePoint> FixedImageSampleContainer;

  virtual void SampleFixedImageDomain( FixedImageSampleContainer & samples) const;

  virtual void SampleFixedImageIndexes( FixedImageSampleContainer & 
                                           samples);

  virtual void SampleFullFixedImageDomain( FixedImageSampleContainer & 
                                           samples);

  FixedImageSampleContainer   m_FixedImageSamples;

  unsigned long               m_NumberOfParameters;
  mutable ParametersType      m_Parameters;

  mutable unsigned long       m_NumberOfFixedImageSamples;
  mutable unsigned long       m_NumberOfMovingImageSamples;

  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_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;


  typedef MultiThreader               MultiThreaderType;

  struct MultiThreaderParameterType
    {
    ImageToImageMetric                   * metric;
    };

  MultiThreaderType::Pointer               m_Threader;
  MultiThreaderParameterType               m_ThreaderParameter;
  mutable unsigned int                     m_ThreaderChunkSize;
  mutable unsigned int                     m_ThreaderSizeOfLastChunk;
  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 );

  void                      GetValueThread( unsigned int threadID ) const;
  virtual inline void       GetValueThreadPreProcess( 
                                       unsigned int threadID,
                                       bool 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                          InitiateGetDerivativeMultiThreaded( void) const;
  static ITK_THREAD_RETURN_TYPE GetDerivativeMultiThreaded( void * arg );
  void                       GetDerivativeThread( unsigned int threadID ) const;
  virtual inline void   GetDerivativeThreadProcessSample( unsigned int threadID,
                                             unsigned long fixedImageSample,
                                             MovingImagePointType mappedPoint,
                                             double movingImageValue) const = 0;
  */

  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);

  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 {};

  void SynchronizeTransforms() const;

  void NumberOfFixedImageSamplesUpdated();

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