/* This code provided as is */

#include <stdio.h>
#include <math.h>

/* Add #define RGBA if your are on a big-endian machine (ie non-Intel) */

#define RGBA
#define FILENAME_SUFFIX 8
#define PI              3.14159265359
#define MFALSE          0
#define MTRUE           !MFALSE

unsigned char iiiAverageSurroundingPixels(float u, float v,
                                          unsigned char *grayimage,
                                          int w, int h);
unsigned char *iiEncodeGrayImage(int w, int h, unsigned char *grayimage);
unsigned char *iiAverageGrayImages (int w, int h, unsigned char **grayimage,
                                    int numberofimages);
unsigned char *iiConvertToGrayscale(unsigned int *rgbimage, int w, int h);
unsigned char *iiReadPGMFile (FILE *fp, int *w, int *h);
unsigned int *iiReadPPMFile (FILE *fp, int *w, int *h);
int iiWritePPMFile (char *filename, unsigned int *pixel, int width, int height);
int iiWritePGMFileFloat (char *filename, float *pixel, int width, int height);
int iiWritePGMFileInt (char *filename, unsigned char *pixel, int width,
                       int height);
int iiWriteHistogramFile (char *filename, unsigned char *pixel, int width,
                          int height);
void iiComputeHistogram (unsigned char *pixel, int width, int height,
                         unsigned int *histogram, unsigned int histsize);
void iiComputeLookup(unsigned int *histogram, float *lookup, int width,
                     int height, unsigned int histsize);
void iiEqualizeImage(char *filename, unsigned char *grayimage, int width,
                     int height);
int iiGetRed(unsigned int pixVal);
int iiGetGreen(unsigned int pixVal);
int iiGetBlue(unsigned int pixVal);

#if 0
unsigned int *floatToColor(float *pixels, int w, int h)
{
  int r,g,b;
  int x,y;
  unsigned int *newPix;
  
  newPix = new unsigned int[w * h];

  for (y = 0; y < h; y++)
  {
    for (x = 0; x < w; x++)
    {
      r  = g = b = (int)(pixels[y * w + x]);
      newPix[y * w + x] = 0;
      newPix[y * w + x] |= (r & 0xFF);
      newPix[y * w + x] |= (g & 0xFF) << 8;
      newPix[y * w + x] |= (b & 0xFF) << 16;
      /*      printf("%f %x\n", pixels[y * w + x], newPix[y * w + x]); */

    }
  }

  return newPix;
}
#endif

/* -----------------------------------------------------------------------------
NAME:
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int iiGetRed(unsigned int pixVal)
{
  return 0xFF & pixVal;
}

/* -----------------------------------------------------------------------------
NAME:
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int iiGetGreen(unsigned int pixVal)
{
  return 0xFF & (pixVal >> 8);
}

/* -----------------------------------------------------------------------------
NAME:
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int iiGetBlue(unsigned int pixVal)
{
  return 0xFF & (pixVal >> 16);
}

/* -----------------------------------------------------------------------------
NAME:     iiReadPGMFile
PURPOSE:  Reads a PGM File into an unsigned char array
PRE:      The fp should be opened
POST:     The fp should be closed
RETURN:   An array of unsigned chars
AUTHOR:   CS Dept BYU
DATE:
MODIFIED: CMS-September 20, 1999
----------------------------------------------------------------------------- */
unsigned char *iiReadPGMFile (FILE *fp, int *w, int *h)
{
  int            max, type, x, y;
  char           buffer[250];
  unsigned char *pixel=NULL;

  if (!fp)
    return(pixel);
  fscanf(fp,"P%d\n",&type);
  if (type != 5)
    return NULL;
  buffer[0] = '#';
  while(buffer[0] == '#')
    fgets(buffer,250,fp);
  sscanf(buffer,"%d%d", w, h);
  pixel = (unsigned char *)malloc((*w) * (*h) * sizeof(unsigned char));
  if (!pixel)
    return(pixel);
  fgets(buffer,250,fp);
  sscanf(buffer,"%d",&max);
  for (y=0; y<(*h); y++) {
    for (x=0; x<(*w); x++) {
      fscanf(fp,"%c",&(pixel[y*(*w)+x]));
    }
  }
  return pixel;
} /* iiReadPGMFile */

