/**
 * The idea behind <TT> fromFloatsToCubes.c </TT> is to read a set of
 * volume data (stored as floating point numbers) at uniform resolution 
 * and output that same resolution represented as cubes.
 * Note that this process will increase the data set size by almost
 * 8 times since all of the interior vertices will be shared amongst
 * 8 cubes.
 * <P>
 * Note that the number of sample points along each dimension is 1 more
 * than the number of cubes.  For example, 4 rows of 129 sample points
 * are needed to make 1 row of 128 cubes.
 * <P>
 * @file       fromFloatToCube.c
 * @author     Robert S Laramee
 * @start date Wed 12 July 00
 */

#include <stdio.h>
#include <stdlib.h>  /*  for atoi(), atof()  */
#include <string.h>

#define INVALID -2
#define TOOBIG  1000

int lengthOfRows    = INVALID;
int heightOfColumns = INVALID;
int depthOfLayers   = INVALID;
int currentOctant   = INVALID;

/**
 * @return the length of one row 
 *         (the number of sample points along the x dimension)
 */
int getLengthOfRows() { return lengthOfRows; }
/**
 * @param the length of one row 
 *        (the number of sample points along the x dimension)
 */
void setLengthOfRows(int newLength) { 

    if ((newLength < 0) || (newLength > TOOBIG)) {
        fprintf(stderr, "***Error, setLengthOfRows(): invalid number of "
		        "rows: %d, exiting...\n", newLength);
        exit(1);
    } else {
	lengthOfRows = newLength;
    }
}
/**
 * @return the height of one column
 *         (the number of sample points along the y dimension)
 */
int getHeightOfColumns() { return heightOfColumns; }
/**
 * @param the height of one column 
 *        (the number of sample points along the y dimension)
 */
void setHeightOfColumns(int newHeight) { 

    if ((newHeight < 0) || (newHeight > TOOBIG)) {
        fprintf(stderr, "***Error, setHeightOfColumns(): invalid number of "
		        "columns: %d, exiting...\n", newHeight);
        exit(1);
    } else {
	heightOfColumns = newHeight; 
    }
}
/**
 * @return the depth of the layers as a collection
 *         (the number of sample points along the z dimension)
 */
int getDepthOfLayers() { return depthOfLayers; }
/**
 * @param the depth of the layers as a collection 
 *        (the number of sample points along the z dimension)
 */
void setDepthOfLayers(int newDepth) { 

    if ((newDepth < 0) || (newDepth > TOOBIG)) {
        fprintf(stderr, "***Error, setDepthOfLayers(): invalid number of "
		        "layers: %d, exiting...\n", newDepth);
        exit(1);
    } else {
	depthOfLayers = newDepth; 
    }
}
/**
 * @return the current octant we are outputting -level 0 is output
 *         1 octant at a time
 */
int getOutputOctant() { return currentOctant; }
/**
 * @param the current octant we are outputting -level 0 is output
 *        1 octant at a time
 */
void setOutputOctant(int newOctant) { 

    if ((newOctant < -2) || (newOctant > 8)) {
        fprintf(stderr, "***Error, setOctant(): invalid octant "
		        ": %d, exiting...\n", newOctant);
        exit(1);
    } else {
	currentOctant = newOctant; 
    }
}
/**
 * @return the number of sample points in one layer
 */
int getSizeOfLayer() { return (getLengthOfRows() * getHeightOfColumns()); }

/**
 * Make sure you declare this structure before its used, even in the
 * function prototypes.
 */
struct cube {
    /* a cube has 8 vertex values */
    float p0;  float p1;  float p2;  float p3;
    float p4;  float p5;  float p6;  float p7;

    int xIndex; /* and 3 index values indicating the x, y, and z */
    int yIndex; /* coordinates (in block space) of its 0th       */
    int zIndex; /* vertex                                        */
};

/* function prototypes */
int  outputVolumeOfCubes(FILE *inputFile, FILE* inputErrFile);
int  outputLayerOfCubes(float layerFront[], float layerBack[],
			float layerFrontErr[],float layerBackErr[],
			int currentLayer);
