itkCRImageToImageMetric.h

Go to the documentation of this file.

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

  NifTK: A software platform for medical image computing.

  Copyright (c) University College London (UCL). All rights reserved.

  This software is distributed WITHOUT ANY WARRANTY; without even
  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  PURPOSE.

  See LICENSE.txt in the top level directory for details.

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

#ifndef itkCRImageToImageMetric_h
#define itkCRImageToImageMetric_h

#include "itkFiniteDifferenceGradientSimilarityMeasure.h"

namespace itk
{
template < class TFixedImage, class TMovingImage > 
class ITK_EXPORT CRImageToImageMetric : 
    public FiniteDifferenceGradientSimilarityMeasure< TFixedImage, TMovingImage>
{
public:

  typedef CRImageToImageMetric                                     Self;
  typedef SimilarityMeasure<TFixedImage, TMovingImage >            Superclass;
  typedef SmartPointer<Self>                                       Pointer;
  typedef SmartPointer<const Self>                                 ConstPointer;
  typedef typename Superclass::FixedImageType::PixelType           FixedImagePixelType;
  typedef typename Superclass::MovingImageType::PixelType          MovingImagePixelType;  
  typedef typename Superclass::MeasureType                         MeasureType;
  typedef std::multimap<FixedImagePixelType, MovingImagePixelType> MapType;
  typedef typename MapType::iterator                               MapIterator;
  
  itkNewMacro(Self);
 
  itkTypeMacro(CRImageToImageMetric, SimilarityMeasure);

  bool ShouldBeMaximized() { return true; };

protected:
  
  CRImageToImageMetric() {};
  virtual ~CRImageToImageMetric() {};

  void ResetCostFunction() 
    { 
      this->m_NumberCounted = 0;
      this->m_MovingSum = 0;
      this->m_Map.clear();      
    }
  
  void AggregateCostFunctionPair(
      FixedImagePixelType fixedValue, 
      MovingImagePixelType movingValue)
    {
      this->m_NumberCounted++;
      this->m_MovingSum += movingValue;
      this->m_Map.insert(std::pair<FixedImagePixelType, MovingImagePixelType>(fixedValue, movingValue));
    }
  
  MeasureType FinalizeCostFunction()
    {
      MeasureType measure = 0;

      double totalMean = 0;
      double totalSigmaSquared = 0;
      
      double conditionalNumber = 0;
      double conditionalSum = 0;
      double conditionalMean = 0;
      double conditionalSigmaSquared = 0;
      
      FixedImagePixelType fixedValue;
      MovingImagePixelType movingValue;
      
      MapIterator mapIterator;
      
      totalMean = this->m_MovingSum / (double)this->m_NumberCounted;
      
      mapIterator   = this->m_Map.begin();
      while(mapIterator != this->m_Map.end())
        {
          fixedValue = (*mapIterator).first;
          
          conditionalNumber = 0;
          conditionalSum = 0;
          conditionalSigmaSquared = 0;
          
          while( (mapIterator != this->m_Map.end()) && ((*mapIterator).first == fixedValue) ) 
          {
            movingValue = (*mapIterator).second; 
              
            conditionalSum          += movingValue;
            conditionalSigmaSquared += movingValue*movingValue;
            totalSigmaSquared       += movingValue*movingValue;
            mapIterator++;
            conditionalNumber++;
          } 
             
          conditionalMean = conditionalSum / conditionalNumber;
          
          conditionalSigmaSquared = (conditionalSigmaSquared / conditionalNumber) - (conditionalMean*conditionalMean);
          
          measure += (conditionalNumber * conditionalSigmaSquared); 
          
        }
      
      totalSigmaSquared  = (totalSigmaSquared / (double)this->m_NumberCounted) - totalMean;
      
      if (totalSigmaSquared == 0)
        {
          measure = 0;
        }
      else
        {
          measure /= ((double)this->m_NumberCounted * totalSigmaSquared);    
        }

      return 1.0-measure;
    }

private:
  CRImageToImageMetric(const Self&);  // purposefully not implemented
  void operator=(const Self&);        // purposefully not implemented
  
  unsigned long int m_NumberCounted;
  double   m_MovingSum;
  MapType m_Map;
  
};

} // end namespace itk

#endif