---

retinalobjs.c

Go to the documentation of this file.

#include <cmath> 
#include <iomanip>
#include <strstream.h>
#include <numeric>
#include <fstream>
#include <utility>

#include "retinalobjs.h"

#ifndef NO_VALGEN_STRINGS
#include <stdio.h>
#include "stringutils.h"
#endif

#ifdef ANSI_COMPATIBLE
/* Optional: provides definition of hypot on systems lacking it */
#include "ind_types.h"
#endif



/******************************************************************************/
/*  Retinal_Object                                                            */
/******************************************************************************/

//#include <typeinfo> /* Uncomment if typeid is enabled below */

string Retinal_Object::stringrep() const
{
  /* Should be modified to use an ostringstream when sstream becomes available;
     also use the "fixed" iomanip instead of the setf below. */
  const char bufsize=100;
  static char buf[bufsize];
  ostrstream os(buf, bufsize);
  os.setf(std::ios::fixed,std::ios::floatfield);
  //os << typeid(*this).name() << " "; /* For debugging */
  if (name()!="") os << name() << " ";
  if (get_active()) {
    os << std::setprecision(1) 
       <<  "cx:"    << std::setw(4) << std::setfill('0') << get_var("cx")
       << " cy:"    << std::setw(4) << std::setfill('0') << get_var("cy")
       << " theta:" << std::setw(5) << std::setfill('0')
       << RADIANS_TO_DEGREES(CONSTRAIN_ANGLE(     get_var("theta")));
  }
  else
    os << "inactive";
      
  /* For debugging or other reasons, print all variables */
  //for(VarMap::const_iterator i=vars.begin(); i!= vars.end(); i++) {
  //  const string& nam= i->first;
  //  if (name != "cx" && nam != "cy" && nam != "theta")
  //      os << " " << nam << ":" << i->second->value();
  //}

  os << std::ends;
  
  return buf;
}


#ifndef NO_VALGEN_STRINGS
Retinal_Object::VarMap RetinalObjectStringArgs::vars(const ParamList& params)
{
  Retinal_Object::VarMap vars;
  
  /* Parse all positional arguments possible, substituting default if arg is missing */
  for (ParamList::const_iterator paramptr=params.begin();
       paramptr!=params.end(); paramptr++) {
    const string& name = *paramptr;
    const string  nextarg = args.top(string(""));
    const unsigned pos = nextarg.find("=");
    const bool    nextarghaseq = (pos<nextarg.length());
    const bool    positionalargsremain = !args.empty() && !nextarghaseq;
    const string  value = (positionalargsremain ? args.next(string("")) : get_default(name));
    set_linked(name,value,vars);
  }
  
  /* Parse any named arguments present, overriding defaults and positional arguments */
  while (!args.empty()) {
    const string   arg   = args.next(string(""));
    const unsigned pivot = arg.find("=");
    if (pivot>=arg.length()) {
      error("Syntax error in retinal object parameter specification");
      return vars;
    }
    const string name  = string(arg,0,pivot);
    const string value = String::strip_quotes(string(arg,pivot+1,arg.length()-pivot-1));
    set_linked(name,value,vars);
  }

  return vars;
}
#endif /* NO_VALGEN_STRINGS */



/******************************************************************************/
/*  Retinal_Composite objects                                                 */
/******************************************************************************/


bool Retinal_Composite::next()
{
  return
    !std::count_if(ISEQ(children),std::not1(std::mem_fun(&ValueGen::next)))
    && Retinal_Object::next();
}



void Retinal_Composite::reset() {
  dominant_child = (accum_type==OneHot && children.size() ? children[0] : 0);
  std::for_each(ISEQ(children),std::mem_fun(&Retinal_Object::reset));
  Retinal_Object::reset();
}



