/*
 *  comparagram.c
 *  
 compile with:
gcc -Wall -O2 -g  -o comparagram comparagram.c PNMImageOffsetable.c -ljpeg -lpng -lm
 *  read in 2 images, output a comparagram image. 
 *  comparagram [-o out.pnm] in.pnm in.pgm in.jpg
 *
 *  here alpha means "certainty", with zero implying
 *  no knowledge of a pixel's actual value (this can
 *  be approximated by saying an alpha-zero pixel is
 *  transparent); for byte wide channels, this value
 *  goes up to 255 for absolutely certain pixels.
 * 
 *  there is a special case when reading in ".pnm"
 *  files with the special comments of the form:
 
# MOffset 0
# NOffset 0
 
 * in which case we assume that any grey channel set
 * to 255 is actually NaN, and likewise for any rgb
 * triplet (255,255,255). this is the "myNaN" convention.
 * this promotes the file from "type" to "type_ALPHA".
 * 
 * maali@eyetap.org 
 * dec 9, 2000
 *
 */


#include <stdio.h>
#include <stdlib.h>

#include <jpeglib.h>
#include <memory.h>
#include <limits.h>
#include <string.h>
#include <math.h>
#include "PNMImageOffsetable.h"
#include "simple_png_read.c"
#include <setjmp.h> // last

#ifdef MAX
#undef MAX
#endif 
#define MAX(a, b)  (((a) > (b)) ? (a) : (b))

#ifdef MIN
#undef MIN
#endif 
#define MIN(a, b)  (((a) < (b)) ? (a) : (b))

//#define REACHED fprintf(stderr, "Reached %s:%d\n",__FILE__,__LINE__); 
#define REACHED

enum image_type {
  IMAGE_TYPE_NONE,
  IMAGE_BYTE_GREY,
  IMAGE_BYTE_GREY_ALPHA,
  IMAGE_BYTE_RGB,
  IMAGE_BYTE_RGB_ALPHA,
  IMAGE_DOUBLE_GREY,
  IMAGE_DOUBLE_GREY_ALPHA,
  IMAGE_DOUBLE_RGB,
  IMAGE_DOUBLE_RGB_ALPHA,
};

enum image_filetype {
  IMAGE_FILE_NONE,
  IMAGE_FILE_PNM, // IMAGE_FILE_PNM_NAN,
  IMAGE_FILE_JPG,
  IMAGE_FILE_PNG
};

struct image {
 enum image_type type;
 int x_size;   int y_size;
 int n_channels;        // 1:g 2:ga 3:rgb 4:rgba
 int bytes_per_channel;
 int x_offset; int y_offset;
 int n_color_channels;  // 1:g,ga 3:rgb,rgba 4:cmyk,cmyka
 void *data;  
};


void image_init(struct image *i)
{
  memset(i, 0, sizeof(struct image));
  i->data=NULL;
}

void image_destroy(struct image *i)
{
  if ((i!=NULL)&&(i->data!=NULL)) free(i->data);
}

int image_get_data_size(struct image *i)
{
  if (i!=NULL) {
    return (i->x_size * i->y_size * i->n_channels * i->bytes_per_channel);
  }
  return 0;
}

void image_copy(struct image *s, struct image *d)
{
  memmove(d,s,sizeof(struct image));
}

void image_deep_copy(struct image *s, struct image *d)
{
  int datasize=0;
  if((s!=d)&&(s!=NULL)&&(d!=NULL)) {
    memmove(d,s,sizeof(struct image));
    datasize=image_get_data_size(s);
    d->data = malloc(datasize);
    memmove(d->data,s->data,datasize);
  } else {
    fprintf(stderr, "image_deep_copy: null or identical pointer arguments.\n");
  }
}

// return true or false
int image_is_pnm_filename(char *fn) 
{
  char *fn_ext=fn+strlen(fn)-4;
  if ((fn_ext[0]=='.') && ((fn_ext[1]=='p')||(fn_ext[1]=='P'))
                       && ((fn_ext[2]=='n')||(fn_ext[2]=='N')||
                           (fn_ext[2]=='g')||(fn_ext[2]=='G')||
                           (fn_ext[2]=='p')||(fn_ext[2]=='G'))
                       && ((fn_ext[3]=='m')||(fn_ext[3]=='M')))
  { return 1; }
  
  return 0;
}


// return true or false
int image_is_png_filename(char *fn) 
{
  char *fn_ext=fn+strlen(fn)-4;
  if ((fn_ext[0]=='.') && ((fn_ext[1]=='p')||(fn_ext[1]=='P'))
                       && ((fn_ext[2]=='n')||(fn_ext[2]=='N'))
                       && ((fn_ext[3]=='g')||(fn_ext[3]=='G')))
  { return 1; }
  
  return 0;
}

// return true or false
int image_is_jpg_filename(char *fn)
{
  char *fn_ext=fn+strlen(fn)-4;
  if ((fn_ext[0]=='.') && ((fn_ext[1]=='j')||(fn_ext[1]=='J'))
                       && ((fn_ext[2]=='p')||(fn_ext[2]=='P'))
                       && ((fn_ext[3]=='g')||(fn_ext[3]=='G')))
  { return 1; }

  fn_ext=fn+strlen(fn)-5;
  if ((fn_ext[0]=='.') && ((fn_ext[1]=='j')||(fn_ext[1]=='J'))
                       && ((fn_ext[2]=='p')||(fn_ext[2]=='P'))
                       && ((fn_ext[3]=='e')||(fn_ext[3]=='E'))
                       && ((fn_ext[4]=='g')||(fn_ext[4]=='G')))
  { return 1; }
  
  return 0;
}

