/*
 * Filename: percept.c
 *
 * Author: Martin Reddy, <reddy@ai.sri.com>  - 20 Dec 2000
 *
 * Purpose:
 *   This program reads an input image (binary PPM format) and processes
 *   this to create another image showing the degree of detail actually
 *   visible to the human eye, assuming that the image is presented 
 *   across a certain extent of our eyes' field of view which defaults
 *   to the entire visual field.
 *
 *   Essentially, we calculate the highest spatial frequency of each
 *   pixel given its eccentricity (degree into the visual periphery) and
 *   a constant user-specified velocity (degrees per second). We then
 *   blur this pixel using an equivalently-sized Gaussian filter.
 *
 *   The model used to evaluate spatial frequency limit is described in:
 *
 *   M. Reddy (1997). "Perceptually Modulated Level of Detail for Virtual
 *      Environments". Phd Thesis. University of Edinburgh.
 *
 * Status:
 *   This code is released as Open Source under the GNU General Public
 *   License (GPL), http://www.opensource.org/licenses/gpl-license.html
 *
 *   (c) 2000,2001 Copyright Martin Reddy. All rights reserved.
 *
 */

#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>

#define BUFSIZE 256
#define max(a,b) (((a)>(b))?(a):(b))
#define min(a,b) (((a)<(b))?(a):(b))

/* just the one global variable... to control output of debug information */

int debug = 0;

/* structure to hold an image and all parameters */

typedef unsigned char uchar;

typedef struct {
  uchar *buffer;                    /* bitmap for the image */
  int   w, h, c;                    /* height, width, and components */
  float hfov, vfov, hoff, voff;     /* field of view and fovea offset (deg) */
} Image;

typedef struct {
  float r, g, b;                    /* red, green, blue components [0..255] */
} Pixel;

/*
 * log_debug
 *
 * Outputs a printf() style message to the log file, percept.log, but
 * only if the debug flag is on. This routine will manage opening the
 * file the first time it is called, and will close the file if passed
 * a NULL format string
 *
 */

void
log_debug( const char *format, ... )
{
  static FILE *fp = NULL;
  va_list     arg;

  if ( ! debug ) return;

  if ( ! format ) {
    if ( fp ) fclose( fp );
    fp = NULL;
    return;
  }

  if ( ! fp ) {
    if ( ! ( fp = fopen( "percept.log", "w" ) ) ) {
      fprintf( stderr, "Cannot create debug log file!\n" );
      debug = 0;
      return;
    }
  }

  va_start( arg, format );
  vfprintf( fp, format, arg );
  fflush( fp );
  va_end( arg );
}

/*
 * alloc_image( w, h )
 *
 * Allocates an Image structure with a specified width and height
 * and allocates the memory for an RGB bitmap to fit this size
 *
 */

Image *
alloc_image( int w, int h )
{
  Image *image = calloc( sizeof(Image), 1 );
  if ( ! image ) return NULL;
  image->w = w;
  image->h = h;
  image->c = 3;
  /* binocular visual field from "Useful Numbers from Vision Science" */
  image->hfov = 200.0;   /* binocular horizontal visual field */
  image->vfov = 135.0;   /* binocular vertical visual field */
  image->hoff = 0.0;
  image->voff = 0.0;
  image->buffer = calloc( w * h, 3 );
  if ( ! image->buffer ) {
    free( image );
    return NULL;
  }
  return image;
}

/*
 * get_pixel( image, x, y )
 *
 * Returns the RGB pixel at the specified (x,y) pixel coordinate,
 * or returns NULL if the coordinate is invalid.
 *
 */

Pixel *
get_pixel( Image *image, int x, int y )
{
  static Pixel pixel;
  int          offset = ( y * image->w * image->c ) + ( x * image->c );
  
  if ( image == NULL || image->buffer == NULL ) return NULL;
  if ( x < 0 || x >= image->w ) return NULL;
  if ( y < 0 || y >= image->h ) return NULL;

  pixel.r = (float) image->buffer[offset];
  pixel.g = (float) image->buffer[offset+1];
  pixel.b = (float) image->buffer[offset+2];

  return &pixel;
}

/*
 * set_pixel( image, x, y, pixel )
 *
 * Takes the specified RGB value and writes it to the pixel at
 * coordinate (x,y) in the image. No action is taken if an invalid
 * pixel coordinate or image are specified.
 *
 */