int  outputRowOfCubes(float frontLayerOfData[], float  backLayerOfDate[], 
		      float frontLayerOfErr[], float backLayerOfErr[],
                      int layerOffset, int currentLayer);
int  readLayerOfFloats(FILE *inputFile,float layer[]);
int  copyLayer(float sourceLayer[], float destinationLater[]);
void printLayerOfFloats(float layer[]);
int  whichOctant(int xCoord, int yCoord, int zCoord);
int readLayerOfError(FILE *inputErrorFile, float error[]);

/**
 * some handy numbers:
 * <PRE>
 * level  number of sample points
 * -----  -----------------------
 *  7     2 x 2 x 2       ( 2^{0} + 1 )
 *  6     3 x 3 x 3       ( 2^{1} + 1 )
 *  5     5 x 5 x 5       ( 2^{2} + 1 )
 *  4     9 x 9 x 9       ( 2^{3} + 1 )
 *  3     17 x 17 x 17    ( 2^{4} + 1 )
 *  2     33 x 33 x 33    ( 2^{5} + 1 )
 *  1     65 x 65 x 65    ( 2^{6} + 1 )
 *  0     129 x 129 x 129 ( 2^{7} + 1 )
 * </PRE>
 * Compile the program using:
 * <PRE>
 * % gcc -Wall fromFloatToErrCube.c -Iinclude -lm -o fromFloatToErrCube.exe
 * </PRE>
 * <P>
 * Start the program using:
 * <PRE>
 * % fromFloatToErrCube.exe [input file] [x] [y] [z] [octant] > [output file] 
 *   e.g.
 * % fromFloatToErrCube.exe ~rlaramee/data/level5/headLev5floats.ascii \
 *   5 5 5 -2 > ~rlaramee/data/level5/headLev5cubes.ascii
 *
 * </PRE>
 * <P>
 * @param argv[1] -the input data file name
 * @param argv[2] -the number of data points along the x dimension
 * @param argv[3] -the number of data points along the y dimension
 * @param argv[4] -the number of data points along the z dimension
 * @param argv[5] -the octant to output, level 0 is output 1 octant
 *                 at a time
 */
int main(int argc , char *argv[]) {

    FILE *inputFile;
    FILE *inputErrorFile;
    char *inputFileName;
    char inputErrorFileName[60];
    char *s;
    int totalOutput = 0;

    /*  check command line arguments  */
    if (argc != 6) {
	fprintf(stderr, "*** Error, usage: fromFloatToCube.exe [inputfile] "
		        "[x] [y] [z] [octant] > [output file] \n");
	exit(1);
    }

    /* We'll build the error file name from the inputFileName,
     * which came from the command line.
     */
    inputFileName = argv[1];    
    strcpy(inputErrorFileName, inputFileName);
    s = strstr(inputErrorFileName,"floats.ascii");
    strcpy(s,"Error.ascii");


    if ( !( inputFile = fopen(inputFileName,"rb"))){
        fprintf(stderr, "***Error, couldn't open input file: %s \n", 
		         inputFileName);
        exit(1);
    }
    
    if ( !( inputErrorFile = fopen(inputErrorFileName,"rb"))){
        fprintf(stderr, "***Error, couldn't open error file: %s \n", 
		         inputErrorFileName);
        exit(1);
    }

    setLengthOfRows(   atoi(argv[2]));
    setHeightOfColumns(atoi(argv[3]));
    setDepthOfLayers(  atoi(argv[4]));    
    setOutputOctant(   atoi(argv[5]));    

    fprintf(stderr, "# fromFloatToErrCube.main() input file is: %s\n", 
	    inputFileName);
    fprintf(stderr, "# fromFloatToErrCube.main() error file is: %s\n", 
	    inputErrorFileName);
    fprintf(stderr, "# fromFloatToErrCube.main() with x, y, z dimensions: "
	            "%d, %d, %d\n", 
	    getLengthOfRows(), getHeightOfColumns(), getDepthOfLayers());
    fprintf(stderr, "# fromFloatToErrCube.main() outputting octant: "
	            "%d, \n", getOutputOctant());

    totalOutput = outputVolumeOfCubes(inputFile,inputErrorFile);

    fclose(inputFile);
    fclose(inputErrorFile);

    fprintf(stderr, "# main(): output: %d cubes \n", totalOutput);
    return 0;
}
/**
 * This function outputs the volume data in cube representation.
 * <P>
 * @param inputFile a handle to the input file
 * @return the total number of cubes output
 */