struct my_error_mgr {
  struct jpeg_error_mgr pub;    /* "public" fields */

  jmp_buf setjmp_buffer;        /* for return to caller */
};

typedef struct my_error_mgr * my_error_ptr;

METHODDEF(void)
my_error_exit (j_common_ptr cinfo)
{
  /* cinfo->err really points to a my_error_mgr struct, so coerce pointer */
  my_error_ptr myerr = (my_error_ptr) cinfo->err;

  /* Always display the message. */
  /* We could postpone this until after returning, if we chose. */
  (*cinfo->err->output_message) (cinfo);

  /* Return control to the setjmp point */
  longjmp(myerr->setjmp_buffer, 1);
}

int image_load_jpg(char *fn, struct image *i)
{
  FILE * infile;
  struct jpeg_decompress_struct cinfo;
  struct my_error_mgr jerr;
  JSAMPARRAY buffer;		/* Output row buffer */
  int row_stride;		/* physical row width in output buffer */
  unsigned char *dest;


  if ((infile = fopen(fn, "rb"))==NULL) {
    fprintf(stderr, "Couldn't open input file '%s'.\n", fn);
    return -1;
  }

  cinfo.err = jpeg_std_error(&jerr.pub);
  jerr.pub.error_exit = my_error_exit;
  if (setjmp(jerr.setjmp_buffer)) {
    jpeg_destroy_decompress(&cinfo);
    return 0;
  }

  jpeg_create_decompress(&cinfo);
  jpeg_stdio_src(&cinfo, infile);
  (void) jpeg_read_header(&cinfo, TRUE);
  (void) jpeg_start_decompress(&cinfo);
  row_stride = cinfo.output_width * cinfo.output_components;
  

  i->x_size     = cinfo.output_width;
  i->y_size     = cinfo.output_height;
  i->x_offset   = 0;
  i->y_offset   = 0;
  i->n_channels = 
    i->n_color_channels = cinfo.out_color_components;
  switch (cinfo.out_color_components) {
    case 1: i->type = IMAGE_BYTE_GREY; break;
    case 3: i->type = IMAGE_BYTE_RGB; break;
    default: i->type = IMAGE_TYPE_NONE; break;
  } 
  i->bytes_per_channel = 1;
  i->data = malloc(cinfo.output_height * row_stride);
  if (i->data == NULL) {
    fprintf(stderr, "image_load_jpg: malloc failed.\n");
    return -1;
  }
  
  buffer=(*cinfo.mem->alloc_sarray)((j_common_ptr) &cinfo, JPOOL_IMAGE, row_stride, 1);
  dest = (unsigned char*)i->data;
  //fprintf(stderr,"size %dx%d, %d channels, rowstride=%d\n", 
  //    i->x_size, i->y_size, i->n_channels, row_stride);
  while (cinfo.output_scanline < cinfo.output_height) {
    jpeg_read_scanlines(&cinfo, buffer, 1);
    memcpy(dest, buffer[0], row_stride);
    dest+=row_stride;
  }

  jpeg_finish_decompress(&cinfo);
  jpeg_destroy_decompress(&cinfo);
  fclose(infile);

  return 0;
}

#define mynan 255
int image_load_pnm(char *fn, struct image *i)
{
  FILE *f;
  char buf[41];
  int bufsize=40, j;
  int use_my_nan=0;
  unsigned char *d, r,g,b;
  double *s;
  PNMImageOffsetable pnm;

  printf("image_load_pnm: reading image file '%s'.\n", fn);
  PNM_InitFromFile(&pnm, fn);
  i->x_size = pnm.xSize;
  i->y_size = pnm.ySize;
  i->n_channels = pnm.depth;
  i->n_color_channels = pnm.depth;
  i->x_offset = pnm.nOffset;
  i->y_offset = pnm.mOffset;

  // test for mynan stuff
  f=fopen(fn,"rb");
  fread(buf, 1, bufsize, f);
  fclose(f);
  buf[40]=0;

  if (strstr(buf, "MOffset")!=NULL) {
    // use myNan convention
    use_my_nan=1;
    printf("image_load_pnm: using myNaN convention.\n");
  } else {
    printf("image_load_pnm: not using myNaN convention.\n");
  }

  switch (i->n_channels+use_my_nan) {
    case 1: i->type = IMAGE_BYTE_GREY; break;
    case 2: i->type = IMAGE_BYTE_GREY_ALPHA; break;
    case 3: i->type = IMAGE_BYTE_RGB; break;
    case 4: i->type = IMAGE_BYTE_RGB_ALPHA; break;
    default: i->type = IMAGE_TYPE_NONE; break;
  } 

  if (!use_my_nan) {
    // g and rgb
    i->data = malloc(pnm.xSize*pnm.ySize*pnm.depth);
    // type mismatch for memcpy
    //memcpy(i->data, pnm.data, pnm.xSize*pnm.ySize*pnm.depth);
      s=(double*)pnm.data;
      d=(unsigned char *)i->data;
      for(j=0; j<pnm.xSize*pnm.ySize*pnm.depth; j++) {
        *d++ = (unsigned char) *s++;
      }
  } else {
    i->n_channels++;
    i->data = malloc(pnm.xSize*pnm.ySize*(i->n_channels));

    if (pnm.depth==1) {
      // gs+mynan -> ga
      s=(double*)pnm.data; 
      d=(unsigned char *)i->data;
      for(j=0; j<pnm.xSize*pnm.ySize; j++) { // per pixel
        if ((*d++=(unsigned char)*s++)==mynan) {
          *d++ = 0;
        } else {
          *d++ = 255;
        }
      }
    } else if (pnm.depth==3) {
      // rgb+mynan -> rgba
      fprintf(stderr,"image_load_pnm: depth=3, converting to rgb_alpha...\n");
      s=(double *)pnm.data; 
      d=(unsigned char *)i->data;
      for(j=0; j<pnm.xSize*pnm.ySize; j++) { // per pixel
        r=(unsigned char)*s++;
        g=(unsigned char)*s++;
        b=(unsigned char)*s++;
        if ((r==mynan)&&(g==mynan)&&(b==mynan))
        {
          *d++ = 0;
          *d++ = 0;
          *d++ = 0;
          *d++ = 0;
        } else {
          *d++ = r;
          *d++ = g;
          *d++ = b;
          *d++ = 255;
        }
      }
    } else {
      printf("image_load_pnm: unsupported pnm depth.\n");
      return -1;
    }
  }  
  return 0;
}

