#define EXPONENT 4.7
#define WA 1.0
#define WB 1.0

#define R1 0.0
#define G1 0.3
#define B1 1.0
#define R2 1.0
#define G2 1.0
#define B2 0.5

#define MAX_NUM_OF_IMAGES 100
/*
 * 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;

  float pix_color_weight[2];
};

double pow(double x, double y);

double certainty_function(double pixelvalue)
{

  // doesn't work:

  //double x = 6*(pixelvalue/255) -3;
 
  //printf("pic=%f  x=%f\n", pixelvalue, x);

  //return pow( 2.7, pow(pixelvalue ,2) );

  double c = .99 - abs(127-pixelvalue)/127;

  if (c < 0) c = 0;

  return c;

}

double light_f(double x)
{

  //return x*x;

  return pow(x, 1/EXPONENT);

}

double light_q(double x)
{
  //return sqrt(x);

  return pow(x, EXPONENT);
 
}


void
usage(void)
{

  fprintf(stderr, "\n\nlight vector merge utility\n");
  fprintf(stderr, "----------------------------------------------------------------\n\n");

  fprintf(stderr, " Use: pnmpwadd -f input-file -o imageout\n");
  fprintf(stderr, "  or: pnmpwadd image1 image2 -o imageout\n\n");
  fprintf(stderr, " the input file has to be a text file, where each line\n"
                  " is as follows:\n\n");
  fprintf(stderr, " <file name>.ppm <red value> <green value> <blue value>\n\n");
  fprintf(stderr, " it should specify names of image files and their weighting\n"
                  " for the composition of the final image.\n"
                  " the weights are per colour\n\n");

  exit(EXIT_SUCCESS);
}

void
parse_commandline(int argc, char ** argv, struct image_params * a_params,
		  struct image_params * b_params, char ** output_filename,
                  char ** input_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 (!strcasecmp(argv[i], "-f"))
        {
          if (argc<=(i+1)) usage();
          *input_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)) &&
     (input_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 );
  }

}


// data from the image vectors list:
struct image_params img_params[MAX_NUM_OF_IMAGES];
FILE *image_file[MAX_NUM_OF_IMAGES];
unsigned char image_pixel[3][MAX_NUM_OF_IMAGES];
int bytes_per_pixel, num_of_images;


void init_img_params()
{
  int i;

  for (i=0; i<MAX_NUM_OF_IMAGES; i++)
  {
   img_params[i].filename = NULL;
   img_params[i].width = 0;
   img_params[i].height = 0;
   img_params[i].max_val = 0;
   img_params[i].type = 0;
   img_params[i].pix_color_weight[0] = 1.0; 
   img_params[i].pix_color_weight[1] = 1.0;
   img_params[i].pix_color_weight[2] = 1.0;

  }
}

void open_img_files()
{
  int i;
 
  printf("opening image files ...\n");

  for (i=0; i<MAX_NUM_OF_IMAGES; i++)
  {
    if ( !( (img_params[i].filename == NULL) || (strlen(img_params[i].filename) == 0) ) )
    {

      if (( image_file[i] = fopen(img_params[i].filename, "r"))==NULL)
      {
        fprintf(stderr, "Unable to open %s.\n", img_params[i].filename);
        exit(EXIT_FAILURE);
      }

      get_image_params(image_file[i], &img_params[i]);

      // check files to be of same size and format:

      if (i > 0)
      {
        if ((img_params[i].width!=img_params[i-1].width) || 
          (img_params[i].height!=img_params[i-1].height) ||
          (img_params[i].max_val!=img_params[i-1].max_val))
        {
          fprintf(stderr, "Images must have the same dimensions and maximum "
                 "colour values.\n");
          exit(EXIT_FAILURE);
        }

        if (img_params[i].type!=img_params[i-1].type)
        {
         fprintf(stderr, "%s is a P%c but %s is a P%c.\n", img_params[i].filename,
              img_params[i].type, img_params[i-1].filename, img_params[i].type);
         exit(EXIT_FAILURE);
        }
      }
      else
      {
        if (img_params[i].type=='5') bytes_per_pixel=1;
        else bytes_per_pixel=3;

      }

    }
  }

}

int
main(int argc, char ** argv)
{

  // data for image pair:
  struct image_params a_params={NULL,0,0,0,0};
  struct image_params b_params={NULL,0,0,0,0};

  char * dest_filename=NULL;
  char * input_filename=NULL;
  FILE * a, * b, * dest, * in_file;
  double p=EXPONENT;  /* exponent for pow */
  double p_inv;  /* 1.0/p */
  double scale_factor; /* to make the result on same interval is input */
  unsigned char a_pixel[3], b_pixel[3];
  unsigned char same_pixel[3]= {0,0,0};
  unsigned char out_pixel[3];
  int i, img;
  int done = 0;
  double weight_sum, pixel_sum;

  int action = 1;

  int img_num = 0;

  char image_file_name[50];
  float pix_color_weight[2]; // rgb weight for a pixel


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

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


  printf("\n\nlight vector merge utility\n");
  printf("------------------------------------------------------------------------------");

  printf("\n\n>> input file is \"%s\"\n>> output file is \"%s\"\n\n", input_filename, dest_filename);


  printf("image files:\n\n");

  if (input_filename)
  {

    init_img_params();

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

    img_num = 0;

    while (!done)
    {
      img_params[img_num].filename = (char *)malloc(50);

      fscanf(in_file, "%s %f %f %f\n", img_params[img_num].filename, 
                                       &(img_params[img_num].pix_color_weight[0]),
                                       &(img_params[img_num].pix_color_weight[1]),
                                       &(img_params[img_num].pix_color_weight[2]) ); 

      if ( (img_params[img_num].filename == NULL) || (strlen(img_params[img_num].filename) == 0) )
      {
        done = 1;
      }
      else
      {
        printf(">> file: \"%s\";  colour weight: r = %f, g = %f, b = %f \n", img_params[img_num].filename, 
                                       (img_params[img_num].pix_color_weight[0]),
                                       (img_params[img_num].pix_color_weight[1]),
                                       (img_params[img_num].pix_color_weight[2]) );
        img_num++;
      }

    }

  }

  num_of_images = img_num;

  printf("\n");

  open_img_files();
  
