/* wyckoff principle */

/* rough program that will serve as a skeletal outline
   for various aspects of the wyckoff principle */

#define K 2
#define k 4
#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]={0,0,0};
  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;
  /* printf("%d %d  %d\n", dims[0], dims[1], dims[2]); */
  do {
    char line[80];
    fgets( line, sizeof(line), ifs );
    /* printf("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" );
      }
    }
    /* printf("dims_scanned:%d\n", dims_scanned); */
  } while( dims_scanned < 3 );
  /* printf("%d %d  %d\n", dims[0], dims[1], dims[2]); */

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

  /* printf("%d %d  %d\n", dims[0], dims[1], dims[2]); */
  if( params->max_val != 255 ){
    fprintf(stderr,"get_frame_params: max value must be 255\n");
    exit( EXIT_FAILURE );
  }

}

float
certainty(double point){

 float x;
 x=6*(point/255)-3;
 return(exp(-(x*x)));
}

double
cement(double a, double b){
/* f was bad choice function name because hard to search for */

 double c, x, result;
 /* result=(k*a+b)/(k+1); */
 result=(k*certainty(a)*a+certainty(b)*b)/(k*certainty(a)+1*certainty(b));
 /* result=(k*certainty(a)*sqrt(a)+certainty(b)*sqrt(b))/(k*certainty(a)+certainty(b));*/
 if (result>255.0) result=255.0;
 /* printf("%f\n",result); */
 return(result);
}
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=2.7;  /* 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 bytes_per_pixel, i;
 
  /* printf("certainty is %f",certainty(0)); */
 
  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);
        for (i=0;i<bytes_per_pixel;i++)
          out_pixel[i]=(unsigned char)
            (
              /* cement((double)a_pixel[i],(double)b_pixel[i])/scale_factor */
              cement((double)a_pixel[i],(double)b_pixel[i])
              /* make sure you scale output on 0..255*/
            );
      fwrite(out_pixel, sizeof(unsigned char), bytes_per_pixel, dest);
    }

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