#ifndef NO_VALGEN_STRINGS
RetinalObjectStringArgs::ParamList Retinal_Composite::paramlist() {
  RetinalObjectStringArgs::ParamList p;
  p.push_back("theta");
  p.push_back("cx");
  p.push_back("cy");
  p.push_back("size_scale");
  p.push_back("scale");
  p.push_back("offset");
  p.push_back("accum_type");
  p.push_back("hot");
  return p;
}
#endif



bool Retinal_Composite::update() const
{
  std::for_each(ISEQ(children),std::mem_fun(&Retinal_Object::update));

  /* Select single dominant child, if any */
  if (accum_type==OneHot && children.size()) {
    const Variable hot = get_var("hot");
    const Variable relative_index = hot-int(hot); /* Value between 0 and 1 */
    const int      child_index    = (int)( children.size()*relative_index );
    assert(child_index<(int)children.size());
    dominant_child = children[child_index];
  }
  else
    dominant_child=0;
  
  /* Cache the result of sub-computations */
  cx                        = get_var("cx");
  cy                        = get_var("cy");
  const Variable mt         = -get_var("theta");
  cosmt                     = cos(mt);
  sinmt                     = sin(mt);
  const Variable size_scale = get_var("size_scale");
  div_size                  = 1/size_scale;


  /* Construct bounding box enclosing all children in local coordinates  */
  bounding_box.set(0,0,0,0); /* Initially empty */
  for(const_iterator i=children.begin(); i!=children.end(); i++)
    bounding_box += (*i)->bounding_box;
  
  /* Convert the bounding box into world coordinates */
  bounding_box.scale(size_scale,size_scale).rotate(mt).translate(cx,cy);

  return Retinal_Object::update();
}


class ActivityAccumulator : public std::binary_function<Retinal_Object*,Retinal_Object*,bool> {
public:
  explicit ActivityAccumulator(const Retinal_Obj::Coordinate xi, const Retinal_Obj::Coordinate yi,
                               Retinal_Composite::AccumulationType accum_type_i)
    : x(xi), y(yi), accum_type(accum_type_i) { }

  Retinal_Obj::Activity operator() (Retinal_Obj::Activity prev, const Retinal_Object* r) {
    const Retinal_Obj::Activity now=r->activation(x,y);
    switch (accum_type) {
    case Retinal_Composite::Max: return std::max(prev,now);
    case Retinal_Composite::Min: return std::min(prev,now);
    case Retinal_Composite::Add: return prev+now;
    case Retinal_Composite::Replace:
      return (r->inside(x,y)? now : prev);
    default: return prev+now;
    }
  }
  
private:
  Retinal_Obj::Coordinate x,y;
  const Retinal_Composite::AccumulationType accum_type;
};



class DefaultActivityAccumulator : public std::binary_function<Retinal_Object*,Retinal_Object*,bool> {
public:
  explicit DefaultActivityAccumulator(Retinal_Composite::AccumulationType accum_type_i)
    : accum_type(accum_type_i) { }

  Retinal_Obj::Activity operator() (Retinal_Obj::Activity prev, const Retinal_Object* r) {
    const Retinal_Obj::Activity now=r->default_activation();
    switch (accum_type) {
    case Retinal_Composite::Max: return std::max(prev,now);
    case Retinal_Composite::Min: return std::min(prev,now);
    case Retinal_Composite::Add: return prev+now;
    case Retinal_Composite::Replace: return prev;
    default: return prev+now;
    }
  }
  
private:
  const Retinal_Composite::AccumulationType accum_type;
};



Retinal_Obj::Activity Retinal_Composite::accum_base() const {
  switch (accum_type) {
    /* Returns activation value assumed to be higher than all others */
  case Min: return  10000000;
    
    /* Returns activation value assumed to be lower  than all others */
  case Max: return -10000000; 

    /* Returns the default activation value of the first child (if any),
       assumed to be a background pattern. */
  case Replace: return children.size()? children[0]->default_activation() : 0;

  default:  return 0;
  }
}