/* -----------------------------------------------------------------------------
NAME:     iiReadPPMFile
PURPOSE:  Reads a PPM file into an unsigned int array
PRE:      The fp should be opened
POST:     The fp should be closed
RETURN:   An array of 
AUTHOR:   CS Dept BYU
DATE:
MODIFIED: CMS-September 10, 1999
----------------------------------------------------------------------------- */
unsigned int *iiReadPPMFile (FILE *fp, int *w, int *h)
{
  int           max, type;
  char          buffer[250];
  unsigned int *pixel = NULL;
  int           rint, gint, bint, x, y;
  unsigned char ruchar, guchar, buchar;

  if (!fp)
    return(pixel);

  fscanf(fp,"P%d\n",&type);
  if (type != 5 && type != 6 && type != 2 && type != 3)
    return NULL;
  buffer[0] = '#';
  while(buffer[0] == '#')
    fgets(buffer,250,fp);
  sscanf(buffer, "%d%d", w, h);
  pixel = (unsigned int *)malloc((*w) * (*h) * sizeof(unsigned int));
  if (!pixel)
    return(pixel);
  fgets(buffer,250,fp);
  sscanf(buffer,"%d",&max);
  for (y=0; y<(*h); y++) {
    for (x=0; x<(*w); x++) {
      if (type == 2 || type == 3) {
        if (type == 3)
          fscanf(fp,"%d%d%d",&rint,&gint,&bint);
        else {
          fscanf(fp,"%d",&rint);
          gint=bint=rint;
        }
        pixel[y*(*w)+x] = 0;
#ifdef RGBA
        pixel[y*(*w)+x] |= (rint&0xFF) << 24;
        pixel[y*(*w)+x] |= (gint&0xFF) << 16;
        pixel[y*(*w)+x] |= (bint&0xFF) << 8;
        pixel[y*(*w)+x] |= 0xFF;
#else
        pixel[y*(*w)+x] |= (rint&0xFF);
        pixel[y*(*w)+x] |= (gint&0xFF) << 8;
        pixel[y*(*w)+x] |= (bint&0xFF) << 16;
        pixel[y*(*w)+x] |= 0xFF << 24;
#endif
      }
      else {
        if (type == 6)
          fscanf(fp,"%c%c%c",&ruchar,&guchar,&buchar);
        else {
          fscanf(fp,"%c",&ruchar);
          guchar=buchar=ruchar;
        }
        pixel[y*(*w)+x] = 0;
#ifdef RGBA
        pixel[y*(*w)+x] |= (ruchar&0xFF) << 24;
        pixel[y*(*w)+x] |= (guchar&0xFF) << 16;
        pixel[y*(*w)+x] |= (buchar&0xFF) << 8;
        pixel[y*(*w)+x] |= 0xFF;
#else
        pixel[y*(*w)+x] |= (ruchar&0xFF);
        pixel[y*(*w)+x] |= (guchar&0xFF) << 8;
        pixel[y*(*w)+x] |= (buchar&0xFF) << 16;
        pixel[y*(*w)+x] |= 0xFF << 24;
#endif
      }
    }
  }
  return pixel;
} /* iiReadPPMFile */

/* -----------------------------------------------------------------------------
NAME:     iiWritePPMFile
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int iiWritePPMFile (char *filename, unsigned int *pixel, int width, int height)
{
  FILE        *fp = NULL;
  int          x, y;
  unsigned int ind;

  /* NOTE All of your values better be <= 255! */

  if (!pixel)
    return 0;
  fp = fopen(filename,"wt");
  if (fp == NULL)
    return -1;
  fprintf(fp,"P6\n");
  fprintf(fp,"%d %d\n",width,height);
  fprintf(fp,"255\n");
  for (y=0; y < height; y++) {
    for (x=0; x<width; x++) {
      ind = y*width+x;
      fprintf(fp,"%c",(char)((pixel[ind])&0xFF));
      fprintf(fp,"%c",(char)((pixel[ind]>>8)&0xFF));
      fprintf(fp,"%c",(char)((pixel[ind]>>16)&0xFF));
    }
  }
  fclose(fp);
  return 0;
} /* iiWritePPMFile */