/*
  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]=img_params[0].max_val;


  // write file header:
  fprintf(dest, "P%c\n%d %d\n%d\n", img_params[0].type, img_params[0].width,
	  img_params[0].height, img_params[0].max_val);


  printf("merging %d images ...\n", num_of_images);


  while (!feof(image_file[0]))
    {
    printf("merging file %d...\n", image_file[0]);

      for (i=0; i<num_of_images; i++)
        fread(image_pixel[i], sizeof(unsigned char), bytes_per_pixel, image_file[i]);

//      if (memcmp(image_pixel[0], image_pixel[1], sizeof(unsigned char)*bytes_per_pixel))
	{
	    for (i=0;i<bytes_per_pixel;i++)
            {

              // weighting with light quantity modeling:
              if (action == 1)
              {

               pixel_sum = 0; weight_sum = 0;
               for (img=0; img<num_of_images; img++)
               {
                 pixel_sum += pow((double)image_pixel[img][i],p)*img_params[img].pix_color_weight[i];
                 weight_sum += img_params[img].pix_color_weight[i];
                 
               }

               out_pixel[i]=(unsigned char) ( pow(pixel_sum / weight_sum, p_inv)/scale_factor); 

/*
	       out_pixel[i]=(unsigned char)
                (
                  pow(
                      (
                       pow((double)image_pixel[0][i],p)*img_params[0].pix_color_weight[i]
                       +
                       pow((double)image_pixel[1][i],p)*img_params[1].pix_color_weight[i]

                      )/(img_params[0].pix_color_weight[i] + img_params[1].pix_color_weight[i])
                  ,p_inv)/scale_factor
                );
*/

              }

              // weighting with no light quantity modeling:
              else if (action == 2)
               out_pixel[i]=(unsigned char)
                (
                    
                      (
                       ((double)a_pixel[i])*WA
                       +
                       ((double)b_pixel[i])*WB
                      )/(WA+WB)
                );

              // weighting with no light quantity modeling but with certainty func:
              else if (action == 3)
               out_pixel[i]=(unsigned char)
                (
                      (
                       (double)a_pixel[i]*certainty_function((double)a_pixel[i])*WA
                       +
                       (double)b_pixel[i]*(certainty_function((double)b_pixel[i]))*WB
                      )/(WA*certainty_function((double)a_pixel[i]) + WB*certainty_function((double)b_pixel[i]))
                );

              // weighting with light quantity modeling and certainty func:
              else if (action == 4)
               out_pixel[i]=(unsigned char)
                (
                  light_f(
                      (
                       light_q((double)a_pixel[i])*certainty_function((double)a_pixel[i])*WA
                       +
                       light_q((double)b_pixel[i])*certainty_function((double)b_pixel[i])*WB
                      )/
                      (WA*certainty_function((double)a_pixel[i]) + 
                       WB*certainty_function((double)b_pixel[i])
                      )
                  )
                );

            }
	    fwrite(out_pixel, sizeof(unsigned char), bytes_per_pixel, dest);
	}

      //else
//	fwrite(same_pixel, sizeof(unsigned char), bytes_per_pixel, dest);
    }


  // close files:
  for (i=0; i<num_of_images; i++)
    fclose(image_file[i]); 
  fclose(dest);

  //printf("c(0) = %f, c(127) = %f, c(256) = %f\n", certainty_function(0), certainty_function(127), certainty_function(256));

  printf("\n\ndone.\n\n");

  exit(EXIT_SUCCESS);
}