Retinal_Object::Activity Retinal_Composite::default_activ() const
{
  //return =0.5; /* To debug Composite bounding boxes */
  return 
    (children.empty() ? 0 
     : (dominant_child? dominant_child->default_activation()
        : std::accumulate(ISEQ(children),accum_base(),DefaultActivityAccumulator(accum_type))));
}


  
Retinal_Object::Activity Retinal_Composite::activ(Coordinate x, Coordinate y) const
{
  const Coordinate dx      = (x-cx)*div_size;
  const Coordinate dy      = (y-cy)*div_size;
  const Coordinate xp      = dx*cosmt-dy*sinmt;
  const Coordinate yp      = dx*sinmt+dy*cosmt;
  
  return (dominant_child ? dominant_child->activation(xp,yp)
          : std::accumulate(ISEQ(children),accum_base(),ActivityAccumulator(xp,yp,accum_type)));
}



const Retinal_Object& Retinal_Composite::mostactive(Coordinate x, Coordinate y) const
{
  if (children.begin()==children.end())
    return *this;
  else {
    Retinal_Object* champ = *children.begin();
    Activity high_act = champ->activation(x,y);
    Activity act;
    for (Retinal_Composite::const_iterator i = children.begin()+1; i!= children.end(); i++) {
      act = (*i)->activation(x,y);
      if (act>high_act) {
        champ = *i;
        high_act=act;
      }
    }
    return *champ;
  }
}



string accumulate_stringreps(string val, Retinal_Object* r)
{  return val + " [" + r->stringrep() + "]";  }



string Retinal_Composite::stringrep() const
{
  return Retinal_Object::stringrep() +
    (dominant_child ? " " + dominant_child->stringrep()
     : std::accumulate(ISEQ(children),string(""),accumulate_stringreps));
}



bool Retinal_Composite::inside(Coordinate x, Coordinate y) const
{
  /* Normal case */
  if (!(accum_type==OneHot && dominant_child))
    return bounding_box.inside(x,y);

  /* When there is one dominant child, ask it instead */
  const Coordinate dx      = (x-cx)*div_size;
  const Coordinate dy      = (y-cy)*div_size;
  const Coordinate xp      = dx*cosmt-dy*sinmt;
  const Coordinate yp      = dx*sinmt+dy*cosmt;
  
  return dominant_child->inside(xp,yp);
}


/******************************************************************************/
/*  Retinal_ManagedComposite objects                                          */
/******************************************************************************/


bool Retinal_ManagedComposite::distance_valid (const Retinal_Object& obj1, const Retinal_Object& obj2)
{
  if (!*min_dist_enforce && !*max_dist_enforce)
    return true;

  const Coordinate dist = Generic::hypot(obj1.get_var("cx") - obj2.get_var("cx"),
                                         obj1.get_var("cy") - obj2.get_var("cy"));

  return ((!*min_dist_enforce || dist >= *min_dist) &&
          (!*max_dist_enforce || dist <= *max_dist) );
}



bool Retinal_ManagedComposite::accumulate_managed_next(bool val, Retinal_Object* r)
{
  /*
    Arbitrary "timeout" value to prevent endless loop; with
    some parameter values it may be impossible to find a valid
    distance.
  */
  const int max_trials=100; 

  /* Generate and test possible positions until a valid one is found */
  for (int trial=0; trial<max_trials; trial++) {
    r->next();
    
    iterator it = children.begin();
    while (*it != r && distance_valid(*r,**it))
      it++;

    if (*it == r) {
      r->set_active(true);
      return val;
    }
  }
  
  r->set_active(false);
  return false;
}



bool Retinal_ManagedComposite::next()
{
  bool val=true;
  for (iterator i = Retinal_Composite::children.begin();
       i != children.end();  i++) {
    val = accumulate_managed_next(val,*i);
  }

  return val && Retinal_Object::next();
}



/******************************************************************************/
/*  Retinal_AnchoredManagedComposite objects                                  */
/******************************************************************************/

