---

kernelfactory.c

Go to the documentation of this file.

#include "kernelfactory.h"
#include "cmdparam.h"


/* Static variables in KernelFactory */
double KernelFactory::blur_radius=0;
double KernelFactory::blur_radius_surround_multiplier=1.6;
double KernelFactory::blur_scale=1.0;
double KernelFactory::blur_range_multiplier=1.0;
int    KernelFactory::blur_type=Blur_CircularAverage;
char   KernelFactory::default_kernel_filename[CMD_MAX_LINE_LENGTH] = "kernel.pgm";



void  KernelFactory::register_params_and_commands( void )
{
  PARAM_II(PARAM_INT,   blur_type,&blur_type,0,int(KernelFactory::max_blur_type),
          "Type of blurring, smoothing, or (in general) convolution to use on the\n"
          "inputs to be presented, when layers_pereye>1.  Available types include\n"
          "the following (prefixed with 'Blur_'):\n"
          "  CircularAverage,SquareAverage,Gaussian,DoG,LoG,GoD,DoGGoD");
  params_define_default_expr("blur_type","Blur_CircularAverage");
  
  PARAM_LL(PARAM_DOUBLE,blur_radius,&blur_radius,0,
          "Default blurring or smoothing radius for the input, in retinal receptors.\n"
          "If zero, the identity function is applied.  This parameter determines the\n"
          "nominal radius, but the various blur_type options each produce a function\n"
          "with a different actual width.  For instance, a Gaussian with sigma blur_radius\n"
          "is much narrower than a uniform circle of radius blur_radius, yet wider\n"
          "than the Laplacian of a Gaussian.  See also blur_type and\n"
          "blur_range_multiplier.\n\n"

          "Compatibility note: in previous versions, setting blur_radius to zero\n"
          "disabled retinal ganglion cell processing, while setting it to a nonzero\n"
          "value enabled it.  Such behavior is now controlled by the layers_pereye\n"
          "parameter instead, and no longer depends on blur_radius.");
  
  PARAM_LL(PARAM_DOUBLE,blur_radius_surround_multiplier,&blur_radius_surround_multiplier,0,
          "Ratio of surround radius to blur_radius.  Defaults to 1.6, which makes the\n"
          "DoG approximate the Laplacian derivative of a Gaussian (see Blur_LoG,\n"
          "\\cite{marr:prslb80,marr:vision}), albeit at a different spatial scale.");
  
  PARAM_LL(PARAM_DOUBLE,blur_range_multiplier,&blur_range_multiplier,0,
          "Each non-constant blur_type (i.e., those other than Blur_SquareAverage\n"
          "and Blur_CircularAverage) has a default maximum circular radius at which\n"
          "to evaluate the blurring function (i.e., clipping).  The default should\n"
          "usually be sufficient, but it can be increased for greater accuracy (or\n"
          "decreased for speed) by making this parameter higher or lower than 1.0,\n"
          "respectively.  Since clipping is done circularly, the shape of an isotropic\n"
          "kernel should remain approximately isotropic even when the multiplier is set\n"
          "far enough below 1.0 that heavy clipping occurs.");
  
  PARAM_LL(PARAM_DOUBLE,blur_scale,&blur_scale,0,
          "Intensity scale for the blurring kernel, if any.  This is primarily useful\n"
          "for blurring types which change the input strength dramatically, such as\n"
          "Blur_DoG or a Blur_Gaussian with a large radius.  This scaling factor is\n"
          "used only when layers_pereye>1, so it can help make the input strength\n"
          "independent of whether blurring is in use.  At present it is the\n"
          "equivalent of gammaaff for a FixedWtRegion, and should probably\n"
          "be phased out in favor of gammaaff for simplicity.");
  
  CONST_I(Blur_CircularAverage,int(Blur_CircularAverage),False,
          "Blur type (see command input_define_convolution) with parameters:\n"
          "  [<radius> [<scale>]]\n\n"
          
          "Replaces each input pixel with the given scale times the arithmetic mean\n"
          "of the pixels in the surrounding circular region with the specified radius.\n\n"
          
          "For a very small radius, the approximation to a circle is poor (e.g.,\n"
          "it is square when radius=1).  That property can make this otherwise\n"
          "isotropic filter anisotropic.  This type of blurring doesn't work\n"
          "particularly well with patterns that have sharp edges, since the\n"
          "blurring pattern also has a sharp edge.  However, this type of blurring\n"
          "is much faster than ones like Blur_Gaussian, since a smaller convolution\n"
          "kernel is needed, and a flat one can sometimes be simpler to reason about.");

  CONST_I(Blur_CircularRandom,int(Blur_CircularRandom),False,
          "Blur type (see command input_define_convolution) with parameters:\n"
          "  [<radius> [<scale>]]\n\n"
          
          "Replaces each input pixel with the convolution of a circular region (with\n"
          "the specified radius) of random values with the input image surrounding\n" 
          "that pixel.  The values are normalized to have a total sum of 1.0 by\n" 
          "default, but if a different scale is specified, the normalized sum\n"
          "is then multiplied by it.\n\n"

          "This type of blurring is a default against which some more effective\n"
          "variety may be compared.");
  
  CONST_I(Blur_SquareAverage,int(Blur_SquareAverage),False,
          "Blur type (see command input_define_convolution) with parameters:\n"
          "  [<radius> [<scale>]]\n\n"

          "Replaces each input pixel with the given scale times the arithmetic mean\n"
          "of the pixels in the surrounding square region with the specified integer\n"
          "radius.\n\n"

          "This type of blurring is not isotropic, i.e. it blurs differently\n"
          "for different directions, so it should usually be avoided.");
  
  CONST_I(Blur_SquareRandom,int(Blur_SquareRandom),False,
          "Blur type (see command input_define_convolution) with parameters:\n"
          "  [<radius> [<scale>]]\n\n"
          
          "Replaces each input pixel with the convolution of a square region (with\n"
          "the specified radius) of random values with the input image surrounding\n" 
          "that pixel.  The values are normalized to have a total sum of 1.0 by\n" 
          "default, but if a different scale is specified, the normalized sum\n"
          "is then multiplied by it.\n\n"

          "This type of blurring is a default against which some more effective\n"
          "variety may be compared.");
  
  CONST_I(Blur_Gaussian,int(Blur_Gaussian),False,
          "Blur type (see command input_define_convolution) with parameters:\n"
          "  [<radius> [<scale> [<radius_range_scale>]]]\n\n"

          "Replaces each input pixel with the convolution of a Gaussian (with the\n"
          "specified radius) with the input image surrounding that pixel.  The\n" 
          "Gaussian is normalized to have a total sum of 1.0 by default, but if a\n" 
          "different scale is specified, the normalized Gaussian is multiplied by\n"
          "it.\n\n"

          "The function is continuous but is only evaluated up to a finite distance from\n"
          "its center in all directions.  If desired, a scale for this maximum radius may\n"
          "be supplied, which is multiplied by the function's radius to get the actual\n"
          "range. By default, a range is used which ensures that the function is nearly\n"
          "zero without having a kernel that is too large, which would slow down the\n"
          "program.  The blur_range_multiplier can be used to increase or decrease this\n"
          "default automatically.\n\n"
          
          "This type of blurring is effective at removing detail at a spatial scale below\n"
          "the radius of the Gaussian, without introducing edge effects like\n"
          "Blur_CircularAverage.  (It is essentially a low-pass filter with a cutoff\n"
          "frequency determined by the blur_radius.)");
  
  CONST_I(Blur_LoG,int(Blur_LoG),False,
          "Blur type (see command input_define_convolution) with parameters:\n"
          "  [<radius> [<scale> [<radius_range_scale>]]]\n\n"

          "Replaces each input pixel with the convolution of the Laplacian derivative\n"
          "of a Gaussian (with radius blur_radius and height 1.0) with the input image\n"
          "surrounding that pixel.  The scale is currently arbitrary, but it can be\n" 
          "adjusted with the scale parameter.  See Blur_Gaussian for a description of\n"
          "the radius_range_scale parameter.\n\n"

          "This type of blurring is nearly identical in shape to a Blur_DoG, and\n"
          "\\cite{marr:vision} has proposed that a DoG is approximating this\n"
          "function.  However, the height and spatial scales differ substantially\n"
          "from a DoG, partially because the DoG can easily be normalized to\n"
          "sum to zero, so it may be difficult to use this blurring type (though\n"
          "no less computationally efficient).");
  
  
  CONST_I(Blur_DoG,int(Blur_DoG),False,
          "Blur type (see command input_define_convolution) with parameters:\n"
          "  [<radius> [<surround_radius> [<scale> [<surround_scale> [<radius_range_scale>]]]]\n\n"
          
          "Replaces each input pixel with the convolution of a Difference of Gaussians\n"
          "(a Mexican-hat whose center Gaussian sums to the given scale and has the given\n"
          "radius and whose surround (negative) Gaussian sums to the given surround_scale\n"
          "and has the given surround_radius) with the input image surrounding that pixel.\n"
          "The surround_radius defaults to be blur_radius*blur_radius_surround_multiplier.\n\n"
          
          "See Blur_Gaussian for a description of the radius_range_scale parameter; the\n"
          "surround radius is used as a base since it is presumably larger than the other.\n\n"
          
          "This type of blurring serves to detect edges and similar structures at a scale\n"
          "somewhat smaller than the radius of the center Gaussian, and resembles\n"
          "the processing done by a population of RF-size-matched retinal ganglion cells.\n"
          "(It is basically a band-pass spatial filter.)");
  
  CONST_I(Blur_GoD,int(Blur_GoD),False,
          "Blur type (see command input_define_convolution) with parameters:\n"
          "  [<radius> [<surround_radius> [<scale> [<surround_scale> [<radius_range_scale>]]]]\n\n"
          
          "This blur_type is the same as Blur_DoG except that the surround Gaussian\n"
          "is taken as positive while the center is negative (OFF-center).  Given that\n"
          "the retina itself is always positive, this filter serves to highlight changes\n"
          "of the opposite sign as Blur_DoG (e.g. a dark area surrounded by light).");
  
  CONST_I(Blur_DoGGoD,int(Blur_DoGGoD),False,
          "Blur type (see command input_define_convolution) with parameters:\n"
          "  [<radius> [<surround_radius> [<scale> [<surround_scale> [<radius_range_scale>]]]]\n\n"
          
          "When this blur_type is used, even-numbered eyes (starting at zero) use\n"
          "Blur_DoG while the rest use Blur_GoD.  Assuming every pair of eyes\n"
          "has identical input patterns, the pair thus approximates a full set of\n"
          "ON-center and OFF-center inputs.  Presumably, num_eyes should be 2 or 4\n"
          "for this blurring type to be meaningful.");
  
  CONST_I(Blur_PGM,int(Blur_PGM),False,
          "Blur type (see command input_define_convolution) with parameters:\n"
          "  [<filename> [<offset> [<scale> [<normalize>]]]]\n\n"
          
          "The given PGM image file (\"kernel.pgm\" by default) is read and\n"
          "used as a convolution kernel.  The value range specified in the file\n"
          "is scaled to a range from 0 to 1.0, then the offset (-0.5 by default)\n"
          "is added. Next, if normalization is enabled (on by default), the\n"
          "matrix is normalized to have a sum of the given scale; otherwise\n"
          "the matrix is simply multiplied by the given scale.");
}



