#include <stdio.h>
#include <math.h>
#include <unistd.h>
#include "string.h"
#ifndef F_OK
# define F_OK        0       /* Test for existence of file */
#endif /* F_OK */

void hsvrgb(float,float,float,int *,int *,int *);
void cmap1(float *,float *,float *,float);
void cmap2(float *,float *,float *,float);

/*
 * Program to convert raw binary data to PGM or PPM files.
 * The approach used here was described in the following article:
 *
 *   J. B. Schneider, P. J. Flynn, and K. L. Shlager, "Animating the
 *   Evolution of a Field," IEEE Antennas and Propagation Magazine,
 *   vol. 38, no. 6, pp. 7--17, 1996.
 *
 * Feel free to do whatever you want with this code provided you obtain
 * consent prior to any redistribution or commercial use.
 *
 * John B. Schneider
 * School of Electrical Engineering and Computer Science
 * Washington State University
 * Pullman, WA 99164-2752
 * schneidj@eecs.wsu.edu
 * (509) 335-4655, FAX: (509) 335-3818
 * http://www.eecs.wsu.edu/~schneidj/
 */

main(){
  float fcol, frow, h, s, v, xpnt, xnorm=-1., xj, eps, xtmp;
  int i, j, irow, icol, nframe=0, ndec=-1, icolor=-1, iside=-1;
  char red, green, blue, gray;
  int ir, ig, ib;
  char base[1024], filein[1024], filout[1024];
  FILE *in, *out;
/*
 * h,s,v    = hue, saturation, and brightness (or "value") for pixel
 * ir,ig,ib = corresponding red, green, and blue values for pixel
 * eps      = "small" number added to data to prevent underflows
 * xnorm    = number used to normalize all input
 * xpnt     = the current field value obtained from the input file
 * xtmp     = temporary variable for holding result
 * i,j      = dummy counters
 * irow,icol= number of rows and columns of data (integer)
 * irow,icol= number of rows and columns of data (float)
 * nframe   = current frame number
 * ndec     = number of decades of compression
 * base     = base input and output file names
 * filein,filout = complete input and output file names
 * iside    = used to indicate two-sided or one-sided data
 * icolor   = used to indicate color or grayscale
 */
  printf("Enter the base name of the raw files: ");
  scanf("%s",base);
  while(xnorm <= 0.) {
    printf("Enter the normalization: ");
    scanf("%f",&xnorm);
  }
  while(ndec < 0) {
    printf("Enter the number of decades (0 = linear scale): ");
    scanf("%d",&ndec);
  }
  while (icolor != 1 && icolor != 2) {
    printf("Do you want color or grayscale output? (1=grayscale; 2=color) ");
    scanf("%d",&icolor);
    if (icolor == 2) {
      while (iside != 1 && iside != 2) {
	printf("Do you want ""one-sided"" output (obtained by taking the ");
	printf("absolute value\nof the input data)?");
	printf(" (1=one-sided; 2=two-sided) ");
	scanf("%d",&iside);
      }
    }
  }
  /* Calculate the "epsilon" that will be added to values if compressing. */
  if (ndec != 0) eps = xnorm/(1.0-1./pow(10.0,ndec))/pow(10.0,ndec);
  /* Loop through all the frames. */
  do {
    nframe++;
    /* Construct the input file name and check that it exists. */
    sprintf(filein,"%s.%d",base,nframe);
    if (access(filein,F_OK)) {
      printf("Input not found: %s",filein);
      printf(" --> My job is done...\n");
      exit();
    }
    /* Construct the output file name and make sure it isn't there. */
    if (icolor == 1) {
      sprintf(filout,"%s.pgm",filein);
    } else {
      sprintf(filout,"%s.ppm",filein);
    }
    if (!(access(filout,F_OK))) {
      printf("Output file for %s already exists.  ",filein);
      printf("Skipping...\n");
    } else {
      if (!(out = fopen(filout,"wb"))) {
	printf("Ouch!  Couldn't open output file: %s\n",filout);
	printf("Terminating...\n");
	break;
      }
      if (!(in = fopen(filein,"rb"))) {
	printf("Ouch!  Couldn't open input file: %s\n",filein);
	printf("Terminating...\n");
	exit();
      }
      /* Read the size of the image. */
      printf("Processing %s\n",filein);
      /* Dimensions. */
      fread(&fcol,sizeof(float),1,in);
      fread(&frow,sizeof(float),1,in);
      irow = frow;
      icol = fcol + 10;
      /* Write the header information. */
      if (icolor == 1) {
	fprintf(out,"P5\n%d %d\n255\n",icol,irow);
      } else {
	fprintf(out,"P6\n%d %d\n255\n",icol,irow);
      }
      /* Loop through all pixels in the image. */
      for(j=1;j<=irow;j++){
	for(i=1;i<=icol;i++){
	  if (i <= 10) {
	    /* Draw the color map if within the first 10 columns. */
	    if (icolor == 1 || iside == 1) {
	      xpnt = ((float)(irow-j))/((float)(irow-1))*xnorm;
	    } else {
	      xpnt = (-2./((float)(irow-1))*(float)(j)
		      +((float)(irow+1))/((float)(irow-1)))*xnorm;
	    }
	  } else {
	    /* If not in first 10 columns, get data from file. */
	    fread(&xpnt,sizeof(float),1,in);
	  }
	  /* If the point is out of range, whack it down to size. */
	  if (fabs(xpnt) > xnorm) xpnt = xpnt/fabs(xpnt)*xnorm;
	  /* Use log scaling if requested. */
	  if (ndec != 0) {
	    xtmp = log10((fabs(xpnt)+eps)/(xnorm+eps))/((float)(ndec))+1.;
	  } else {
	    xtmp = fabs(xpnt/xnorm);
	  }
	  if (icolor == 1) {
	    /* Grayscale output. */
	    gray = (int)(255*xtmp);
	    fwrite(&gray,sizeof(char),1,out);
	  } else {
	    if (iside == 1) {
	      /* Use color map #1 for one-sided data. */
	      cmap1(&h,&s,&v,xtmp);
	    } else {
	      /* Reconstruct sign and use color map #2 for two-sided data. */
	      if (xpnt < 0.) xtmp = -xtmp;
	      cmap2(&h,&s,&v,xtmp);
	    }
	    /* Convert (h,s,v) to (r,g,b) and write output. */
	    hsvrgb(h,s,v,&ir,&ig,&ib);
	    red = ir; green = ig; blue = ib;
	    fwrite(&red,sizeof(char),1,out);
	    fwrite(&green,sizeof(char),1,out);
	    fwrite(&blue,sizeof(char),1,out);
	  }
	}
      }
      fclose(in); fclose(out);
    }
  } while (1);
}