Retinal_Object::Activity Retinal_AnchoredManagedComposite::activ(Coordinate x, Coordinate y) const
{
  return (dominant_child ? dominant_child->activation(x,y)
          : std::accumulate(ISEQ(children),accum_base(),ActivityAccumulator(x,y,accum_type)));
}



string Retinal_AnchoredManagedComposite::stringrep() const
{
  return name() + 
    (dominant_child ? " " + dominant_child->stringrep()
     : std::accumulate(ISEQ(children),string(""),accumulate_stringreps));
}



/******************************************************************************/
/*                                                                            */
/*  Miscellaneous leaf objects                                                */
/*  (See corresponding constants for documentation)                           */
/*                                                                            */
/*  Each usually caches parameter values that it needs (once) in update() and */
/*  then uses them (many times) in activ().  If it uses a bounding box, it    */
/*  usually computes the bounding box in update().                            */
/*                                                                            */
/******************************************************************************/

#ifndef NO_VALGEN_STRINGS
RetinalObjectStringArgs::ParamList Retinal_CircularGaussian::paramlist() {
  RetinalObjectStringArgs::ParamList p;
  p.push_back("cx");
  p.push_back("cy");
  p.push_back("xsigma");
  p.push_back("scale");
  p.push_back("offset");
  return p;
}
#endif


bool   Retinal_CircularGaussian::update() const
{
  /* Cache the result of sub-computations */
  cx                        = get_var("cx");
  cy                        = get_var("cy");
  const Variable xsigma     = get_var("xsigma");
  div_sigmasq               = 1/(xsigma*xsigma);

  /* Compute bounding box */
  const Coordinate boundrad=bound_mult*xsigma;
  bounding_box.set(cx-boundrad,cy-boundrad,
                   cx+boundrad,cy+boundrad);

  return Retinal_Object::update();
}


Retinal_Object::Activity Retinal_CircularGaussian::activ(Coordinate x, Coordinate y) const
{
  const Coordinate dx  = x-cx;
  const Coordinate dy  = y-cy;

  return exp( -(dx*dx+dy*dy)*div_sigmasq);
}



#ifndef NO_VALGEN_STRINGS
RetinalObjectStringArgs::ParamList Retinal_Gaussian::paramlist() {
  RetinalObjectStringArgs::ParamList p;
  p.push_back("theta");
  p.push_back("cx");
  p.push_back("cy");
  p.push_back("xsigma");
  p.push_back("ysigma");
  p.push_back("scale");
  p.push_back("offset");
  return p;
}
#endif


bool Retinal_Gaussian::update() const
{
  /* Cache the result of sub-computations */
  cx                        = get_var("cx");
  cy                        = get_var("cy");
  const Variable t          = get_var("theta");
  cost                      = cos(t);
  sint                      = sin(t);
  const Variable xsigma     = get_var("xsigma");
  const Variable ysigma     = get_var("ysigma");
  div_xsigma                = 1/xsigma;
  div_ysigma                = 1/ysigma;
  
  /* Compute bounding box */
  const Variable xrad=bound_mult*xsigma;
  const Variable yrad=bound_mult*ysigma;
  bounding_box.set(-xrad,-yrad,xrad,yrad).rotate(t).translate(cx,cy);

  return Retinal_Object::update();
}


Retinal_Object::Activity Retinal_Gaussian::activ(Coordinate x, Coordinate y) const
{
  const Coordinate dx    = x-cx;
  const Coordinate dy    = y-cy;
  const Coordinate xp    = ( dx * cost + dy * sint)*div_xsigma;
  const Coordinate yp    = (-dx * sint + dy * cost)*div_ysigma;

  return exp(-(xp*xp + yp*yp));
}