/* -----------------------------------------------------------------------------
NAME:
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int iiWritePGMFileFloat (char *filename, float *pixel, int width, int height)
{
  int          x, y;
  unsigned int ind;
  FILE        *fp = fopen(filename,"wt");

  /* NOTE All of your values better be <= 255! */

  if (!fp)
    return -1;
  fprintf(fp,"P5\n");
  fprintf(fp,"%d %d\n",width,height);
  fprintf(fp,"255\n");
  for (y=0; y < height; y++) {
    for (x=0; x< width; x++) {
      ind = y*width+x;
      fprintf(fp,"%c",(unsigned char)(pixel[ind]));
    }
  }
  fclose(fp);
  return 0;
} /* iiWritePGMFileFloat */

/* -----------------------------------------------------------------------------
NAME:
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int iiWritePGMFileInt (char *filename, unsigned char *pixel, int width,
                       int height)
{
  FILE        *fp = NULL;
  int          x, y;
  unsigned int ind;

  /* NOTE All of your values better be <= 255! */

  fp = fopen(filename,"wt");
  if (!fp)
    return(MFALSE);
  fprintf(fp,"P5\n");
  fprintf(fp,"%d %d\n",width,height);
  fprintf(fp,"255\n");
  for (y=0; y < height; y++) {
    for (x=0; x< width; x++) {
      ind = y*width+x;
      fprintf(fp,"%c",pixel[ind]);
    }
  }
  fclose(fp);
  return(MTRUE);
} /* iiWritePGMFileInt */

/* -----------------------------------------------------------------------------
NAME:     iiComputeHistogram
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
void iiComputeHistogram (unsigned char *pixel, int width, int height,
                         unsigned int *histogram, unsigned int histsize)
{
  unsigned int arraysize;
  unsigned int i;

  memset(histogram, 0, histsize*sizeof(unsigned int));

  arraysize = width*height;
  for (i=0;i<arraysize;i++) {
    if (pixel[i] < histsize)
      (histogram[pixel[i]])++;
  }
} /* iiComputeHistogram */
/* -----------------------------------------------------------------------------
NAME:     iiWriteHistogramFile
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int iiWriteHistogramFile (char *filename, unsigned char *pixel, int width,
                          int height)
{
  FILE        *fp = NULL;
  unsigned int arraysize;
  unsigned int i;
  unsigned int histogram[256];

  /* NOTE All of your values better be <= 255! */

  fp = fopen(filename,"w");
  if (!fp)
    return(MFALSE);

  iiComputeHistogram(pixel, width, height, histogram, 256);

  for (i=0;i<256;i++)
    fprintf(fp, "%d\n", histogram[i]);

  fclose(fp);
  return(MTRUE);
} /* iiWriteHistogramFile */

