// Kevin Beason
// http://www.cs.fsu.edu/~beason/
// 4/16/02
//
// This is my data set manipulation tool. You can decimate data sets and
// extract slices for 2d-viewing. It's pretty fast.
// 
// To compile:
// 
//   gcc -o deci deci.cxx
	
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <fstream.h>
#include <iostream.h>
#include <string.h>

unsigned int readVal(unsigned char *ptr, int valSize, int bigEndian=0){
  int i;
  unsigned int c = 0;

  if (! bigEndian){
    for (i=valSize-1; i>=0; i--){
      c = (c << 8) | ptr[i];
    }
  }
  else {
    for (i=0; i<valSize; i++){
      c = (c << 8) | ptr[i];
    }
  }

  if (c==0xf830)
    c=0;
  
  return c;
}

int main(int argc, char* argv[]){
  char previewFileName[128] = "out.pgm";
  const int maxDim=512;
  const int maxValSize=4;
  
  int stats[66] = { 0 };
  int i, j, k, m;
  unsigned char inbuf[maxDim*maxDim*maxValSize];
  unsigned char outbuf[maxDim*maxDim*maxValSize];
  unsigned int c;
  int numBytes;
  int numc=0, min=HUGE, max=0;
  int rows=0, cols=0, planes=0;

  // arguments
  int header = 0;
  int valSize = 1;
  char inFileName[128] = {0};
  char outFileName[128] = {0}; 
  int inDims[3] = { 0 };
  int outDims[3] = { 0 };
  int outFactor = 0;
  int slice = 0;
  int bigEndian = 0;

  // read arguments
  if (argc <= 1){
    cerr << "\nUsage: " << argv[0] << " -i <filename> <x> <y> <z> "
	 << "[ -n numBytes ] [ -b (bigEndian) ]\n"
	 << " [ -h header ] [ -o <filename> <intfactor> ] [ -s <slice> ]\n"
	 << endl;
    exit(1);
  }
  
  int argn;
  for (argn=1; argn<argc; argn++){
    if (!strcmp(argv[argn], "-i")){
      strcpy(inFileName, argv[++argn]);
      inDims[0] = atoi(argv[++argn]);
      inDims[1] = atoi(argv[++argn]);
      inDims[2] = atoi(argv[++argn]);
    }
    else if (!strcmp(argv[argn], "-o")){
      strcpy(outFileName, argv[++argn]);
      outFactor = atoi(argv[++argn]);
      if (outFactor <= 0){
	cerr << "Error: factor option to -o switch must be an integer > 0."
	     << endl;
	exit(1);
      }
      outDims[0] = inDims[0]/outFactor;
      outDims[1] = inDims[1]/outFactor;
      outDims[2] = inDims[2]/outFactor;
    }
    else if (!strcmp(argv[argn], "-n")){
      valSize = atoi(argv[++argn]);
    }
    else if (!strcmp(argv[argn], "-b")){
      bigEndian = 1;
    }
    else if (!strcmp(argv[argn], "-h")){
      header = atoi(argv[++argn]);
    }
    else if (!strcmp(argv[argn], "-s")){
      slice = atoi(argv[++argn]);
    }
  }
  if (! (outFactor || slice) ){
    cerr << "Error: Must use at least one of -o or -s switches.\n";
    exit(1);
  }
  
  // open input file
  FILE *inFile = fopen(inFileName, "rb");
  if (!inFile){
    cerr << "Could not open " <<  inFileName
         << " for reading.\n";
    exit(1);
  }

  // open output file
  FILE *outFile;
  if (outFactor){
    outFile = fopen(outFileName, "wb");
    if (!outFile){ 
      cerr << "Could not open " << outFileName
	   << " for writing.\n";
      exit(1);
    }
  }

  // open preview file
  ofstream preview;  
  if (slice){
    preview.open(previewFileName);
    if (!preview){
      cerr << "Could not open " << previewFileName
	   << " for writing.\n";
    }
  }

  // skip header
  if (header){
    i = maxDim*maxDim*maxValSize;
    fread(inbuf, 1, header%i, inFile);
    for (j=0; j<header/i; j++)
      fread(inbuf, 1, i, inFile);
  }

  // loop through data
  for (k=0; k<inDims[2]; k++){
    
    // read slice
    cerr << "Reading slice " << k+1 << "... ";
    cerr << fread(inbuf, valSize, inDims[0]*inDims[1], inFile)
	 << " values read.\n";
    
    // decimate
    if (outFactor && k%outFactor==0){
      numBytes=0;
      rows=0;
      for (j=0; j<inDims[1]; j+=outFactor){
	rows++;
	cols=0;
	for (i=0; i<inDims[0]; i+=outFactor){
	  cols++;
	  for (m=0; m<valSize; m++){
	    outbuf[numBytes] = inbuf[ (j*inDims[0]+i)*valSize + m ];
	    numBytes++;
	  }
	}
      }
      fwrite(outbuf, valSize, numBytes/valSize, outFile);
      planes++;
    }
    
    // write slice
    if (k==slice-1){

      // determine max value
      for (j=0; j<inDims[1]; j++)
	for (i=0; i<inDims[0]; i++){
	  c = 0;
	  for (m=0; m<valSize; m++){
	    c = readVal(&inbuf[(j*inDims[0]+i)*valSize], valSize, bigEndian);
	  }
	  if (max<c)
	    max = c;
	}

      // write ppm header
      preview << "P2\n\n" << inDims[0] << " " << inDims[1] << "\n"
	      << max << "\n";
      
      // write
      cerr << "Writing.\n";
      for (j=0; j<inDims[1]; j++)
	for (i=0; i<inDims[0]; i++){
	  c = readVal(&inbuf[(j*inDims[0]+i)*valSize], valSize, bigEndian);

	  // write data
	  preview << c << " ";

	  // stats
	  numc++;
	  stats[c/46]++;
	  if (min>c) min = c;
	  if (max<c) max = c;
	}

      // quit loop if not writing new one
      if (!outFactor)
	break;
    }
  }

  fclose(inFile);
  fclose(outFile);
  preview.close();

  cerr << "Wrote " << numc << " values to " << outFileName << "." << endl;
  cerr << "Stats:\n";
  for (i=0;i<66;i++)
    cerr << " i=" << stats[i];
  cerr << "\nMin: " << min << endl << "Max: " << max << endl;
  cerr << "Decimated file set is " << cols << "x" << rows << "x" << planes
       << endl;
  cerr << endl;
  
  return(0);
}