#define K 2
#define k 4

/*

Take the Wyckoff challenge!

This is a fun way to learn about PhotoRendering (Photoquantigraphic Rendering).

obtain the two "mannfamily" test images
called "Dark exposure"
from http://wearcam.org/mannfamily_dark.pgm
and "Light exposure"
from http://wearcam.org/mannfamily_light.pgm
and combine the two images together
to render a single image of improved dynamic range and tonal fidelity, as
described in http://wearcam.org/wyckoff.html (don't use the gif files
on the www page, use the pgm files, because the C-code below doesn't work
with gif files).

Experiment with different ways of combining the images.

Try first multiplying the first by four (to lighten it up enough to
"register" it with the second) and then adding the two, e.g. 4a+b
where a is the first and b is the second.  Don't forget to scale the
result back down or you will overrun the 0..255 image interval.  When
you scale down, what you're doing is a weighted average:
(4a+1b)/(4+1), where the weights are 4 and 1.

Then try using weights that depend on the image pixel value.  The
weighting should be higher for pixel values near 127, and lower for
pixel values near 0 or 255.  this can be done by defining a gaussian
certainty function, c, such that c=exp(-x^2), or c=pow(e,pow(x,2)) in "C",
where e is the naperian base (2.7....) and x=6*(pixelvalue/255)-3
makes it roll off nicely at "3 sigma" either side (e.g. so pixel value
0 has certainty very small as does pixel value 255, but pixel value 127
has a high certaintly value).  The reason for doing this is so that extreme
pixel values (which are under or over exposed) are weighted less while
properly or moderate exposed values more.  Implement the "Certainty"
function as a function call.  Test it to see if you get
certainty(0)=1.2341e-04, and certainty(127)=0.9999

So now try calculating the weighted sum: (4*c(a)*a + c(b)*b)/(4*c(a)+c(b))

Finally, if you feel adventurous, then try computing the quantity
f((4*c(a)*finv(a) + c(b)*finv(b))/(4*c(a)+c(b)))
where f is a nonlinear but monotonic function of positive 2nd derivative.
A good one to try is f(x)=x*x (the "squaring" function),
in which case finv()=sqrt().

Notice how the nonlinearity improves the photorendering.

--Prof. Steve Mann, 1997

*/

/*
 *derived from pnmbmath.c, steve, for wyckoff principle
 *to compile: gcc -Wall pnmwyckoff.c -o pnmwyckoff -lm
 *to run: pnmwyckoff darkpicture.ppm lightpicture.ppm -o test.ppm
 *
 *(from ndf's diff.c, with help of sbeck@media.mit.edu on parsing input images)
 *ndf, steve, sbeck
 *mathematical theory and original concept by steve, 1992
 */

#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 [-mark] 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 * mark)
{
  int i;

  if (argc<3) usage();

  for (i=1;i<argc;i++)
    {
      if (!strcasecmp(argv[i], "-mark"))
	*mark=1;
      else 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 );
    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 );

  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 */
  unsigned char a_pixel[3], b_pixel[3];
  unsigned char same_pixel[3]= {0,0,0};
  unsigned char out_pixel[3];
  int mark=0;
  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, &mark);

  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))
        {
	  if (!mark)
	    for (i=0;i<bytes_per_pixel;i++)
	      out_pixel[i]=(unsigned char)
                (
                  f(                /* if is the nonlinear function */
                  /* assume we know first image is 1/4 the exposure of second
                     so it needs to be increased by k=4 times to "register"
                     it with the second (lighter) image */
                     k*certainty((double)a_pixel[i])*finv((double)a_pixel[i])
                     +
                       certainty((double)b_pixel[i])*finv((double)b_pixel[i])
                   )/scale_factor  /* make sure you scale output on 0..255/*
                );
	  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);
}