int image_load_png(char *fn, struct image *i)
{
  REACHED;
  simple_png_read(fn, 
    &(i->n_color_channels), 
    ((unsigned char**)&(i->data)), 
    &(i->x_size), &(i->y_size),
    &(i->x_offset), &(i->y_offset));
  i->bytes_per_channel = 1;

  // "simple" always adds alpha layer
  i->n_channels = i->n_color_channels+1; 
  //fprintf(stderr, "image_load_png: n_color_channels=%d\n", 
  //  i->n_color_channels);
  switch(i->n_color_channels) {
    case 1:
    case 2: i->type = IMAGE_BYTE_GREY_ALPHA; break;
    case 3:
    case 4: i->type = IMAGE_BYTE_RGB_ALPHA; break;
    default: i->type = IMAGE_TYPE_NONE;
      fprintf(stderr, "image_load_png: "
          "warning, unsupported png type loaded.\n");
  }
  REACHED;
  return 0;
}

// return zero on success
int image_read(char*filename, struct image *i)
{
  FILE *fp=NULL;
  unsigned char buf[8];
  int bufsize=sizeof(buf);
  int retval=0;
  enum image_filetype ift=IMAGE_FILE_NONE;

  // init
  image_init(i);
  memset(buf, 0, bufsize);
  fprintf(stderr, "image_read: trying to load '%s'\n", filename);

  // read file signature (magic number)
  if ((fp=fopen(filename, "rb"))==NULL) {
    fprintf(stderr, "image_read: Couldn't open file '%s'.\n", filename);
    exit(-1);
    return -1;
  }
  retval=fread(buf, 1, bufsize, fp);

  // try to detect filetype
  if ((buf[0]==0xff)&&(buf[1]==0xd8)) 
  { 
    fprintf(stderr, "image_read: filetype JPG (%s)\n", filename); 
    ift=IMAGE_FILE_JPG;
    
  }
  else if ((buf[0]==137)&&(buf[1]==80)&&(buf[2]==78)&&(buf[3]==71)
         &&(buf[4]==13 )&&(buf[5]==10)&&(buf[6]==26)&&(buf[7]==10))
  { 
    fprintf(stderr, "image_read: filetype PNG (%s)\n", filename); 
    ift=IMAGE_FILE_PNG;
  }
  else if ((buf[0]=='P')&&(buf[1]=='5')) 
  { 
    fprintf(stderr, "image_read: "
        "filetype PGM \"rawbits\" image data (%s)\n", filename); 
    ift=IMAGE_FILE_PNM; // need to check for myNaN
  }
  else if ((buf[0]=='P')&&(buf[1]=='6')) 
  { 
    fprintf(stderr, "image_read: "
        "filetype PPM \"rawbits\" image data (%s)\n", filename); 
    ift=IMAGE_FILE_PNM; // need to check for myNaN
  }
  else if ((buf[0]=='P')&&((buf[1]=='9')||(buf[1]=='A'))) 
  { 
    fprintf(stderr, "image_read: "
        "filetype PLM \"rawbits\" image data (%s)\n", filename); 
  }
  else {
    // detect type from file name.. this shouldn't be
    // here (either), it should only be for output...
    if (image_is_png_filename(filename)) {
      fprintf(stderr, "image_read: "
          "filetype PNG (%s)\n", filename); 
      ift=IMAGE_FILE_PNG;
    } else if (image_is_pnm_filename(filename)) {
      fprintf(stderr, "image_read: "
          "filetype PNM \"rawbits\" image data (%s)\n", filename);
      ift=IMAGE_FILE_PNM; // need to check for myNaN
    } else if (image_is_jpg_filename(filename)) {
      fprintf(stderr, "image_read: "
          "filetype JPG image data (%s)\n", filename);
      ift=IMAGE_FILE_JPG; 
    } else {
      fprintf(stderr, "image_read: "
          "Unknown filetype for file '%s'.\n", filename);
      return -1;
    }
  }

  // now 'ift' has the image filetype
  //rewind(fp);
  fclose(fp);
  switch (ift) {
    case IMAGE_FILE_PNM: image_load_pnm(filename, i); break;
    case IMAGE_FILE_JPG: image_load_jpg(filename, i); break;
    case IMAGE_FILE_PNG: image_load_png(filename, i); break;
    default: 
      fprintf(stderr, "image_read: Unknown filetype\n"); 
      return -1; 
    break;
  }
  
  fprintf(stderr, "image_read: "
      "size: %dx%d, channels: %d, offset: %dx%d, alpha: %s\n",
      i->x_size, i->y_size, i->n_color_channels, i->x_offset, i->y_offset,
      (i->n_channels-i->n_color_channels)?"yes":"no");

  return 0;
}