void
set_pixel( Image *image, int x, int y, Pixel *pixel )
{
  int offset = ( y * image->w * image->c ) + ( x * image->c );
  if ( pixel == NULL || image == NULL || image->buffer == NULL ) return;
  if ( x < 0 || x >= image->w ) return;
  if ( y < 0 || y >= image->h ) return;
  image->buffer[offset]   = (uchar) pixel->r;
  image->buffer[offset+1] = (uchar) pixel->g;
  image->buffer[offset+2] = (uchar) pixel->b;
}

/*
 * read_image( filename )
 * 
 * Takes the filename for a binary PPM image, reads this image into
 * memory, and returns an Image structure all filled out. Returns
 * NULL if an error was encountered
 *
 */

Image *
read_image( char *filename )
{
  FILE  *fp = fopen( filename, "r" );
  Image *image;
  char  buffer[BUFSIZE];
  int   w, h;

  if ( fp == NULL ) {
    fprintf( stderr, "Couldn't open image %s\n", filename );
    return NULL;
  }

  fgets( buffer, BUFSIZE, fp );
  if ( strcmp( buffer, "P6\n" ) != 0 ) {
    fprintf( stderr, "File is not a binary PPM %s\n", filename );
    fclose( fp );
    return NULL;
  }

  do { 
    fgets( buffer, BUFSIZE, fp );
  } while ( buffer[0] == '#' );

  sscanf( buffer, "%d %d", &w, &h );
  fgets( buffer, BUFSIZE, fp );

  if ( ( image = alloc_image( w, h ) ) == NULL ) {
    fprintf( stderr, "Couldn't allocate memory for image\n" );
    fclose( fp );
    return NULL;
  }
 
  if ( fread( image->buffer, 1, w*h*3, fp ) <= 0 ) {
    fclose( fp );
    return NULL;
  }

  fclose( fp );
  return image;
}

/*
 * write_image( image, filename, vel, no_ecc )
 *
 * Takes an Image structure and writes out the RGB bitmap to a
 * binary PPM image using the specified filename. The vel and no_ecc
 * parameters are used to create a comment in the PPM header describing
 * the velocity and eccentricity settings appropriate to this image.
 *
 */

void
write_image( Image *image, char *filename, float vel, int no_ecc )
{
  FILE *fp = fopen( filename, "w" );
  if ( ! fp ) return;
  fprintf( fp, "P6\n# percept: vel %g, ecc %s, hfov %g vfov %g",
	   vel, (no_ecc)?"off":"on", image->hfov, image->vfov );
  fprintf( fp, ", hoff %g, voff %g\n", image->hoff, image->voff );
  fprintf( fp, "%d %d\n255\n", image->w, image->h );
  if ( fwrite( image->buffer, 1, image->w * image->h * image->c, fp ) <= 0 )
    fprintf( stderr, "Error while writing image!\n" );
  fclose( fp );
}

/*
 * get_arg_float( argc, argv, opt, value )
 * get_arg_string( argc, argv, opt, value )
 * get_arg( argc, argv, opt )
 *
 * A number of functions for parsing command line options.
 *
 */

void
get_arg_float( int argc, char **argv, char *opt, float *value )
{
  int i;
  for ( i = 1; i < argc; i++ ) {
    if ( strcmp( argv[i], opt ) == 0 && i < argc-1 ) {
      *value = (float) atof( argv[++i] );
      return;
    }
  }
}

void
get_arg_string( int argc, char **argv, char *opt, char **value )
{
  int i;
  for ( i = 1; i < argc; i++ ) {
    if ( strcmp( argv[i], opt ) == 0 && i < argc-1 ) {
      *value = argv[++i];
      return;
    }
  }
}

int
get_arg( int argc, char **argv, char *opt )
{
  int i;
  for ( i = 1; i < argc; i++ )
    if ( strcmp( argv[i], opt ) == 0 )
      return 1;
  return 0;
}

/*
 * gaussian2d( xvalue, yvalue, xsize, ysize )
 *
 * calculates a 2D Gaussian scaling factor. You specify a location
 * (xvalue, yvalue) inside an (xsize, ysize) area, and the function
 * will return a scaling factor from 0 to 1.0, where 1.0 is at the
 * center of the area.
 *
 */

float
gaussian2d( float xvalue, float yvalue, float xsize, float ysize )
{
  float xgauss, ygauss;
  xsize /= 3.0;  ysize /= 3.0;
  xgauss = (float) exp( -( xvalue * xvalue ) / ( xsize * xsize ) );
  ygauss = (float) exp( -( yvalue * yvalue ) / ( ysize * ysize ) );
  return xgauss * ygauss;
}