KernelFactory::KernelMatrix KernelFactory::create( int index, double wss, StringArgs arglist )
{
  KernelMatrix kernel;
  const int type = arglist.next(blur_type);
  string errors;
  
  /* Types are grouped into classes based on arguments accepted */
  switch (type) {
    
  case Blur_PGM: {
    const string filename  = arglist.next(string(default_kernel_filename));
    const double offset    = arglist.next(-0.5);
    const double scale     = arglist.next(blur_scale);
    const int    normalize = arglist.next(1);
    KernelMatrix new_kernel;
    
    if (!pgm_read(filename, new_kernel)) {
      ipc_notify(IPC_ONE,IPC_ERROR, "Can't read file: %s", filename.c_str());
      break;
    }
    
    if (new_kernel.ncols()%2 != 1 || new_kernel.nrows()%2 != 1){
      ipc_notify(IPC_ONE,IPC_ERROR,
                 "The convolution matrix must have an odd number of columns and rows (%d,%d)",
                 new_kernel.ncols(),new_kernel.nrows());
      break;
    }
    new_kernel += offset;
    
    if (!normalize)
      new_kernel *= scale;
    else {
      const double kernelsum = mat::sum(new_kernel);
      if (kernelsum==0)
        ipc_notify(IPC_ONE,IPC_ERROR,"Can't normalize blurring kernel since it sums to zero");
      else new_kernel *= fabs(scale/kernelsum);
    }
    
    kernel = new_kernel;
    break;
  }
  
  
  case Blur_SquareAverage: case Blur_CircularAverage: case Blur_Gaussian: case Blur_LoG: {
    
    const double rad   = arglist.next(blur_radius*wss);
    const double scale = arglist.next(blur_scale);
    
    switch (type) {
    case Blur_SquareAverage:   kernel = create(RadialFunction::Constant, rad,rad-0.5,errors,false); break;
    case Blur_CircularAverage: kernel = create(RadialFunction::Constant, rad,rad-0.5,errors); break;
    case Blur_Gaussian: {
      const double rscale = arglist.next(4.7*blur_range_multiplier);
      kernel = create(RadialFunction::Gaussian, rad, rad*rscale,errors); break;
    }
    case Blur_LoG: {
      const double rscale = arglist.next(3.2*blur_range_multiplier);
      kernel = create(RadialFunction::LoG, rad, rad*rscale,errors,true,false); break;
    }
    }
    kernel *= scale;
    break;
  }
  
  case Blur_DoG: case Blur_GoD: case Blur_DoGGoD: {
    /* Determine sign of kernel; DoGGoD uses 1 for even, -1 for odd */
    const double scalesign   = (  type == Blur_DoG    ?  1.0           :
                                ( type == Blur_GoD    ? -1.0            :
                                 (type == Blur_DoGGoD ?  1.0-2*(index%2) :
                                  1.0)));
    
    const double rad         = arglist.next((blur_radius)*wss);
    
    const double srad        = arglist.next((blur_radius_surround_multiplier*rad/wss)*wss);
    const double scale       = arglist.next(blur_scale);
    const double sscale      = arglist.next(scale);
    const double range_scale = arglist.next(4.7*blur_range_multiplier);
    const double maxrad      = srad*range_scale;
    
#ifdef USE_MTL_MATRIX   
    kernel  = create(RadialFunction::Gaussian,  rad, maxrad,errors);
    mtl::scale(kernel,scale*scalesign);
    mtl::add( mtl::scaled(create(RadialFunction::Gaussian, srad, maxrad,errors), -1*sscale*scalesign), kernel);
#else
    kernel = create(RadialFunction::Gaussian, rad,  maxrad,errors) * scale*scalesign  -
             create(RadialFunction::Gaussian, srad, maxrad,errors) * sscale*scalesign;
#endif
    break;
  }

  case Blur_CircularRandom: case Blur_SquareRandom: {
    const double rad   = arglist.next(blur_radius*wss);
    const double scale = arglist.next(blur_scale);
    
    switch (type) {
    case Blur_SquareRandom:   kernel = create(RadialFunction::Random, rad,rad-0.5,errors,false); break;
    case Blur_CircularRandom: kernel = create(RadialFunction::Random, rad,rad-0.5,errors); break;
    }
    kernel *= scale;
    break;
  }
    
  default: ipc_notify(IPC_ONE,IPC_ERROR,"Unknown blur_type %d",type); break;
  }

  if (errors!="")
    ipc_notify(IPC_ONE,IPC_ERROR,errors.c_str());

  return kernel;
}

---