void swap_int(int*a, int*b) 
{
  int tmp=*a;
  *a=*b;
  *b=tmp;
}

// the convention here is that if we
// return i0>i1 this implies null set 
// for the intersection
void intersect_1d(int a0, int a1, int b0, int b1, int *i0, int *i1)
{
  fprintf(stderr, "intersect_1d: "
      "[%5d,%5d] intersect [%5d,%5d]", a0,a1,b0,b1);
  *i0 = *i1 = 0;
  if(a0>a1) { swap_int(&a0,&a1); }
  if(b0>b1) { swap_int(&b0,&b1); }

  /*  a:            a0|<------>|a1 
   *  b:                            b0|<------->|b1
   *  a intersect b : null
   *  
   *  a:                             a0|<------->|a1
   *  b:             b0|<------->|b1
   *  a intersect b : null
   */
  if((a1<b0)||(a0>b1)) { *i0=1; *i1=0; } 

  /*  a:              a0|<--------------->|a1 
   *  b:                           b0|<------->|b1
   *  a intersect b : {x:b0<x<a1}
   */
  else if((a0<=b0)&&(a1<=b1)) { *i0=b0; *i1=a1; }

  /*  a:             a0|<--------------->|a1 
   *  b:        b0|<------->|b1
   *  a intersect b : {x:a0<x<b1}
   */
  else if((a0>=b0)&&(a1>=b1)) { *i0=a0; *i1=b1; }

  /*  a:             a0|<--------------->|a1 
   *  b:        b0|<--------------------------->|b1
   *  a intersect b : {x:a0<x<a1}
   */
  else if((a0>=b0)&&(a1<=b1)) { *i0=a0; *i1=a1; }

  /*  a:        a0|<--------------------------->|a1 
   *  b:               b0|<------------->|b1
   *  a intersect b : {x:b0<x<b1}
   */
  else if((a0<=b0)&&(a1>=b1)) { *i0=b0; *i1=b1; }

  else { 
    fprintf(stderr,
      "intersect_1d: "
      "impossible condition (a0,a1,b0,b1)=(%d,%d,%d,%d)\n",
      a0,a1,b0,b1);
  }
  fprintf(stderr, "=[%5d,%5d]\n", *i0, *i1);

  return;
}

void 
comparagram_setup(
    struct image *a, struct image *b,
    int *astart, int *askip,
    int *bstart, int *bskip,
    int *icols0, int *icols1,
    int *irows0, int *irows1)
{
  *astart=0; 
  *bstart=0;
  *askip=0; 
  *bskip=0;

  intersect_1d(
    a->x_offset, a->x_size+a->x_offset, 
    b->x_offset, b->x_size+b->x_offset,
    icols0, icols1);

  intersect_1d(
    a->y_offset, a->y_size+a->y_offset, 
    b->y_offset, b->y_size+b->y_offset,
    irows0, irows1);

  // handle no overlap case
  if ((*irows0>*irows1)||(*icols0>*icols1)) {
    fprintf(stderr, 
      "comparagram_setup: no intersecting points in \n"
      "comparagram_setup: (%d,%d)-(%d,%d) and (%d,%d)-(%d,%d)\n",
      a->x_offset, a->y_offset, a->x_size+a->x_offset, a->y_size+a->y_offset,
      b->x_offset, b->y_offset, b->x_size+b->x_offset, b->y_size+b->y_offset);
    return;
  }

  // compute number of bytes to skip to start of image
  // and between rows of image data
  *astart  = (*irows0 - a->y_offset)*(a->x_size);
  *astart += (*icols0 - a->x_offset);
  *astart *= a->n_channels;

  *bstart  = (*irows0 - b->y_offset)*(b->x_size);
  *bstart += (*icols0 - b->x_offset);
  *bstart *= b->n_channels;

  *askip = ((a->x_size)-(*icols1-*icols0))*(a->n_channels);
  *bskip = ((b->x_size)-(*icols1-*icols0))*(b->n_channels);

  return;
}