/*
 * calc_pixel( image, x, y, deg_size, value )
 *
 * Calculates the RGB value for a pixel at (x, y) given a visibility
 * limit of deg_size. This degree extent is used to calculate the
 * size of a Gaussian kernel to apply a low-pass filter at that
 * pixel location. The result is to blur the pixel by smoothing across
 * a number of pixels in propotion to the visiblity limit
 *
 */

void
calc_pixel( Image *image, int x, int y, float deg_size, Pixel *value )
{
  float xspanf = deg_size * image->w / image->hfov;
  float yspanf = deg_size * image->h / image->vfov;
  int   xspan  = (int) ceil( xspanf / 2.0f );
  int   yspan  = (int) ceil( yspanf / 2.0f );
  float total  = 0.0, gauss;
  int   i, j;
  Pixel *pixel;

  value->r = 0.0f;
  value->g = 0.0f;
  value->b = 0.0f;

  log_debug( "calc_pixel(%d, %d), deg_size = %1.4f, span(%1.1f, %1.1f)\n",
	     x, y, deg_size, xspanf, yspanf );

  for ( i = -yspan; i <= yspan; i++ ) {
    for ( j = -xspan; j <= xspan; j++ ) {
      pixel = get_pixel( image, x+j, y+i );
      if ( ! pixel ) continue;
      gauss = gaussian2d( (float) j, (float) i, xspanf, yspanf );
      pixel->r *= gauss;
      pixel->g *= gauss;
      pixel->b *= gauss;
      total += gauss;
      value->r += pixel->r;
      value->g += pixel->g;
      value->b += pixel->b;
    }
  }

  value->r = min( max( value->r / total, 0.0f ), 255.0f );
  value->g = min( max( value->g / total, 0.0f ), 255.0f );
  value->b = min( max( value->b / total, 0.0f ), 255.0f );

  log_debug( "gaussian(%d, %d; %1.4f deg) gives (%1.0f, %1.0f, %1.0f)\n",
	     x, y, deg_size, value->r, value->g, value->b );
}

/*
 * usage()
 *
 * Displays the program usage to stdout and exits
 *
 */

void
usage( void )
{
  puts( "Percept V1.0 - Martin Reddy, 30 Dec 2000\n" );
  puts( "usage: percept <ppm_image> [options]" );
  puts( "where," );
  puts( "  -vfov <deg> = vertical field of view (deg)" );
  puts( "  -hfov <deg> = horizontal field of view (deg)" );
  puts( "  -voff <deg> = vertical offset of sweet spot (deg)" );
  puts( "  -hoff <deg> = horizontal offset of sweet spot (deg)" );
  puts( "  -vel <degs> = velocity (deg/s)" );
  puts( "  -noecc = assume constant eccentricity of 0 deg throughout" );
  puts( "  -debug = output debug details to percept.log" );
  puts( "  -o <name> = output ppm filename (default is percept.ppm)" );
  puts( "" );
  exit( 0 );
}

/*
 * show_percentage( image, x, y )
 *
 * Displays a progress report with a percentage done readout
 *
 */

void
show_percentage( Image *image, int x, int y )
{
  static float last_perc = -2.0f;
  float pixels = (float) image->w * image->h;
  float perc = 100.0f * ( ( y * image->w ) + x ) / pixels;
  if ( perc > last_perc + 1.0f ) {
    fprintf( stderr, "Processing (%1.0f%%)...  \r", perc );
    last_perc = perc;
  }
}

/*
 * threshold_spatial_frequency( ecc, vel )
 *
 * Calculates the threshold spatial frequency of the human visual
 * system to a stimulus at a set eccentricity (deg) and moving at a
 * specified velocity (deg/s). The result is in units of c/deg.
 *
 */

float
threshold_spatial_frequency( float ecc, float vel )
{
  float motionSF = 60.0f, cortMag = 1.0f;

  if ( vel > 0.825f )
    motionSF = -27.78f * log10( vel ) + 57.69f;
  
  if ( ecc > 5.79f ) {
    float temp = ( 0.3f * ecc ) + 1.0f;
    cortMag = 7.49f / ( temp * temp );
  }

  return motionSF * cortMag;
}

/*
 * max_size_at_pixel( image, x, y, vel, no_ecc )
 *
 * Calculates the maxium size of stimulus (deg) that can be resolved
 * at the specified pixel coordinate. A velocity parameter can be
 * specified, and the eccentricity parameter can be selectively ignored.
 *
 */

