cmtkGaussianKernel.h

Go to the documentation of this file.

/*
//
//  Copyright 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: 2672 $
//
//  $LastChangedDate: 2010-12-15 14:40:41 -0800 (Wed, 15 Dec 2010) $
//
//  $LastChangedBy: torstenrohlfing $
//
*/

#ifndef __cmtkGaussianKernel_h_included_
#define __cmtkGaussianKernel_h_included_

#include <cmtkconfig.h>

#include <Base/cmtkUnits.h>
#include <Base/cmtkMathUtil.h>

#include <vector>


namespace
cmtk
{

template<class TFloat=double>
class GaussianKernel
{
public:
  typedef GaussianKernel<TFloat> Self;

  static std::vector<TFloat> GetSymmetricKernel( const Units::GaussianSigma& sigma , 
                                                 const TFloat maxError = 1e-5  )
  {
    const double normFactor = 1.0/(sqrt(2*M_PI) * sigma.Value());
    const size_t radius = Self::GetRadius( sigma, normFactor, maxError );
    
    std::vector<TFloat> kernel( 2 * radius + 1 );
    for ( size_t i = 0; i <= radius; ++i )
      {
      kernel[radius-i] = kernel[radius+i] = normFactor * exp( -MathUtil::Square( 1.0 * i / sigma.Value() ) / 2 );
      }

    return kernel;
  }

  static std::vector<TFloat> GetHalfKernel( const Units::GaussianSigma& sigma , 
                                            const TFloat maxError = 1e-5  )
  {
    const double normFactor = 1.0/(sqrt(2*M_PI) * sigma.Value());
    const size_t radius = static_cast<size_t>( Self::GetRadius( sigma, normFactor, maxError ) );
    
    std::vector<TFloat> kernel( radius + 1 );
    for ( size_t i = 0; i <= radius; ++i )
      {
      kernel[i] = normFactor * exp( -MathUtil::Square( 1.0 * i / sigma.Value() ) / 2 );
      }
    
    return kernel;
  }

private:
  static TFloat GetRadius(  const Units::GaussianSigma& sigma, const TFloat normFactor, const TFloat maxError )
  {
    if ( maxError >= normFactor ) // if normFactor is less than max error, then we really need no kernel at all
      return 0;
    else
      return static_cast<size_t>( sqrt( -2.0 * log( maxError / normFactor ) ) * sigma.Value() );
  }
};

} // namespace cmtk


#endif // #ifndef __cmtkGaussianKernel_h_included_