int outputVolumeOfCubes(FILE *inputFile,FILE *inputErrorFile) {

    int debug       = 1;
    int z           = 0;
    int totalRead   = 0;
    int totalOutput = 0;
    float layerOfData1[getSizeOfLayer()];
    float layerOfData2[getSizeOfLayer()];
    float layerOfError1[getSizeOfLayer()];
    float layerOfError2[getSizeOfLayer()];

    if (debug == 1) 
	fprintf(stderr, "# fromFloatToCube.outputVolumeOfCubes()\n");

    /** FOR EACH LAYER */
    for (z = 0; z < (getDepthOfLayers() - 1); z++) {

	if ( z == 0 ) {  /* the 1st time, read 2 layers */

	    totalRead = totalRead + readLayerOfFloats(inputFile, layerOfData1);
	    readLayerOfError(inputErrorFile, layerOfError1);

	    totalRead = totalRead + readLayerOfFloats(inputFile, layerOfData2);
	    readLayerOfError(inputErrorFile, layerOfError2);

	} else {

	    copyLayer(layerOfData2, layerOfData1);
	    copyLayer(layerOfError2, layerOfError1);

	    totalRead = totalRead + readLayerOfFloats(inputFile, layerOfData2);
	    readLayerOfError(inputErrorFile,layerOfError2);
	} /* end if (z == 0) */

	/** output every layer */
	totalOutput = totalOutput + outputLayerOfCubes(layerOfData1, 
						       layerOfData2,
						       layerOfError1,
						       layerOfError2,
						       z);

    } /* end FOR EACH LAYER */
    if (debug == 1) fprintf(stderr, "# fromFloatToCube.outputVolumeOfCubes()"
			    " returning: %d\n", totalRead);
    return totalOutput;
}

/**
 * @param frontLayerOfFloats the front layer of scalar values
 * @param backLayerOfFloats  the back layer of scalar values
 * @param currentLayer       the current layer number of cubes we are 
 *                           outputting
 * @return the total number of cubes output in the layer
 */
int outputLayerOfCubes(float frontLayerOfFloats[], float backLayerOfFloats[], 
		       float frontLayerOfError[],  float backLayerOfError[], 
		       int currentLayer) {

    int debug       = 0;
    int y           = 0;
    int rowOffset   = 0;
    int totalOutput = 0;

    if (debug == 1) fprintf(stdout, "# fromFloatToCube.outputLayerOfCubes()" 
			    " frontLayerOfFloats = \n");
    if (debug == 1) printLayerOfFloats(frontLayerOfFloats);
    if (debug == 1) fprintf(stdout, "# fromFloatToCube.outputLayerOfCubes()" 
			    " backLayerOfFloats = \n");
    if (debug == 1) printLayerOfFloats(backLayerOfFloats);
 
    /* FOR EACH ROW - 1 (# sample points in column = # of cubes + 1) */
    for (y = 0; y < (getHeightOfColumns() - 1); y++) {

	/* index into every row */
	rowOffset = y * getLengthOfRows();
	
	/* output every row */
	totalOutput = totalOutput + 
	    outputRowOfCubes(frontLayerOfFloats, backLayerOfFloats,
			     frontLayerOfError,  backLayerOfError, 
			     rowOffset, currentLayer);

    } /* end FOR EACH ROW */
    if (debug > 0) fprintf(stdout, "\n");
    return totalOutput;
}