float
max_size_at_pixel( Image *image, int x, int y, float vel, int no_ecc )
{
  float xps, yps, xdoff, ydoff, ecc, sf;

  if ( no_ecc ) 
    ecc = 0.0;
  else {
    /* calculate the pixel coordinate of the sweet spot */
    xps   = ( image->w / 2.0f ) + ( image->hoff * image->w ) / image->hfov;
    yps   = ( image->h / 2.0f ) + ( image->voff * image->h ) / image->vfov;
    /* calculate the degree offset from the sweet spot */
    xdoff = ( x - xps ) * image->hfov / image->w;
    ydoff = ( y - yps ) * image->vfov / image->h;
    /* calculate the degree distance from the sweet spot (eccentricity) */
    ecc   = (float) sqrt( ( xdoff * xdoff ) + ( ydoff * ydoff ) );
  }

  /* calculate the threshold spatial frequency at this ecc and vel */
  sf    = threshold_spatial_frequency( ecc, vel );

  log_debug( "pixel(%d, %d), ecc = %1.2f, vel = %1.2f, spatfreq = %1.4f\n",
       x, y, ecc, vel, sf );

  /* return the threshold stimulus size in degrees */
  return ( sf == 0.0f ) ? 0.0f : 1.0f / ( sf * 2.0f );
}

/*
 * main( argc, argv )
 *
 * The main program which parses the command line arguments, reads the
 * input image, applies the perceptually modulated Gaussian blurring,
 * and outputs the transformed image.
 *
 */

int
main( int argc, char **argv )
{
  Image *in_image, *out_image;
  float vel = 0.0f, size;
  Pixel pixel;
  int   x, y, no_ecc = 0;
  char  *out_name = "percept.ppm";

  /* read the input image */

  if ( argc < 2 ) usage();

  if ( ( in_image = read_image( argv[1] ) ) == NULL ) {
    fprintf( stderr, "Could not read the inputs image\n" );
    exit( 1 );
  }

  /* parse the user's command line arguments */

  if ( get_arg( argc, argv, "-noecc" ) )
    no_ecc = 1;
  if ( get_arg( argc, argv, "-debug" ) )
    debug = 1;

  get_arg_float( argc, argv, "-hfov", &in_image->hfov );
  get_arg_float( argc, argv, "-vfov", &in_image->vfov );
  get_arg_float( argc, argv, "-hoff", &in_image->hoff );
  get_arg_float( argc, argv, "-voff", &in_image->voff );
  get_arg_float( argc, argv, "-vel", &vel );

  get_arg_string( argc, argv, "-o", &out_name );

  if ( get_arg( argc, argv, "-calc" ) ) {
    float ecc = 0.0, sf;
    get_arg_float( argc, argv, "-calc", &ecc );
    sf = threshold_spatial_frequency( ecc, vel );
    printf( "threshold sf at %g deg (%g deg/s) = %g c/deg\n", ecc, vel, sf );
    exit( 0 );
  }

  /* create some header debug output */

  log_debug( "Debug output for percept utility\n\n" );
  log_debug( "image dimensions = %d x %d pixels\n", in_image->w, in_image->h );
  log_debug( "field of view = %1.1f x %1.1f deg\n", in_image->hfov,
	     in_image->vfov );
  log_debug( "sweet spot offset = %1.1f x %1.1f deg\n", in_image->hoff,
	     in_image->voff );
  log_debug( "eccentricity degradation = %s\n", ( no_ecc ) ? "off" : "on" );
  log_debug( "global velocity = %1.1f deg/s\n", vel );
  log_debug( "pixel frequency = %1.4f x %1.4f c/deg\n\n",
	     in_image->hfov / in_image->w,in_image->vfov / in_image->h );

  /* create the output bitmap */

  if ( ! ( out_image = alloc_image( in_image->w, in_image->h ) ) ) {
    fprintf( stderr, "Could not allocate output bitmap\n" );
    exit( 1 );
  }

  /* go through each pixel and work out how much we need to blur it */

  for ( y = 0; y < in_image->h; y++ ) {
    for ( x = 0; x < in_image->w; x++ ) {
      show_percentage( in_image, x, y );
      size = max_size_at_pixel( in_image, x, y, vel, no_ecc );
      calc_pixel( in_image, x, y, size, &pixel );
      set_pixel( out_image, x, y, &pixel );
    }
  }

  /* close the log file */

  log_debug( NULL );

  /* write out the processed image */
  
  write_image( out_image, out_name, vel, no_ecc );

  printf( "Saved image to %s.\n", out_name );
  return 0;
}

/* EOF: percept.c */