#ifndef NO_VALGEN_STRINGS
RetinalObjectStringArgs::ParamList Retinal_Rectangle::paramlist() {
  RetinalObjectStringArgs::ParamList p;
  p.push_back("cx");
  p.push_back("cy");
  p.push_back("xsigma");
  p.push_back("ysigma");
  p.push_back("scale");
  p.push_back("offset");
  return p;
}
#endif


bool Retinal_Rectangle::update() const
{
  /* Compute bounding box. */
  const Coordinate xrad   = get_var("xsigma");
  const Coordinate yrad   = get_var("ysigma");
  const Coordinate cx_    = get_var("cx");
  const Coordinate cy_    = get_var("cy");
  bounding_box.set(-xrad,-yrad,xrad,yrad).translate(cx_,cy_);

  return Retinal_Object::update();
}



#ifndef NO_VALGEN_STRINGS
RetinalObjectStringArgs::ParamList Retinal_SineGrating::paramlist() {
  RetinalObjectStringArgs::ParamList p;
  p.push_back("theta");
  p.push_back("freq");
  p.push_back("phase");
  p.push_back("scale");
  p.push_back("offset");
  p.push_back("cx");
  p.push_back("cy");
  return p;
}
#endif


bool Retinal_SineGrating::update() const
{
  /* Cache the result of sub-computations */
  const Variable t          = get_var("theta");
  cost                      = cos(t);
  sint                      = sin(t);
  phase                     = get_var("phase");
  freq                      = get_var("freq");
  cx                        = get_var("cx");
  cy                        = get_var("cy");

  return Retinal_Object::update();
}


Retinal_Object::Activity Retinal_SineGrating::activ(Coordinate x, Coordinate y) const
{  return sin( freq*((x-cx)*sint-(y-cy)*cost) + phase);  }


  
#ifndef NO_VALGEN_STRINGS  
RetinalObjectStringArgs::ParamList Retinal_Gabor::paramlist() {
  RetinalObjectStringArgs::ParamList p;
  p.push_back("theta");
  p.push_back("cx");
  p.push_back("cy");
  p.push_back("xsigma");
  p.push_back("ysigma");
  p.push_back("freq");
  p.push_back("phase");
  p.push_back("scale");
  p.push_back("offset");
  return p;
}
#endif


bool Retinal_Gabor::update() const
{
  /* Cache the result of sub-computations */
  cx                        = get_var("cx");
  cy                        = get_var("cy");
  const Variable t          = get_var("theta");
  cost                      = cos(t);
  sint                      = sin(t);
  phase                     = get_var("phase");
  freq                      = get_var("freq");
  const Variable xsigma     = get_var("xsigma");
  const Variable ysigma     = get_var("ysigma");
  div_xsigmasq              = 1/(xsigma*xsigma);
  div_ysigmasq              = 1/(ysigma*ysigma);

  /* Compute bounding box */
  const Variable xrad=bound_mult*xsigma;
  const Variable yrad=bound_mult*ysigma;
  bounding_box.set(-xrad,-yrad,xrad,yrad).rotate(t).translate(cx,cy);
  
  return Retinal_Object::update();
}


Retinal_Object::Activity Retinal_Gabor::activ(Coordinate x, Coordinate y) const
{
  const Coordinate dx       = x-cx;
  const Coordinate dy       = y-cy;
  const Coordinate xp       = ( dx * cost + dy * sint);
  const Coordinate yp       = (-dx * sint + dy * cost);

  const Activity   gaussian = exp(-(xp*xp*div_xsigmasq + yp*yp*div_ysigmasq));
  const Activity   grating  = 0.5+0.5*cos( freq*yp + phase );
  
  return grating * gaussian;
}



#ifndef NO_VALGEN_STRINGS  
RetinalObjectStringArgs::ParamList Retinal_FuzzyLine::paramlist() {
  RetinalObjectStringArgs::ParamList p;
  p.push_back("theta");
  p.push_back("cx");
  p.push_back("cy");
  p.push_back("ysigma");
  p.push_back("center_width");
  p.push_back("scale");
  p.push_back("offset");
  return p;
}
#endif