void 
comparagram_g(unsigned int **cgram, struct image *a, struct image *b)
{
  int i,j;
  register unsigned char 
    *aptr, ag,
    *bptr, bg; 
  int irows0, irows1, icols0, icols1;
  int astart, bstart, askip, bskip;

  fprintf(stderr, "comparagram_g: begin\n");
  comparagram_setup(a, b,
      &astart, &askip, &bstart, &bskip,
      &icols0, &icols1, &irows0, &irows1);
  for (i=0; i<255; i++) {
    memset(cgram[i], 0, 256*sizeof(unsigned char));
  }

  // set aptr,bptr to start pixels
  aptr=(unsigned char*)a->data+astart;
  bptr=(unsigned char*)b->data+bstart;
  for (j=irows0; j<irows1; j++)   // per row of a intersect b
  {
    for (i=icols0; i<icols1; i++) // per pixel of a intersect b
    {
      ag=*aptr++; bg=*bptr++; // get values
      (cgram[ag][bg])++;      // count value
    }
    // set aptr,bptr to start of next row
    aptr+=askip;
    bptr+=bskip;
  }
}

void 
comparagram_rgb(unsigned int **cgram_r, unsigned int **cgram_g,
                unsigned int **cgram_b, struct image *a, struct image *b)
{
  int i,j;
  register unsigned char 
    *aptr, ar, ag, ab,
    *bptr, br, bg, bb; 
  int irows0, irows1, icols0, icols1;
  int astart, bstart, askip, bskip;

  fprintf(stderr, "comparagram_rgb: begin\n");

  comparagram_setup( a, b, 
      &astart, &askip, &bstart, &bskip,
      &icols0, &icols1, &irows0, &irows1);
  for (i=0; i<255; i++) {
    memset(cgram_r[i], 0, 256*sizeof(unsigned char));
    memset(cgram_g[i], 0, 256*sizeof(unsigned char));
    memset(cgram_b[i], 0, 256*sizeof(unsigned char));
  }

  // set aptr,bptr to start pixels
  aptr=(unsigned char*)a->data+astart;
  bptr=(unsigned char*)b->data+bstart;
  for (j=irows0; j<irows1; j++)   // per row of a intersect b
  {
    fprintf(stderr,"doing row %d of input images\r",j);
    for (i=icols0; i<icols1; i++) // per pixel of a intersect b
    {
      ar=*aptr++; br=*bptr++; // get reds
      ag=*aptr++; bg=*bptr++; // get greens
      ab=*aptr++; bb=*bptr++; // get blue

      (cgram_r[ar][br])++;      // count reds
      (cgram_g[ag][bg])++;      // count greens
      (cgram_b[ab][bb])++;      // count blues
    }
    // set aptr,bptr to start of next row
    aptr+=askip;
    bptr+=bskip;
  }
}

void 
comparagram_ga(unsigned int **cgram, struct image *a, struct image *b)
{
  int i,j;
  register unsigned char 
    *aptr, aa, ag, 
    *bptr, ba, bg;
  int irows0, irows1, icols0, icols1;
  int astart, bstart, askip, bskip;


  fprintf(stderr, "comparagram_ga: begin\n");

  comparagram_setup( a, b,
      &astart, &askip, &bstart, &bskip,
      &icols0, &icols1, &irows0, &irows1);
  for (i=0; i<255; i++) {
    memset(cgram[i], 0, 256*sizeof(unsigned char));
  }

  // set aptr,bptr to start pixels
  aptr=(unsigned char*)a->data+astart;
  bptr=(unsigned char*)b->data+bstart;
  for (j=irows0; j<irows1; j++)   // per row of a intersect b
  {
    for (i=icols0; i<icols1; i++) // per pixel of a intersect b
    {
      aa=*aptr++; ba=*bptr++; // get alphas
      ag=*aptr++; bg=*bptr++; // get values

      if (aa==0||ba==0) { // zero alpha
      } else {
        (cgram[ag][bg])++; // count value
      } 
    }
    // set aptr,bptr to start of next row
    aptr+=askip;
    bptr+=bskip;
  }
}


void 
comparagram_rgba(unsigned int **cgram_r, unsigned int **cgram_g,
                 unsigned int **cgram_b, struct image *a, struct image *b)
{
  int i,j;
  register unsigned char 
    *aptr, ar, ag, ab, aa, 
    *bptr, br, bg, bb, ba; 
  int irows0, irows1, icols0, icols1;
  int astart, bstart, askip, bskip;

  fprintf(stderr, "comparagram_rgba: begin\n");

  comparagram_setup(a, b,
      &astart, &askip, &bstart, &bskip,
      &icols0, &icols1, &irows0, &irows1);
  for (i=0; i<255; i++) {
    memset(cgram_r[i], 0, 256*sizeof(unsigned char));
    memset(cgram_g[i], 0, 256*sizeof(unsigned char));
    memset(cgram_b[i], 0, 256*sizeof(unsigned char));
  }

  // set aptr,bptr to start pixels
  aptr=(unsigned char*)a->data+astart;
  bptr=(unsigned char*)b->data+bstart;
  for (j=irows0; j<irows1; j++)   // per row of a intersect b
  {
    for (i=icols0; i<icols1; i++) // per pixel of a intersect b
    {
      ar=*aptr++; br=*bptr++; // get reds
      ag=*aptr++; bg=*bptr++; // get greens
      ab=*aptr++; bb=*bptr++; // get blue
      aa=*aptr++; ba=*bptr++; // get alphas

      if (aa==0||ba==0) { 
        // zero alpha
      } else {
        (cgram_r[ar][br])++; // count reds
        (cgram_g[ag][bg])++; // count greens
        (cgram_b[ab][bb])++; // count blues
      } 
    }
    // set aptr,bptr to start of next row
    aptr+=askip;
    bptr+=bskip;
  }
}