/* -----------------------------------------------------------------------------
NAME:     iiConvertToGrayscale
PURPOSE:  Converts an RGB image to grayscale.
PRE:      The image must be read into the image array.
POST:     The grayscale image should be freed.
RETURN:   The image in grayscale
AUTHOR:   CMS
DATE:     September 11, 1999
MODIFIED:
----------------------------------------------------------------------------- */
unsigned char *iiConvertToGrayscale(unsigned int *rgbimage, int w, int h)
{
  unsigned char *grayimage = NULL;
  unsigned int   ind;
  unsigned char  r, g, b;
  int            x, y;

  if (!rgbimage)
    return(NULL);

  grayimage = (unsigned char *)malloc(w*h*sizeof(unsigned char));
  if (!grayimage)
    return(NULL);

  for (y=0; y<h; y++) {
    for (x=0; x<w; x++) {
      ind = y*w+x;
      r = (unsigned char)(rgbimage[ind])&0xFF;
      g = (unsigned char)(rgbimage[ind]>>8)&0xFF;
      b = (unsigned char)(rgbimage[ind]>>16)&0xFF;
      grayimage[ind] = (unsigned char)((float)r*0.299+(float)g*0.587+
                                       (float)b*0.114);
      if (grayimage[ind] > 255)
        grayimage[ind] = 255;
    }
  }
  return(grayimage);
} /* iiConvertToGrayscale */
/* -----------------------------------------------------------------------------
NAME:     iiComputeLookup
PURPOSE:  Computes a lookup table for equalizing the giving histogram.
PRE:      The histogram must be computed.
POST:     None.
RETURN:   None.
AUTHOR:   CMS
DATE:     September 18, 1999
MODIFIED:
----------------------------------------------------------------------------- */
void iiComputeLookup(unsigned int *histogram, float *lookup, int width,
                     int height, unsigned int histsize)
{
  unsigned int i;
  float        cumprob;
  float        imagearea;

  cumprob = 0.0;
  memset(lookup, 0, histsize*sizeof(float));
  imagearea = (float)width*(float)height;
  for (i=0; i<histsize; i++) {
    cumprob += (float)histogram[i]/imagearea;
    lookup[i] = cumprob*((float)histsize-1.001);
  }
} /* iiComputeLookup */
/* -----------------------------------------------------------------------------
NAME:     iiEqualizeImage
PURPOSE:  Equalizes the histogram of the given image and modifies the image
          to represent this equalization.
PRE:      The image must exist.
POST:     None.
RETURN:   None.
AUTHOR:   CMS
DATE:     September 18, 1999
MODIFIED:
----------------------------------------------------------------------------- */
void iiEqualizeImage(char *filename, unsigned char *grayimage, int width,
                     int height)
{
  unsigned int      histogram[256];
  float             lookup[256];
  unsigned long int imagearea;
  unsigned long int i;
  FILE             *fp;

  fp = fopen(filename,"w");

  imagearea = width*height;
  iiComputeHistogram(grayimage, width, height, histogram, 256);
  iiComputeLookup(histogram, lookup, width, height, 256);
  if (fp) {
    for (i=0;i<256;i++)
      fprintf(fp, "%1.3lf\n", lookup[i]);
  }
  for (i=0;i<imagearea;i++)
    grayimage[i] = (unsigned char)lookup[grayimage[i]];

  if (fp)
    fclose(fp);
} /* iiEqualizeImage */
/* -----------------------------------------------------------------------------
NAME:     iiAverageGrayImages
PURPOSE:  Averages an array of grayscale images
PRE:
POST:
RETURN:   An array of unsigned chars
AUTHOR:
DATE:
MODIFIED: CMS-September 26, 1999
----------------------------------------------------------------------------- */
unsigned char *iiAverageGrayImages (int w, int h, unsigned char **grayimage,
                                    int numberofimages)
{
  unsigned char     *pixel = NULL;
  float             *imagesum = NULL, myfloat;
  unsigned long int  i, j, imagearea;

  imagearea = w*h;
  pixel = (unsigned char *)malloc(imagearea*sizeof(unsigned char));
  if (!pixel)
    return(pixel);
  memset(pixel, 0, imagearea*sizeof(unsigned char));
  imagesum = (float *)malloc(imagearea*sizeof(float));
  if (!imagesum)
    return(pixel);
  memset(imagesum, 0, imagearea*sizeof(float));

  for (i=0;i<numberofimages;i++) {
    if (grayimage[i]) {
      for (j=0;j<imagearea;j++) {
        imagesum[j] += (float)grayimage[i][j];
      }
    }
  }

  for (i=0;i<imagearea;i++) {
    myfloat = imagesum[i]/(float)numberofimages;
    pixel[i] = (unsigned char)(myfloat > (255.0-1.0e-30)) ? 255.0 : myfloat;
  }

  if (imagesum)
    free(imagesum);

  return pixel;
} /* iiAverageGrayImages */
/* -----------------------------------------------------------------------------
NAME:     iiiAverageSurroundingPixels
PURPOSE:  Given x and y locations on an image and the image, this function
          interpolates the grayscale value for the given x and y on the image.
          If the x and y are off the image, 0 is returned.
PRE:      None.
POST:     None.
RETURN:   The grayscale value.
AUTHOR:   CMS
DATE:     September 28, 1999
MODIFIED:
----------------------------------------------------------------------------- */
unsigned char iiiAverageSurroundingPixels(float u, float v,
                                          unsigned char *grayimage,
                                          int w, int h)
{
  long int ind1 = -1, ind2 = -1, ind3 = -1, ind4 = -1;
  int      i, j;

  if ((u < 0) || (u >= w) || (v < 0) || (v >= h))
    return(0);
  i = u;
  j = v;
  if ((i-1) >= 0)
    ind1 = j*w+i-1;
  if (((j-1) >= 0) && ((i-1) >= 0))
    ind2 = (j-1)*w+i-1;
  if ((j-1) >= 0)
    ind3 = (j-1)*w+i;
  if (((j-1) >= 0) && ((i+1) < w))
    ind4 = (j-1)*w+i+1;
  if ((ind1 != -1) && (ind2 != -1) && (ind3 != -1) && (ind4 != -1))
    return((unsigned char)((float)(grayimage[ind1]+grayimage[ind2]+
                                   grayimage[ind3]+grayimage[ind4])/4.0));
  else if (ind4 == -1)
    return((unsigned char)((float)(grayimage[ind1]+grayimage[ind2]+
                                   grayimage[ind3])/3.0));
  else
    return(0);
} /* iiiAverageSurroundingPixels */
/* -----------------------------------------------------------------------------
NAME:     iiEncodeGrayImage
PURPOSE:  Encodes a gray image.
PRE:      The image must be set up.
POST:     Must free the finalimage.
RETURN:   the encoded image.
AUTHOR:   CMS
DATE:     September 28, 1999
MODIFIED:
----------------------------------------------------------------------------- */
unsigned char *iiEncodeGrayImage(int w, int h, unsigned char *grayimage)
{
  int                u, v, cury, curx;
  unsigned long int  imagearea, ind;
  unsigned char     *pixel = NULL, averageval;

  imagearea = w*h;
  pixel = (unsigned char *)malloc(imagearea*sizeof(unsigned char));
  if (!pixel)
    return(pixel);
  memset(pixel, 0, imagearea*sizeof(unsigned char));

  for (u=0;u<w;u++) {
    for (v=0;v<h;v++) {
      ind = v*w+u;
      if ((u==0) || (v==0))
        pixel[ind] = grayimage[ind];
      else {
        averageval = iiiAverageSurroundingPixels(u, v, grayimage, w, h);
        pixel[ind] = grayimage[ind]-averageval;
      }
    }
  }
  return(pixel);
} /* iiEncodeGrayImage */
/* -----------------------------------------------------------------------------
NAME:     main
PURPOSE:  This is the main function, buckeroo!
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int main(int argc, char *argv[])
{
  unsigned char  *grayimage = NULL;
  unsigned char  *finalimage = NULL;
  FILE           *fp = NULL;
  int             w, h;

  if (argc != 3) {
    printf("Usage: encode <infilename> <outfilename>\n");
    return(1);
  }

  fp = fopen(argv[1], "r");
  if (!fp)
    return(2);

  grayimage = iiReadPGMFile(fp, &w, &h);
  if (fp)
    fclose(fp);
  finalimage = iiEncodeGrayImage(w, h, grayimage);

/* CMS_COMMENT Mon Sep 27 08:27:04 MDT 1999
  finalimage = iiAverageGrayImages(w, h, grayimage, numberoffiles);
  iiWriteHistogramFile("out1.his", grayimage, w, h);
  iiEqualizeImage("out.cum", grayimage, w, h);
  iiWriteHistogramFile("out2.his", grayimage, w, h);
CMS_COMMENT */

  iiWritePGMFileInt(argv[2], finalimage, w, h);
  if (grayimage)
    free(grayimage);
  if (finalimage)
    free(finalimage);
  return(0);
} /* main */