bool Retinal_FuzzyLine::update() const
{
  /* Cache the result of sub-computations */
  cx                        = get_var("cx");
  cy                        = get_var("cy");
  centerw                   = get_var("center_width");
  const Variable t          = get_var("theta");
  cost                      = cos(t);
  sint                      = sin(t);
  const Variable ysigma     = get_var("ysigma");
  div_ysigmasq              = 1/(ysigma*ysigma);

  return Retinal_Object::update();
}


Retinal_Object::Activity Retinal_FuzzyLine::activ(Coordinate x, Coordinate y) const
{
  const Coordinate dx                 = x-cx;
  const Coordinate dy                 = y-cy;
  const Coordinate distance_from_line = fabs(dy*cost - dx*sint);
  const Coordinate gaussian_x_coord   = (distance_from_line - centerw/2);
  return (gaussian_x_coord<=0 ? 1.0
          : exp(-gaussian_x_coord*gaussian_x_coord*div_ysigmasq));
}



#ifndef NO_VALGEN_STRINGS  
RetinalObjectStringArgs::ParamList Retinal_FuzzyDisc::paramlist() {
  RetinalObjectStringArgs::ParamList p;
  p.push_back("cx");
  p.push_back("cy");
  p.push_back("ysigma");
  p.push_back("center_width");
  p.push_back("scale");
  p.push_back("offset");
  return p;
}
#endif


bool Retinal_FuzzyDisc::update() const
{
  /* Cache the result of sub-computations */
  cx                        = get_var("cx");
  cy                        = get_var("cy");
  centerw                   = get_var("center_width");
  const Variable ysigma     = get_var("ysigma");
  div_ysigmasq              = 1/(ysigma*ysigma);

  return Retinal_Object::update();
}


Retinal_Object::Activity Retinal_FuzzyDisc::activ(Coordinate x, Coordinate y) const
{
  const Coordinate dx                 = x-cx;
  const Coordinate dy                 = y-cy;
  const Coordinate distance_from_line = sqrt(dx*dx+dy*dy);
  const Coordinate gaussian_x_coord   = (distance_from_line - centerw/2);
  return (gaussian_x_coord<=0 ? 1.0
          : exp(-gaussian_x_coord*gaussian_x_coord*div_ysigmasq));
}



#ifndef NO_VALGEN_STRINGS  
RetinalObjectStringArgs::ParamList Retinal_FuzzyRing::paramlist() {
  RetinalObjectStringArgs::ParamList p;
  p.push_back("cx");
  p.push_back("cy");
  p.push_back("ysigma");
  p.push_back("center_width");
  p.push_back("radius");
  p.push_back("scale");
  p.push_back("offset");
  return p;
}
#endif


bool Retinal_FuzzyRing::update() const
{
  /* Cache the result of sub-computations */
  cx                        = get_var("cx");
  cy                        = get_var("cy");
  centerw                   = get_var("center_width");
  radius                    = get_var("radius");
  const Variable ysigma     = get_var("ysigma");
  div_ysigmasq              = 1/(ysigma*ysigma);

  return Retinal_Object::update();
}


Retinal_Object::Activity Retinal_FuzzyRing::activ(Coordinate x, Coordinate y) const
{
  const Coordinate dx                 = x-cx;
  const Coordinate dy                 = y-cy;
  const Coordinate distance_from_line = fabs(radius - sqrt(dx*dx+dy*dy));
  const Coordinate gaussian_x_coord   = (distance_from_line - centerw/2);
  return (gaussian_x_coord<=0 ? 1.0
          : exp(-gaussian_x_coord*gaussian_x_coord*div_ysigmasq));
}



/******************************************************************************/
/*  Retinal_PGM objects                                                       */
/******************************************************************************/

#ifndef NO_VALGEN_STRINGS

