cmtkVoxelMatchingCorrRatio.cxx

Go to the documentation of this file.

/*
//
//  Copyright 1997-2009 Torsten Rohlfing
//
//  Copyright 2004-2010 SRI International
//
//  This file is part of the Computational Morphometry Toolkit.
//
//  http://www.nitrc.org/projects/cmtk/
//
//  The Computational Morphometry Toolkit is free software: you can
//  redistribute it and/or modify it under the terms of the GNU General Public
//  License as published by the Free Software Foundation, either version 3 of
//  the License, or (at your option) any later version.
//
//  The Computational Morphometry Toolkit is distributed in the hope that it
//  will be useful, but WITHOUT ANY WARRANTY; without even the implied
//  warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License along
//  with the Computational Morphometry Toolkit.  If not, see
//  <http://www.gnu.org/licenses/>.
//
//  $Revision: 2022 $
//
//  $LastChangedDate: 2010-07-21 15:26:03 -0700 (Wed, 21 Jul 2010) $
//
//  $LastChangedBy: torstenrohlfing $
//
*/

#include "cmtkVoxelMatchingCorrRatio.h"

namespace
cmtk
{


template<Interpolators::InterpolationEnum I>
typename VoxelMatchingCorrRatio<I>::ReturnType
VoxelMatchingCorrRatio<I>::Get () const
{
  double invSampleCount = 1.0 / HistogramI.SampleCount();
  // initialize variable for the weighted sum of the sigma^2 values over all
  // reference intensity classes.
  double sumSigmaSquare = 0;
  // run over all bins, i.e., reference classes
  for ( unsigned int j = 0; j < NumBinsX; ++j ) 
    {
    // are there any values in the current class?
    if ( HistogramI[j] ) 
      {
      // compute mean floating value for this reference class
      double mu = SumJ[j] / HistogramI[j];
      // compute variance of floating values for this reference class
      double sigmaSq = ( mu*mu*HistogramI[j] - 2.0*mu*SumJ[j] + SumJ2[j] ) / HistogramI[j]; 
      // update sum over all classes with weighted sigma^2 for this class.
      sumSigmaSquare += (invSampleCount * HistogramI[j]) * sigmaSq;
      }
    }
  
  // compute (supposedly) correlation ratio
  typename Self::ReturnType cr = static_cast<typename Self::ReturnType>( 1.0 - (1.0 /  SigmaSqJ ) * sumSigmaSquare );
  
  sumSigmaSquare = 0;
  for ( unsigned int i = 0; i < NumBinsY; ++i ) 
    {
    if ( HistogramJ[i] ) 
      {
      double mu = SumI[i] / HistogramJ[i];
      double sigmaSq = ( mu*mu*HistogramJ[i] - 2.0*mu*SumI[i] + SumI2[i] ) / HistogramJ[i]; 
      // update sum over all classes with weighted sigma^2 for this class.
      sumSigmaSquare += (invSampleCount * HistogramJ[i]) * sigmaSq;
      }
    }
  
  // add reverse correlation ratio
  cr += static_cast<typename Self::ReturnType>(1.0 - (1.0 /  SigmaSqI ) * sumSigmaSquare);
  
  return cr;
}

template class VoxelMatchingCorrRatio<Interpolators::LINEAR>;

template class VoxelMatchingCorrRatio<Interpolators::NEAREST_NEIGHBOR>;

} // namespace cmtk