int main(int argc, char *argv[])
{
  FILE *out, *raw;
  //char *outfile_template = "out%04d.pnm";
  //char *rawfile_template = "out%04d.u32.raw";
  char *outfile_template = "comparagram.pnm";
  char *rawfile_template = "comparagram.u32.rgb.raw";

  int noutfile=0;
  struct image *i0,  *i1, i0actual, i1actual;
  int i,m,n, l, img;
  unsigned int **cgram, *cgram_data, *cgram_data_rgb,
               cgram_size=256, cmax, cmax_r, cmax_g, cmax_b;
  unsigned int L=3, M=256, N=256;
  unsigned int **cgram_r, **cgram_g, **cgram_b, ***cgram_rgb;
  unsigned char *cgram_byte;
  double cgram_double;
  long int csum;
  double dcmax;
  double thresh = .002;  /* threshold for log scale */
  char outfile[256], rawfile[256]; 

  REACHED;
  i0=&i0actual;
  i1=&i1actual;
  
  if ((argc<3)||(argc>3)) { 
    printf("usage: comparagram f1.jpg f2.jpg\n" 
    "This program currently accepts jpeg, png and pnm images.\n"
    "Just use two input images, and the output is in '%s'.\n"
    //"and a tonally adjusted rendition of the second file (in the\n"
    //"tonal range of the first image) will be stored in 'out.pnm'.\n"
    "Offsets and alpha values of zero are handled correctly.\n"
    "Supported types (all channels are 'unsigned char'):\n"
    "  (certainty is stored in the alpha channel, 0=no knowledge,\n"
    "   and the c channel is always added to png's; it is also added\n"
    "   to pnm inputs if \"MOffset\" is found near the top of the file)\n"
    "  RGB  -- from jpeg, pnm\n"
    "  G    -- from jpeg(?), pnm\n"
    "  RGBC -- from png (c channel always added), pnm (if with myNaN)\n"
    "  GC   -- from png (''), pnm ('')\n",
    outfile_template
    );
    return 1; 
  }

  // buffer for comparagram: one row of pointers to a matrix common block
  // let's be pigs and mallocate amount for color even if greyscale:
  REACHED;
  cgram = (unsigned int**) malloc(cgram_size*sizeof(unsigned int *));
  cgram_r = (unsigned int**)malloc(cgram_size*sizeof(unsigned int *));
  cgram_g = (unsigned int**)malloc(cgram_size*sizeof(unsigned int *));
  cgram_b = (unsigned int**)malloc(cgram_size*sizeof(unsigned int *));
  cgram_rgb = (unsigned int***) malloc(3*cgram_size*sizeof(unsigned int *));
  cgram_data = (unsigned int*)
    malloc(cgram_size*cgram_size*sizeof(unsigned int));
  cgram_data_rgb=(unsigned int*)
    malloc(3*cgram_size*cgram_size*sizeof(unsigned int));
  cgram_byte = (unsigned char*)
    malloc(3*cgram_size*cgram_size*sizeof(unsigned char));

  REACHED;
  for (i=0;i<cgram_size; i++) { cgram[i] = (int*)(cgram_data+i*cgram_size); } 

  REACHED;
  for (i=0;i<cgram_size; i++) { 
    cgram_r[i] = (int*)(cgram_data_rgb+i*cgram_size); 
    cgram_g[i] = (int*)(cgram_data_rgb+i*cgram_size + cgram_size*cgram_size); 
    cgram_b[i] = (int*)(cgram_data_rgb+i*cgram_size + 2*cgram_size*cgram_size); 
  } 
  cgram_rgb[0] = cgram_r;  // common block cgram_rgb so we can
  cgram_rgb[1] = cgram_g;  // index it like cgram_rgb[l][m][n]
  cgram_rgb[2] = cgram_b;  // from individual color components

  // start reading images
  REACHED;
  image_read(argv[img=1], i0);

  do {
  REACHED;
    image_init(i1);
    image_read(argv[++img], i1);

    // clear comparagrams (we use 3x the need RAM if greyscale)
    memset(cgram_data, 0, cgram_size*cgram_size*sizeof(unsigned int));
    memset(cgram_data_rgb, 0, 3*cgram_size*cgram_size*sizeof(unsigned int));
    memset(cgram_byte, 0, 3*cgram_size*cgram_size*sizeof(unsigned char));
    if(i0->type!=i1->type) {
      fprintf(stderr, "comparagram: "
          "pixmap type mismatch on files '%s' and '%s'\n", 
          argv[img-1], argv[img]);
      return -1;
    }

    if((i0->x_size!=i1->x_size)||(i0->y_size!=i1->y_size)) {

      fprintf(stderr, "comparagram: "
          "image size mismatch on files '%s' and '%s', continuing...\n", 
	  argv[img-1], argv[img]);
      //return -1;
    }
   
  REACHED;
    switch(i0->type) { 
      case IMAGE_BYTE_GREY:       comparagram_g   (cgram, i0, i1); break;
      case IMAGE_BYTE_GREY_ALPHA: comparagram_ga  (cgram, i0, i1); break;
      case IMAGE_BYTE_RGB: 
        comparagram_rgb (cgram_r,cgram_g,cgram_b, i0, i1); break;
      case IMAGE_BYTE_RGB_ALPHA:  
        comparagram_rgba(cgram_r,cgram_g,cgram_b, i0, i1); break;
      default: 
        printf("comparagram: filetype not supported ('%s') %d/%d.\n", 
          argv[img], i0->type, img);
        break;
    }

    // get byte sized comparametric output
  REACHED;
    switch(i0->type) { 
    case IMAGE_BYTE_GREY:   
    case IMAGE_BYTE_GREY_ALPHA:
  REACHED;
    for (i=0, csum=0, cmax=0, dcmax=0; i<cgram_size*cgram_size; i++) { 
      cmax=MAX(cmax, cgram_data[i]);  //cgram is pointer to pointer
      dcmax=MAX(dcmax, log((double)cgram_data[i])); //cgram_data is contiguous
      csum+=cgram_data[i];
    }
    if (cmax==0) {
      fprintf(stderr, "comparagram: "
          "no data returned from comparagram function\n");
    }
    fprintf(stderr, "comparagram: "
        "csum: %ld, cmax: %u, dcmax: %f\n", csum, cmax, dcmax);
    for (i=0; i<cgram_size*cgram_size; i++) { 
      //cgram_byte[i] = 255-(int)(255.0*(cgram_data[i]/(double)cmax)) ; 
        cgram_double = (double) cgram_data[i]; // integer to double
        cgram_byte[i] = 
          (unsigned char) ((log(cgram_double+thresh)-log(thresh))/ 
                           (log((double)cmax+thresh)-log(thresh))*255);
    }
    break;
    case IMAGE_BYTE_RGB:
    case IMAGE_BYTE_RGB_ALPHA:
      REACHED;
      cmax_r=0;
      cmax_g=0;
      cmax_b=0;
      for (m=0; m<M; m++) { 
        for (n=0; n<N; n++) { 
           cmax_r=MAX(cmax_r, cgram_r[m][n]);  //cgram is pointer to pointer
           cmax_g=MAX(cmax_g, cgram_g[m][n]);  //cgram is pointer to pointer
           cmax_b=MAX(cmax_b, cgram_b[m][n]);  //cgram is pointer to pointer
         }
      }
      if ((cmax_r==0)&&(cmax_g==0)&&(cmax_b==0)) 
      {
        fprintf(stderr, "comparagram: "
          "no data returned from comparagram function\n");
      }

      for (l=0; l<L; l++) { 
        for (m=0; m<M; m++) { 
          for (n=0; n<N; n++) { 
            cgram_double = (double) cgram_rgb[l][m][n]; // integer to double
            cgram_byte[3*(m*256+n)+l] = (unsigned char)
            //cgram_byte[m][n][l] = ..... would be easier to understand
              ((log(cgram_double+thresh)-log(thresh))/
              (log((double)cmax+thresh)-log(thresh))*255);
          }
        }
      }
    REACHED;
    break;
    default: 
      printf("comparagram: "
          "filetype not supported ('%s') %d/%d.\n", 
          argv[img], i0->type, img);
    }// end switch

    // write output
    REACHED;
    noutfile++;
    sprintf(outfile, outfile_template, noutfile);
    fprintf(stderr, "comparagram: ===> writing output file '%s'\n", outfile);
    if((out=fopen(outfile, "wb"))==NULL) {
      fprintf(stderr, "comparagram: error opening output file '%s'\n", outfile);
      return -1;
    }
    switch(i0->type) { 
    case IMAGE_BYTE_GREY:   
    case IMAGE_BYTE_GREY_ALPHA:
      fprintf(out, "P5\n256 256\n255\n");
      fwrite(cgram_byte, 1, cgram_size*cgram_size, out);
    break;
    case IMAGE_BYTE_RGB_ALPHA:
    case IMAGE_BYTE_RGB:
      fprintf(out, "P6\n%d %d\n255\n",M,N);
      fwrite(cgram_byte, 3, cgram_size*cgram_size, out);
    break;
    default: printf("comparagram: "
               "filetype not supported ('%s') %d/%d.\n", 
               argv[img], i0->type, img);
  }
  fclose(out);

  // write raw file
  REACHED;
    sprintf(rawfile, rawfile_template, noutfile);
    if((raw=fopen(rawfile, "wb"))==NULL) {
      fprintf(stderr, "comparagram: "
          "error opening output file '%s'\n", 
          rawfile);
      return -1;
    }
    switch(i0->type) {
    case IMAGE_BYTE_GREY:
    case IMAGE_BYTE_GREY_ALPHA:
      fprintf(stderr, "comparagram: "
        "writing raw 256x256 (unsigned int) matrix\n"
        "comparagram: output file '%s'\n", 
        rawfile);
      fwrite(cgram_data, sizeof(int), cgram_size*cgram_size, raw);
      break;
    case IMAGE_BYTE_RGB:
    case IMAGE_BYTE_RGB_ALPHA:
      fprintf(stderr, "comparagram: "
        "writing raw 3x256x256 (unsigned int) matrix)\n" 
        "comparagram: output file '%s'\n",
        rawfile);
      fwrite(cgram_data_rgb, sizeof(int), 3*cgram_size*cgram_size, raw);
      break;
    default:
      fprintf(stderr, "comparagram: no raw file written.\n");
    } 
 
    fclose(raw);

/* this section is out-of-sync with the changes made
 * for color input... will not work right now 
 */
//#define DO_TONAL_ALIGNMENT
#if defined(DO_TONAL_ALIGNMENT)
    { // output tonally aligned input image
      int i, j, k;
      unsigned char *p=i1->data, pr, pg, pb, pa;
      unsigned char lut[256];
      //int l; unsigned char lut2[256];

      for(i=0;i<256;i++) { lut[i]=0; }
      for(i=0;i<256;i++) {       // ref. image 
        for(j=0,k=0;j<256;j++) { // cor. image
          if (k<cgram[j][i]) { 
            k=cgram[j][i];
            //if (j>lut[i]) //force semimonitonic
            lut[i]=j;
          }
	}
        if(i<255) lut[i+1]=lut[i];
      }
      printf("comparametric map (no conditioning):\n");
      for(i=0;i<256;i++) { 
        if(i&&!(i%16)) { printf("\n"); } 
        printf(" %3d",lut[i]); 
      }
      printf("\n");

      /*
      // make lut semimonotonic REALLY LAME!!
      if (lut[255]==0) lut[255]=255;
      if (lut[0]>lut[1]) lut[0]=lut[1];
      for(i=0;i<256;i++) { if (lut[i]!=0) { k=i; break; } } 
      for(i=1;i<256;i++) { if (lut[255-i]!=255) { l=255-i+1; break; } } 
      while (1)
      {
        j=1;
        for (i=0; i<254; i++) { 
          if (lut[i]-5>lut[i+1]) { 
            fprintf(stderr,"i=%3d: %d>>%d\n", i, lut[i], lut[i+1]);
            j=0; 
            break; } 
        }
        if (j) break;
        if(k==l) break;
        for (i=k; i<l; i++) { 
          if (lut[i]>lut[i+1]) 
            lut2[i]=i
              (unsigned char)(((int)lut[i-1]+2*(int)lut[i]+(int)lut[i+1])/4);
          else
            lut2[i]=lut[i];
        }
        for (i=k; i<l; i++) { lut[i]=lut2[i]; }
        lut2[0]=255;
      }
      
      if(lut2[0]==255){
        printf("comparametric map (after conditioning):\n");
        for(i=0;i<256;i++) { 
          if(i&&!(i%16)) { printf("\n"); } 
          printf(" %3d",lut[i]); 
        }
        printf("\n");
      }
    */

      out=fopen("out.pnm", "wb");
      fprintf(stderr, 
          "comparagram: output '%s' is '%s'\n"
          "comparagram: tonally aligned with '%s'\n",
          "out.pnm", argv[img], argv[img-1]);
      switch(i1->type)
      {
      case IMAGE_BYTE_GREY:
      fprintf(out, "P2\n%d %d\n255\n#\n", i1->x_size, i1->y_size);
        for (j=0; j<i1->y_size; j++)
        {
          for (i=0; i<i1->x_size; i++)
          {
            k++; fprintf(out, "%3d ", lut[*p++]);
            if (!(k%12)) { fprintf(out, "\n"); }
          }
        }
        break;
      case IMAGE_BYTE_GREY_ALPHA:
      for(i=0;i<255;i++) { if (lut[i]==255) { lut[i]=254; } }
      lut[255]=255;
      fprintf(out, "P5\n# MOffset %d\n# NOffset %d\n%d %d\n255\n", 
          i1->y_offset, i1->x_offset, i1->x_size, i1->y_size);
        for (j=0; j<i1->y_size; j++)
        {
          for (i=0; i<i1->x_size; i++)
          {
            pg = *p++;
            pa = *p++;
            if (pa==0) { pg = 255; }
            fwrite(lut+pg, 1, 1, out);
          }
        }
        break;
      case IMAGE_BYTE_RGB:
      fprintf(out, "P3\n%d %d\n255\n#\n", i1->x_size, i1->y_size);
        for (j=0; j<i1->y_size; j++)
        {
          for (i=0; i<i1->x_size; i++)
          {
            k++; fprintf(out, "%3d ", lut[*p++]);
            k++; fprintf(out, "%3d ", lut[*p++]);
            k++; fprintf(out, "%3d ", lut[*p++]);
            if (!(k% 3)) { fprintf(out, "  "); }
            if (!(k%12)) { fprintf(out, "\n"); }
          }
        }
        break;
      case IMAGE_BYTE_RGB_ALPHA:
      for(i=0;i<255;i++) { if (lut[i]==255) { lut[i]=254; } }
      lut[255]=255;
      fprintf(out, "P6\n# MOffset %d\n# NOffset %d\n%d %d\n255\n", 
          i1->y_offset, i1->x_offset, i1->x_size, i1->y_size);
        for (j=0; j<i1->y_size; j++)
        {
          for (i=0; i<i1->x_size; i++)
          {
            pr = *p++;
            pg = *p++;
            pb = *p++;
            pa = *p++;

            if (pa==0) { pr = pg = pb = 255; }
            fwrite(lut+pr, 1, 1, out);
            fwrite(lut+pg, 1, 1, out);
            fwrite(lut+pb, 1, 1, out);
          }
        }
        break;
      default:
      }
    }
    fclose(out);
    fprintf(stderr, 
      "cutpaste: xv \"%s\" \"%s\" \"%s\" \"%s\"\n",
      "out0001.pnm", argv[img], "out.pnm", argv[img-1]);
#endif

    // next pair
    image_destroy(i0);
    i0=i1;
  } while (img<argc-1);
  //image_destroy(i1);

  return 0;
}