string pgm_full_pathname(const string& filename, const string& paths="")
{
  StringParser::arglist words;
  const StringParser p;
  p.parse(paths,words);
  
  StringParser::arglist::iterator i=words.begin();
  FILE *file=NULL;
  string name=filename;
  
  /* First try the current directory, then try all the possible paths */
  for (;;) {
    file=fopen(name.c_str(),"r");
    if (file) {
      fclose(file);
      return name;
    }
    if (i==words.end())
      return filename; /* Couldn't find file */
    name = (*i++)+filename;
  }
#ifdef __GNUC__
  return ""; /* Never reached; just silences warning */
#endif
}


RetinalObjectStringArgs::ParamList Retinal_PGM::paramlist() {
  RetinalObjectStringArgs::ParamList p;
  p.push_back("theta");
  p.push_back("cx");
  p.push_back("cy");
  p.push_back("size_scale");
  p.push_back("autoscale");
  p.push_back("transparent");
  p.push_back("scale");
  p.push_back("offset");
  return p;
}


Retinal_PGM::Retinal_PGM( const string& name_val, RetinalObjectStringArgs& sa,
                          Retinal_Composite* parent)
  : Retinal_Object(name_val), ascale(1.0), aoffset(0.0), transparent_level(-1)
{
  const string filename = sa.parsed_next("image_filename");
  const string paths    = sa.parsed_get_default("image_paths");
  basename=pgm_full_pathname(filename, paths);
  const std::pair<int,int> size = pnm_size(basename);
  const int width  = size.first;
  const int height = size.second;
  
  /* Read the image itself into memory if valid */
  if (width<=0 || height<=0) {
    sa.error("Could not read image file `"+basename+"'");
    return;
  }
  set_active(pgm_read(basename,image));
  if (!get_active())
    sa.error("Problem reading image file `"+basename+"'; not in .pgm or .pbm format?");
  border=mat::edge_average(image);

  /*
    Compute special default parameter values
    
    Images smaller than the retina will always have at least one
    border showing, so we assume that all borders are ok.  If that's
    true, the picture can be placed and rotated just like any other
    object.  Larger images are assumed to represent the entire retina,
    so their position range is cropped to fully cover the retina at
    all times, never showing a border.  All such images default to
    being upright (with the width and height assuming that
    orientation), because rotating them would show a border unless the
    positions were cropped to the axis-aligned rectangle inscribed
    within the largest circle fitting in the ellipse inscribed within
    the image border.  (Even if that wouldn't be as complicated to
    calculate as to describe, it would certainly limit the usable area
    of the picture.)

    Note: The defaults will not work if the user specifies a size_scale.
  */
  if (parent) {
    double dummy;
    const string size_scale_str = sa.parsed_get_default("size_scale");
    const double size_scale     = sa.parse(size_scale_str,dummy);
    const double parentwidth    = parent->nominal_width();
    const double parentheight   = parent->nominal_height();
    const double excesswidth    =  width*size_scale - parentwidth;
    const double excessheight   = height*size_scale - parentheight;

    if (excesswidth>=0 && excessheight>=0) {
      /* The formulas for the Uncorrelation bounds are not yet correct. */
      const string cx_default     = "Uncorrelate &uncorrelation " + String::stringrep(parentwidth/2-excesswidth/2) + " " + String::stringrep(parentwidth/2+excesswidth/2) + " Random " + String::stringrep(parentwidth/2)
        + " " + String::stringrep(excesswidth/2);
      const string cy_default     = "Uncorrelate &uncorrelation " + String::stringrep(parentheight/2-excessheight/2) + " " + String::stringrep(parentheight/2+excessheight/2) + " Random " + String::stringrep(parentheight/2)
        + " " + String::stringrep(excessheight/2);
      
      sa.set_default("cx",    cx_default);
      sa.set_default("cy",    cy_default);
      sa.set_default("theta", "PI/2");
    }
  }

  /* Process the user-supplied arguments using the constructed set of defaults */
  merge_vars(sa.vars(paramlist()));
  
  /* Set up transparency, if any */
  const Variable transparent = get_var("transparent");
  if (transparent>0)
    transparent_level=transparent;
  else if (pgm_white_is_transparent(basename))
    transparent_level=0.999; /* Slightly larger than 254/255 (0.996) */

  //cout << "[" << "Transparent: " << transparent_level << "(" << transparent << ")]" << endl;

  /* To debug the computed default values */
  //cout << "Retinal_PGM -- file " << filename << " will use cx '"
  //         << vargen_stringrep("cx") << "' and cy '"
  //         << vargen_stringrep("cy") << "'" << endl;
  
  /* Automatic scaling */
  const Variable autoscale  = get_var("autoscale");
  if (autoscale>0) {
    const double maxval = *(std::max_element(MSEQ(image)));
    const double minval = *(std::min_element(MSEQ(image)));
    const double range  = maxval-minval;
    if (range>0) {
      ascale  = 1/range;
      aoffset = minval;
      border  = border*ascale+aoffset;
      //cout << "Autoscaling to " << range << " (range " << minval << ".." << maxval << ")" << endl;
    }
  }
}
#endif /* NO_VALGEN_STRINGS */