/**
 * This method outputs the 8 interpolated data values.  We have to be
 * careful to output the cube vertices in the correct order.
 * <PRE>
 *              P3_____________P2   This is the cube representation
 *              /|            /|    for Bob's cube's.  It is the
 *             / |           / |    same as the VTK's.
 *    y     P7____________P6/  |    
 *    ^      |   |         |   |    
 *    |      |   |         |   |    
 *    |      |   |         |   |
 *    |--->  | P0|_________|___|P1
 *   /    x  |  /          |  /
 *  /        | /           | /
 * z         |_____________|/
 *          P4             P5
 * </PRE>
 * <P>
 * @param frontLayerOfData the front layer of vertex scalars
 * @param backLayerOfData  the back layer of vertex scalars
 * @param offsetIntoLayer  the number of rows already processed * 
 *                         the length of each row
 * @param currentLayer     the current layer number of cubes we are 
 *                         outputting
 * @return the total number of cubes output in that row
 */
int outputRowOfCubes(float frontLayerOfData[],  float backLayerOfData[], 
		     float frontLayerOfError[], float backLayerOfError[],
		     int offsetIntoLayer, int currentLayer) {

    int debug     = 0;
    int x         = 0;
    int y         = 0;
    int z         = 0;
    int numOutput = 0;

    if (debug == 1) fprintf(stdout, "# fromFloatToCube.outputRowOfCubes() "
			    "offsetIntoLayer = %d, currentLayer = %d \n", 
			    offsetIntoLayer, currentLayer);

    /* FOR EACH cube (# of sample points in row = # cubes in row + 1 */
    for (x = 0; x < (getLengthOfRows() - 1); x++) {

	y = offsetIntoLayer / getLengthOfRows();
	z = currentLayer;

	/* only output if the cube is in the given octant */
	if ( (whichOctant(x, y, z) == getOutputOctant()) ||
	     (getOutputOctant() < 0) 
	   ) {

	    /* print a cube */

	    fprintf(stdout, "%f %f %f %f %f %f %f %f\n" 
			    "%f %f %f %f %f %f %f %f\n",

     /* p0 */ frontLayerOfData[offsetIntoLayer + x ],
     /* p1 */ frontLayerOfData[offsetIntoLayer + x + 1 ],
     /* p2 */ frontLayerOfData[offsetIntoLayer + x + getLengthOfRows() + 1 ],
     /* p3 */ frontLayerOfData[offsetIntoLayer + x + getLengthOfRows() ],

     /* p4 */  backLayerOfData[offsetIntoLayer + x ],
     /* p5 */  backLayerOfData[offsetIntoLayer + x + 1 ],
     /* p6 */  backLayerOfData[offsetIntoLayer + x + getLengthOfRows() + 1],
     /* p7 */  backLayerOfData[offsetIntoLayer + x + getLengthOfRows() ],
     /* Now the error */
     /* p0 */ frontLayerOfError[offsetIntoLayer + x ],
     /* p1 */ frontLayerOfError[offsetIntoLayer + x + 1 ],
     /* p2 */ frontLayerOfError[offsetIntoLayer + x + getLengthOfRows() + 1 ],
     /* p3 */ frontLayerOfError[offsetIntoLayer + x + getLengthOfRows() ],

     /* p4 */  backLayerOfError[offsetIntoLayer + x ],
     /* p5 */  backLayerOfError[offsetIntoLayer + x + 1 ],
     /* p6 */  backLayerOfError[offsetIntoLayer + x + getLengthOfRows() + 1],
     /* p7 */  backLayerOfError[offsetIntoLayer + x + getLengthOfRows() ]);
		  
	    fprintf(stdout, "%d %d %d\n\n",  /* x, y, z indices */
		    x, y, z);
	    
	    numOutput++;

	} /*  end if whichOctant()  */
    }     /*  end FOR EACH XROW     */

    if (debug == 1) fprintf(stdout, "\n");
    return numOutput;
}
/**
 * This method decides which octant the cube belongs to  
 * The decision(s) are based on the VTK representation of the cube.
 * This method is used to split up headLevel0cubes.ascii into
 * 8 binary cube files -1 for each octant.
 * @param cube    -the current cube
 * @return octant -the octant in which the cube belongs
 */
