cmtkJointHistogram.h

Go to the documentation of this file.

/*
//
//  Copyright 1997-2009 Torsten Rohlfing
//
//  Copyright 2004-2011 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: 2752 $
//
//  $LastChangedDate: 2011-01-17 11:33:31 -0800 (Mon, 17 Jan 2011) $
//
//  $LastChangedBy: torstenrohlfing $
//
*/

#ifndef __cmtkJointHistogram_h_included_
#define __cmtkJointHistogram_h_included_

#include <cmtkconfig.h>

#include <Base/cmtkJointHistogramBase.h>
#include <Base/cmtkUniformVolume.h>
#include <Base/cmtkTypedArray.h>
#include <Base/cmtkMathUtil.h>
#include <Base/cmtkTypes.h>
#include <Base/cmtkHistogram.h>

#include <System/cmtkSmartPtr.h>

#include <vector>
#include <algorithm>

namespace
cmtk
{


template<class T>
class JointHistogram : 
  public JointHistogramBase 
{
protected:
  size_t NumBinsX;

  Types::DataItem BinWidthX;

  Types::DataItem BinOffsetX;

  size_t NumBinsY;

  Types::DataItem BinWidthY;

  Types::DataItem BinOffsetY;

  std::vector<T> JointBins;

  size_t m_TotalNumberOfBins;
  
public:
  typedef JointHistogram<T> Self;

  typedef SmartPointer<Self> SmartPtr;

  JointHistogram () 
  {
    this->m_TotalNumberOfBins = NumBinsX = NumBinsY = 0; 
    BinWidthX = BinWidthY = 1.0;
    BinOffsetX = BinOffsetY = 0.0;
  }

  JointHistogram( const size_t numBinsX, const size_t numBinsY, const bool reset = true ) 
  {
    this->m_TotalNumberOfBins = NumBinsX = NumBinsY = 0; 
    BinWidthX = BinWidthY = 1.0;
    BinOffsetX = BinOffsetY = 0.0;
    this->Resize( numBinsX, numBinsY, reset );
  }

  Self* Clone () const 
  {
    return new Self( *this );
  }

  void Resize( const size_t numberOfBinsX, const size_t numberOfBinsY, const bool reset = true )
  {
    this->NumBinsX = numberOfBinsX;
    this->NumBinsY = numberOfBinsY;
    this->m_TotalNumberOfBins = this->NumBinsX * this->NumBinsY;

    this->JointBins.resize( this->m_TotalNumberOfBins );
    
    if ( reset ) 
      this->Reset();
  }

  size_t GetNumBinsX() const 
  { 
    return this->NumBinsX;
  }

  size_t GetNumBinsY() const
  { 
    return this->NumBinsY;
  }

  void SetRangeX( const Types::DataItemRange& range ) 
  {
    this->BinOffsetX = range.m_LowerBound;
    this->BinWidthX = range.Width() / ( this->NumBinsX - 1 );
  }

  void SetRangeY( const Types::DataItemRange& range ) 
  {
    this->BinOffsetY = range.m_LowerBound;
    this->BinWidthY = range.Width() / ( this->NumBinsY - 1 );
  }
  
  void SetRangeCenteredX( const Types::DataItemRange& range ) 
  {
    this->BinWidthX = range.Width() / (this->NumBinsX - 1);
    this->BinOffsetX = - this->BinWidthX / 2;
  }

  void SetRangeCenteredY( const Types::DataItemRange& range ) 
  {
    this->BinWidthY = range.Width() / (this->NumBinsY - 1);
    this->BinOffsetY = - this->BinWidthY / 2;
  }
  
  const Types::DataItemRange GetRangeX() const 
  {
    return Types::DataItemRange( this->BinOffsetX, this->BinOffsetX + this->BinWidthX * ( this->NumBinsX - 1) );
  }
  
  const Types::DataItemRange GetRangeY() const 
  {
    return Types::DataItemRange( this->BinOffsetY, this->BinOffsetY + this->BinWidthY * ( this->NumBinsY - 1) );
  }
  
  void Reset () 
  {
    std::fill( this->JointBins.begin(), this->JointBins.end(), 0 );
  }
  
  size_t ValueToBinX ( const Types::DataItem value ) const 
  {
    const size_t binIndex = static_cast<size_t>( (value - BinOffsetX) / BinWidthX );
    return std::max<size_t>( 0, std::min<size_t>( NumBinsX-1, binIndex ) );
  }
  
  size_t ValueToBinY ( const Types::DataItem value ) const 
  {
    const size_t binIndex = static_cast<int>( (value - BinOffsetY) / BinWidthY );
    return std::max<size_t>( 0, std::min<size_t>( NumBinsY-1, binIndex ) );
  }

  Types::DataItem BinToValueX ( const size_t bin ) const 
  {
    return BinOffsetX + (bin+0.5) * BinWidthX;
  }

  Types::DataItem BinToValueY ( const size_t bin ) const 
  {
    return BinOffsetY + (bin+0.5) * BinWidthY;
  }

  T ProjectToX ( const size_t indexX ) const 
  {
    T project = 0;
    for ( size_t j=0; j < NumBinsY; ++j )
      project += this->JointBins[indexX + j * NumBinsX];
    
    return project;
  }
  
  T ProjectToY ( const size_t indexY ) const 
  {
    T project = 0;
    size_t offset = indexY * NumBinsX;
    for ( size_t i = 0; i < NumBinsX; ++i )
      project += this->JointBins[i + offset];
    
    return project;
  }

  Histogram<T>* GetMarginalX() const;

  Histogram<T>* GetMarginalY() const;

  T GetBin( const size_t indexX, const size_t indexY ) const 
  {
    return this->JointBins[ indexX + NumBinsX * indexY ];
  }
  
  T SampleCount () const 
  {
    T sampleCount = 0;
    
    for ( size_t idx = 0; idx < this->m_TotalNumberOfBins; ++idx )
      {
      sampleCount += this->JointBins[idx];
      }
    
    return sampleCount;
  }

  T GetMaximumBinValue () const 
  {
    T maximum = 0;
    
    size_t idx = 0;
    for ( size_t i=0; i<NumBinsY; ++i ) 
      {
      for ( size_t j=0; j<NumBinsX; ++j, ++idx )
        maximum = std::max( maximum, this->JointBins[idx] );
      }
    
    return maximum;
  }
  
  void GetMarginalEntropies ( double& HX, double& HY ) const;

  double GetJointEntropy() const;
  
  void Increment ( const size_t sampleX, const size_t sampleY ) 
  {
    ++this->JointBins[sampleX+sampleY*NumBinsX];
  }

  void Increment 
  ( const size_t sampleX, const size_t sampleY, const T weight ) 
  {
    this->JointBins[ sampleX + sampleY * NumBinsX ] += weight;
  }

  void Decrement ( const size_t sampleX, const size_t sampleY ) 
  {
    --this->JointBins[sampleX+sampleY*NumBinsX];
  }
  
  void Decrement
  ( const size_t sampleX, const size_t sampleY, const Types::DataItem weight ) 
  {
    this->JointBins[ sampleX + sampleY * NumBinsX ] -= static_cast<T>( weight );
  }
  
  void AddJointHistogram ( const Self& other ) 
  {
    for ( size_t idx = 0; idx < this->m_TotalNumberOfBins; ++idx )
      this->JointBins[idx] += other.JointBins[idx];
  }
  
  void AddHistogramRow( const Histogram<T>& other, const size_t sampleY, const float weight = 1 ) 
  {
    size_t idx = this->NumBinsX * sampleY;
    for ( size_t i = 0; i<NumBinsX; ++i, ++idx )
      {
      this->JointBins[idx] += static_cast<T>( weight * other[i] );
      }
  }
  
  void AddHistogramColumn( const size_t sampleX, const Histogram<T>& other, const float weight = 1 ) 
  {
    size_t idx = sampleX;
    for ( size_t j = 0; j<NumBinsY; ++j, idx += NumBinsX )
      this->JointBins[idx] += static_cast<T>( weight * other[j] );
  }
  
  void RemoveJointHistogram ( const Self& other ) 
  {
    for ( size_t idx = 0; idx < this->m_TotalNumberOfBins; ++idx )
      {
      this->JointBins[idx] -= other.JointBins[idx];
      }
  }
  
  void NormalizeOverX( const double normalizeTo = 1.0 ) 
  {
    for ( size_t j = 0; j < NumBinsY; ++j ) 
      {
      const T project = this->ProjectToY( j );
      if ( project > 0 )
        {
        const double factor = normalizeTo / project;
        for ( size_t i = 0; i < NumBinsX; ++i )
          this->JointBins[ i + NumBinsX * j ] = static_cast<T>( this->JointBins[ i + NumBinsX * j ] * factor );
        }
      }
  }
  
  void NormalizeOverY( const double normalizeTo = 1.0 ) 
  {
    for ( size_t i = 0; i < NumBinsX; ++i ) 
      {
      const T project = this->ProjectToX( i );
      if ( project > 0 ) 
        {
        const double factor = normalizeTo / project;
        for ( size_t j = 0; j < NumBinsY; ++j )
          this->JointBins[ i + NumBinsX * j ] = static_cast<T>( this->JointBins[ i + NumBinsX * j ] * factor );
        }
      }
  }
  
  /* Return the index of the bin with the maximum value for one row.
   *\param j Index of the row.
   *\return The index of the bin with the maximum value in row j.
   */
  size_t GetMaximumBinIndexOverX( const size_t j ) const 
  {
    size_t offset = j * NumBinsX;
    
    size_t maxIndex = 0;
    T maxValue = this->JointBins[ offset ];
    
    for ( size_t i = 1; i < NumBinsX; ++i ) 
      {
      offset++;
      if ( this->JointBins[ offset ] > maxValue ) 
        {
        maxValue = this->JointBins[ offset ];
        maxIndex = i;
        }
      }
    return maxIndex;
  }
  
  /* Return the index of the bin with the maximum value for one column.
   *\param j Index of the column.
   *\return The index of the bin with the maximum value in column j.
   */
  size_t GetMaximumBinIndexOverY( const size_t i ) const 
  {
    size_t offset = i;
    
    size_t maxIndex = 0;
    T maxValue = this->JointBins[ offset ];
    
    for ( size_t j = 1; j < NumBinsY; ++j ) 
      {
      offset += NumBinsX;
      if ( this->JointBins[ offset ] > maxValue ) 
        {
        maxValue = this->JointBins[ offset ];
        maxIndex = j;
        }
      }
    return maxIndex;
  }

  double GetMutualInformation( const bool normalized = false ) const 
  {
    double hX, hY, hXY;
    this->GetMarginalEntropies( hX, hY );
    hXY = this->GetJointEntropy();
    if ( hXY > 0 )
      if ( normalized ) 
        return (hX + hY) / hXY;
      else
        return (hX + hY) - hXY;
    else
      return 0;
  }

  double GetCorrelationRatio() const 
  {
    T sampleCount = this->SampleCount();
    if ( ! sampleCount ) return 0;
    
    double sigSquare = 0, m = 0;
    for ( size_t j = 0; j < NumBinsY; ++j ) 
      {
      sigSquare += ( MathUtil::Square(j) * ProjectToY(j) );
      m += (j * ProjectToY(j));
      }
    
    double invSampleCount = 1.0 / sampleCount;
    
    m *= invSampleCount;
    (sigSquare *= invSampleCount) -= MathUtil::Square(m);
    
    double eta = 0;
    for ( size_t i = 0; i < NumBinsX; ++i ) 
      {
      if ( ProjectToX(i) > 0 ) 
        {
        double sigSquare_i = 0, m_i = 0;
        for ( size_t j = 0; j < NumBinsY; ++j ) 
          {
          sigSquare_i += (MathUtil::Square(j) * this->JointBins[ i + j * NumBinsX ]);
          m_i += (j * this->JointBins[ i + j * NumBinsX ]);
          }
        m_i *= (1.0 / ProjectToX(i) );
        (sigSquare_i *= (1.0 / ProjectToX(i))) -= MathUtil::Square(m_i);
        eta += (sigSquare_i * ProjectToX(i));
        }
      }
    
    return eta / (sigSquare * sampleCount); 
  }
};


} // namespace cmtk

#endif // #ifndef __cmtkJointHistogram_h_included_