---

radialfunction.h

Go to the documentation of this file.

#ifndef __RADIAL_FUNCTION_H__
#define __RADIAL_FUNCTION_H__

#include <cmath>
#include <string>
using std::string;

#include "shuffledrand.h"
#include "mathconstants.h"


namespace RadialFunction {

  template<class T>
  inline T Constant( T, T )
    {  return 1.0;  }

  template<class T>
  inline T Random( T, T )
  {  return shuffled_rand();  }

  template<class T>
  inline T Gaussian( T x_sq, T sig_sq)
    {  return exp( -x_sq/(2*(Math::pi)*sig_sq) );  }

  template<class T>
  inline T LoG( T x_sq, T sig_sq) { 
    const T exponent = -x_sq/(2*sig_sq);
    return 1/((Math::pi)*sig_sq*sig_sq)*(1.0+exponent) * exp(exponent);
  }

  
  template<class matrix_type, class radial_function, class radius_type>
    matrix_type matrix
    (/* Function pointer or function object of the form T fn(T,T) */
     radial_function radial_fn,
     radius_type radius,
     radius_type max_radius_,
     string& errors,
     bool circular=true,
     bool normalize=true
     )
  {
    typedef typename matrix_type::value_type value_type;
    typedef typename matrix_type::size_type  size_type;
    
    const radius_type radius_sq      = radius*radius;
    /* Uses floating-point max_radius for circular radius, integer for other tasks */
    const radius_type max_radius_sq  = max_radius_*max_radius_;
    const size_type   max_radius     = size_type(floor((0.5+(max_radius_))));
    const size_type   kernelwidth    = std::max(size_type(1),size_type(max_radius*2+1));
    matrix_type kernel(kernelwidth,kernelwidth);
    
    /* Cropping is counterproductive for very small radii */
    if (max_radius<=1) circular=false;
    
    /* Create kernel matrix */
    value_type kernelsum=0;
    /*   ki and kj cannot be size_type since that might be unsigned */
    int ki = -max_radius;
    for (size_type i=0; i < kernelwidth; ki++,i++) {
      int kj = -max_radius;
      for (size_type j=0; j < kernelwidth; kj++,j++) {
        const size_type ki_sq = ki*ki;
        const size_type kj_sq = kj*kj;
        const size_type r_sq  = ki_sq + kj_sq;
        
        if (circular && (r_sq > max_radius_sq))
          kernel[i][j] = 0;
        else {
          const value_type val = (radial_fn)(r_sq,radius_sq);
          kernelsum += val;
          kernel[i][j] = val;
        }
      }
    }
    //cout << "Unnormalized sum of blurring kernel values is " <<kernelsum;
    
    /* Normalize kernel matrix */
    if (normalize) {
      if (kernelsum==0)
        errors += "Can't normalize a kernel that sums to zero";
      else
        kernel *= 1/kernelsum;
    }
    
    return kernel;
  }

} /* namespace */


#endif /* __RADIAL_FUNCTION_H__ */

---