#define EXPONENT 4.7

/* the following two are no longer used */
#define WA 1.0
#define WB 1.0

#define R1 1
#define G1 1
#define B1 1
#define R2 1
#define G2 0
#define B2 1

/*
 * WEIGHTED PHOTOQUANTIGRAPHIC SUM OF TWO IMAGES
 *
 * See also, http://wearcam.org/lightspaces/index.html
 * etc...
 *
 * Related references: Proc. IEEE Nov. 1998, etc.. http://wearcam.org
 *
 * photoquantigraphic image compositing add program
 * photoquantigraphically adds two images
 * with specification of weighting
 *
 * bugs or shortcomings: still need to allow WA and WB to be RGB sets,
 * e.g. weights for each of red, green, and blue separately
 *
 * pnmbmath.c - law of composition on 2 images
 * to compile: gcc -Wall pnmbmath.c -o pnmbmath -lm
 * to run: pnmbmath v080.ppm v081.ppm -o pork.ppm
 *
 * steve; from ndf's diff.c, with help of sbeck on parsing input images
 */

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

struct image_params {

  char * filename;
  int width;
  int height;
  /* max_val is the maximum colour component value - often 255 */
  int max_val;
  /* type is either '5' to indicate a P5 file or '6' to
     indicate a P6 file */
  int type;
};

double pow(double x, double y);

void
usage(void)
{
  fprintf(stderr, "Use: pnmbmath image1 image2 [-o imageout]\n");
  exit(EXIT_SUCCESS);
}

void
parse_commandline(int argc, char ** argv, struct image_params * a_params,
		  struct image_params * b_params, char ** output_filename)
{
  int i;

  if (argc<3) usage();

  for (i=1;i<argc;i++)
    {
      if (!strcasecmp(argv[i], "-o"))
	{
	  if (argc<=(i+1)) usage();
	  *output_filename=(char *)strdup(argv[++i]);
	}
      else if (a_params->filename==NULL)
	a_params->filename=(char *)strdup(argv[i]);
      else if (b_params->filename==NULL)
	b_params->filename=(char *)strdup(argv[i]);
      else usage();
    }

  if ((a_params->filename==NULL) || (b_params->filename==NULL)) usage();
}

void
get_image_params(FILE * ifs, struct image_params * params)
{
  int c;
  int dims[3];
  int dims_scanned;
  if( ifs == NULL ) { /* check for file pointers to non existant files */
    fprintf(stderr,"get_frame_params: non-existant file\n");
    exit(EXIT_FAILURE);
  }

  if (fgetc(ifs)!='P') /* pnm file must begin with P, e.g. P5 or P6 */
    {
      fprintf(stderr, "get_frame_params: input does not begin with ``P''"
	              " thus not a pnm image.\n");
      exit(EXIT_FAILURE);
    }
  c=fgetc(ifs);
  if (c!='5' && c!='6')
    {
      fprintf(stderr, "get_frame_params: input file must begin with P5 or P6.\n");
      exit(EXIT_FAILURE);
    }
  params->type=c;

  /* enter a loop
   * skip comments which start with #
   * and look for the dimensions
   */
  dims_scanned = 0;
  do {
    char line[80];
    fgets( line, sizeof(line), ifs );
    fprintf(stderr,"got line [%s]\n", line );
    if( *line != '#' ) {
      char *p = strtok(line," \n" );
      while( p && (*p != '#') ) {
        sscanf(p,"%d", &dims[dims_scanned] );
        dims_scanned++;
        p = strtok( NULL, " \n" );
      }
    }
    fprintf(stderr,"dims_scanned:%d\n", dims_scanned);
  } while( dims_scanned < 3 );

  params->width = dims[0];
  params->height = dims[1];
  params->max_val = dims[2];

  if( params->max_val != 255 ){
    fprintf(stderr,"get_frame_params: max value must be 255\n");
    exit( EXIT_FAILURE );
  }

}

int
main(int argc, char ** argv)
{
  struct image_params a_params={NULL,0,0,0,0};
  struct image_params b_params={NULL,0,0,0,0};
  char * dest_filename=NULL;
  FILE * a, * b, * dest;
  double p=EXPONENT;  /* exponent for pow */
  double p_inv;  /* 1.0/p */
  double scale_factor; /* to make the result on same interval is input */
  double wa, wb; /* weights for image a and image b */
  unsigned char a_pixel[3], b_pixel[3];
  unsigned char same_pixel[3]= {0,0,0};
  unsigned char out_pixel[3];
  int bytes_per_pixel, i;

  p_inv=1.0/p;
  scale_factor = pow(2.0,p_inv);

  parse_commandline(argc, argv, &a_params, &b_params, &dest_filename);

  if ((a=fopen(a_params.filename, "r"))==NULL)
    {
      fprintf(stderr, "Unable to open %s.\n", a_params.filename);
      exit(EXIT_FAILURE);
    }

  if ((b=fopen(b_params.filename, "r"))==NULL)
    {
      fprintf(stderr, "Unable to open %s.\n", b_params.filename);
      exit(EXIT_FAILURE);
    }

  if (dest_filename==NULL) dest=stdout;
  else if ((dest=fopen(dest_filename, "w"))==NULL)
    {
      fprintf(stderr, "Unable to open %s.\n", dest_filename);
      exit(EXIT_FAILURE);
    }

  get_image_params(a, &a_params);
  get_image_params(b, &b_params);

  if ((a_params.width!=b_params.width) || (a_params.height!=b_params.height)
      || (a_params.max_val!=b_params.max_val))
    {
      fprintf(stderr, "Images must have the same dimensions and maximum "
	      "colour values.\n");
      exit(EXIT_FAILURE);
    }

  if (a_params.type!=b_params.type)
    {
      fprintf(stderr, "%s is a P%c but %s is a P%c.\n", argv[1],
	      a_params.type, argv[2], b_params.type);
      exit(EXIT_FAILURE);
    }

  if (a_params.type=='5') bytes_per_pixel=1;
  else bytes_per_pixel=3;

  *out_pixel=out_pixel[1]=out_pixel[2]=a_params.max_val;

  fprintf(dest, "P%c\n%d %d\n%d\n", a_params.type, a_params.width,
	  a_params.height, a_params.max_val);

  while (!feof(a))
    {
      fread(a_pixel, sizeof(unsigned char), bytes_per_pixel, a);
      fread(b_pixel, sizeof(unsigned char), bytes_per_pixel, b);
      if (memcmp(a_pixel, b_pixel, sizeof(unsigned char)*bytes_per_pixel))
	{
	  for (i=0;i<bytes_per_pixel;i++) {
            if (i==0) {wa=R1; wb=R2;}
            if (i==1) {wa=G1; wb=G2;}
            if (i==2) {wa=B1; wb=B2;}
            if ((wa>0) || (wb>0)) {
	      out_pixel[i]=(unsigned char)
                (
                  pow(
                      (
                       pow((double)a_pixel[i],p)*wa
                       +
                       pow((double)b_pixel[i],p)*wb
                      )/(wa+wb)
                  ,p_inv)/scale_factor
                );
            }
            else out_pixel[i]=0.0;
            if (out_pixel[i] > 255.0) out_pixel[i]=255.0;
            if (out_pixel[i] < 0.0) out_pixel[i]=0.0;
          }
	  fwrite(out_pixel, sizeof(unsigned char), bytes_per_pixel, dest);
	}
      else
	fwrite(same_pixel, sizeof(unsigned char), bytes_per_pixel, dest);
    }

  fclose(a); fclose(b); fclose(dest);
  exit(EXIT_SUCCESS);
}