int whichOctant(int xCoord, int yCoord, int zCoord) {

    /* the -1 is important.  We've told the program that the
     * length of a row is 129, however, the length of a row
     * is 128 (at level 0) for this computation based
     * on the location of vertex 0
     */
    float maximumBlockCoord = (getLengthOfRows() - 1);

    if ((xCoord < (maximumBlockCoord/2.0)) &&   /* octant 0 */
	(yCoord < (maximumBlockCoord/2.0)) &&
	(zCoord < (maximumBlockCoord/2.0)) ) {
	return 0;

    } else if ((xCoord >= (maximumBlockCoord/2.0)) &&  /* octant 1 */
	       (yCoord <  (maximumBlockCoord/2.0)) &&
	       (zCoord <  (maximumBlockCoord/2.0)) ) {
	return 1;

    } else if ((xCoord >= (maximumBlockCoord/2.0)) &&  /* octant 2.0 */
	       (yCoord >= (maximumBlockCoord/2.0)) &&
	       (zCoord <  (maximumBlockCoord/2.0)) ) {
	return 2;
	
    } else if ((xCoord <  (maximumBlockCoord/2.0)) &&  /* octant 3 */
	       (yCoord >= (maximumBlockCoord/2.0)) &&
	       (zCoord <  (maximumBlockCoord/2.0)) ) {
	return 3;
	
    } else if ((xCoord <  (maximumBlockCoord/2.0)) &&  /* octant 4 */
	       (yCoord <  (maximumBlockCoord/2.0)) &&
	       (zCoord >= (maximumBlockCoord/2.0)) ) {
	return 4;
	
    } else if ((xCoord >= (maximumBlockCoord/2.0)) &&  /* octant 5 */
	       (yCoord <  (maximumBlockCoord/2.0)) &&
	       (zCoord >= (maximumBlockCoord/2.0)) ) {
	return 5;
	
    } else if ((xCoord >= (maximumBlockCoord/2.0)) &&  /* octant 6 */
	       (yCoord >= (maximumBlockCoord/2.0)) &&
	       (zCoord >= (maximumBlockCoord/2.0)) ) {
	return 6;
	
    } else if ((xCoord <  (maximumBlockCoord/2.0)) &&  /* octant 7 */
	       (yCoord >= (maximumBlockCoord/2.0)) &&
	       (zCoord >= (maximumBlockCoord/2.0)) ) {
	return 7;
	
    } else {
	fprintf(stderr, "# *** Error, build_mr_v3.whichOctant()");
	return -1;
    }
}
/**
 * @param sourceLayer      the layer of floats containing the data to be copied
 * @param destinationLayer the layer to store the copied data in
 * @return numCopied -the number of scalars copied from the source layer to
 *                    the destination layer
 */
int copyLayer(float sourceLayer[], float destinationLayer[]) {

    int i         = 0;
    int numCopied = 0;

    for (i = 0; i < getSizeOfLayer(); i++) {

	destinationLayer[i] = sourceLayer[i];
	numCopied++;
    }

    return numCopied;
}
/**
 * @param inputFile a handle to the input file
 * @param data      the array to store the data in
 */
int readLayerOfFloats(FILE *inputFile, float data[]) {

    int i = 0;

    for (i = 0; i < getSizeOfLayer(); i++) {

	fscanf(inputFile, "%f", &data[i]);
    }
    return getSizeOfLayer();
}
/**
 * @param inputErrorFile a handle to the input error file
 * @param data      the array to store the data in
 */
int readLayerOfError(FILE *inputErrorFile, float error[]) {

    int i = 0;
    int size = getSizeOfLayer();

    for (i = 0; i < size && 1==fscanf(inputErrorFile, "%f", &error[i]); i++); 

    return i;
}
/**
 * @param layer a layer of data to print
 */
void printLayerOfFloats(float layer[]){

    int x = 0;
    int y = 0;
    int offset = 0;

    fprintf(stderr,"printLayerofFloats() \n");
    fprintf(stdout, "[ ");

    /* FOR EACH COLUMN */
    for (y = 0; y < getHeightOfColumns(); y++) {

	offset = y * getLengthOfRows();

	/* FOR EACH ROW */
	for (x = 0; x < getLengthOfRows(); x++) {

	    fprintf(stdout, "%f ", layer[offset + x]);

	} /* end FOR EACH ROW */
	fprintf(stdout, "\n");
    } /* end FOR EACH COLUMN */
    fprintf(stdout, " ]\n");
}