itkPIUImageToImageMetric.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 itkPIUImageToImageMetric_h
#define itkPIUImageToImageMetric_h

#include "itkHistogramSimilarityMeasure.h"

namespace itk
{
template < class TFixedImage, class TMovingImage > 
class ITK_EXPORT PIUImageToImageMetric : 
    public HistogramSimilarityMeasure< TFixedImage, TMovingImage>
{
public:

  typedef PIUImageToImageMetric                                           Self;
  typedef HistogramSimilarityMeasure<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 typename Superclass::HistogramType                              HistogramType;
  typedef typename HistogramType::FrequencyType                           HistogramFrequencyType;
  typedef typename HistogramType::Iterator                                HistogramIteratorType;
  typedef typename HistogramType::MeasurementVectorType                   HistogramMeasurementVectorType;
  typedef typename HistogramType::IndexType                               HistogramIndexType;
  typedef typename HistogramType::SizeValueType                           HistogramSizeValueType;

  itkNewMacro(Self);
 
  itkTypeMacro(PIUImageToImageMetric, HistogramSimilarityMeasure);

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

  MeasureType FinalizeCostFunction()
    {
      MeasureType piu = NumericTraits<MeasureType>::Zero;
      HistogramFrequencyType movingFrequency;
      HistogramIndexType index(2);
       
      HistogramSizeValueType sf = this->m_Histogram->GetSize(0);
      HistogramSizeValueType sm = this->m_Histogram->GetSize(1);
       
      HistogramFrequencyType totalFrequency = this->m_Histogram->GetTotalFrequency();

      for (unsigned int f = 0; f < sf; f++)
        {
          // Reset to zero.
          MeasureType value = NumericTraits<MeasureType>::Zero;
          MeasureType mean = NumericTraits<MeasureType>::Zero;
          MeasureType sum = NumericTraits<MeasureType>::Zero;
          MeasureType sumSquared = NumericTraits<MeasureType>::Zero;
          MeasureType stdDev = NumericTraits<MeasureType>::Zero;
          HistogramFrequencyType totalMovingFrequency = NumericTraits<HistogramFrequencyType>::Zero;
           
          for (unsigned int m = 0; m < sm; m++)
            {
              index[0] = f;
              index[1] = m;
               
              movingFrequency = this->m_Histogram->GetFrequency(index);
              value = m+1;
              sum        += (movingFrequency * value);
              sumSquared += (movingFrequency * (value*value));
              totalMovingFrequency += movingFrequency;
               
            }
           
          if (totalMovingFrequency > 0)
            {
              mean = sum/(double)totalMovingFrequency;
              stdDev = sqrt( (sumSquared - ((sum*sum)/(double)totalMovingFrequency))/((double)totalMovingFrequency - 1));
              piu += ((stdDev/mean) * ((double)totalMovingFrequency/(double)totalFrequency));
           
            }
        }
      return piu;                                           
    }  

private:
  PIUImageToImageMetric(const Self&); // purposefully not implemented
  void operator=(const Self&);        // purposefully not implemented  
};

} // end namespace itk

#endif