bool Retinal_PGM::update() const
{
  /* Cache the result of sub-computations */
  const Variable mt         = -get_var("theta");
  cx                        =  get_var("cx");
  cy                        =  get_var("cy");
  cosmt                     = cos(mt);
  sinmt                     = sin(mt);
  const Variable size_scale = get_var("size_scale");
  div_size                  = 1/size_scale;
  /* Image is stored normally but rotated 90 degrees for drawing;
     see Retinal_PGM::activ.  Thus the bounds are switched here. */
  xoff                      = image.nrows()/2;
  yoff                      = image.ncols()/2;
    
  /* Compute bounding box */
  bounding_box.set(-xoff-1,-yoff-1,+xoff+1,+yoff+1); /* Generous to avoid truncation */
  bounding_box.scale(size_scale,size_scale).rotate(mt).translate(cx,cy);

  return Retinal_Object::update();
}



Retinal_Object::Activity Retinal_PGM::pixel_value(Coordinate x, Coordinate y, bool& inbounds) const
{
  const Coordinate dx = (x-cx)*div_size;
  const Coordinate dy = (y-cy)*div_size;
  const Coordinate xp = dx*cosmt-dy*sinmt+xoff;
  const Coordinate yp = dx*sinmt+dy*cosmt+yoff;
      
  /*
    Even though the image is stored in the usual orientation, it is
    rotated for drawing so that the top points right for theta=0; that
    way, it will be upright given the usual PI/2 orientation default.
  */
  const int r = image.nrows()-int(xp+0.5)-1;
  const int c = image.ncols()-int(yp+0.5)-1;
 
  inbounds = ( r >= 0 && c >= 0 && r < int(image.nrows()) && c < int(image.ncols()) );
  return inbounds? image[r][c] : 0;
}



Retinal_Object::Activity Retinal_PGM::activ(Coordinate x, Coordinate y) const
{
  bool inbounds;
  Retinal_Object::Activity pixval=pixel_value(x,y,inbounds);
  return inbounds? (ascale * (pixval-aoffset)) : border;
  //: 0.75; /* To debug image boundaries */ 
}


bool Retinal_PGM::inside(Coordinate x, Coordinate y) const
{
  /* Quick check */
  if (!bounding_box.inside(x,y)) return false;

  /* W/o transparency, bbox determines inside() */
  if (transparent_level<0) return true;
  
  /* W/ transparency, a slower, pixel-by-pixel check against
     transparency threshold is needed */
  bool inbounds;
  Retinal_Object::Activity pixval=pixel_value(x,y,inbounds);
  return (inbounds && !(pixval>transparent_level));
}

---