/**
 * 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()  */

#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);
int  outputLayerOfCubes(float layerFront[], float layerBack[], 
			 int currentLayer);
int  outputRowOfCubes(float frontLayerOfData[], float  backLayerOfDate[], 
		       int layerOffset, int currentLayer);
int  readLayerOfFloats(FILE *inputFile, float layer[]);
int  copyLayer(float sourceLayer[], float destinationLater[]);
void printLayerOfFloats(float layer[]);
void printCube(struct cube *c_ptr);
int  whichOctant(int xCoord, int yCoord, int zCoord);

/**
 * 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 fromFloatToCube.c -Iinclude -lm -o fromFloatToCube.exe
 * </PRE>
 * <P>
 * Start the program using:
 * <PRE>
 * % fromFloatToCube.exe [input file] [x] [y] [z] [octant] > [output file] 
 *   e.g.
 * % fromFloatToCube.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;
    char *inputFileName;
    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);
    }

    inputFileName = argv[1];    

    if ( !( inputFile = fopen(inputFileName,"rb"))){
        fprintf(stderr, "***Error, couldn't open input file: %s \n", 
		         inputFileName);
        exit(1);
    }
    
    setLengthOfRows(   atoi(argv[2]));
    setHeightOfColumns(atoi(argv[3]));
    setDepthOfLayers(  atoi(argv[4]));    
    setOutputOctant(   atoi(argv[5]));    

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

    totalOutput = outputVolumeOfCubes(inputFile);

    if (fclose(inputFile)) {
	fprintf(stderr, "***Error fromFloatToCube.main(): closing input "
		"file \n");
    } else {
	fprintf(stderr, "# fromFloatToCube.main(): closed input file \n");
    }
    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) {

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

    if (debug == 1) 
	fprintf(stdout, "# 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);
	    totalRead = totalRead + readLayerOfFloats(inputFile, layerOfData2);

	} else {

	    copyLayer(layerOfData2, layerOfData1);
	    totalRead = totalRead + readLayerOfFloats(inputFile, layerOfData2);
	} /* end if (z == 0) */

	/** output every layer */
	totalOutput = totalOutput + outputLayerOfCubes(layerOfData1, 
						       layerOfData2, 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[], 
		       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, rowOffset,
			     currentLayer);

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

/**
 * This method outputs the 8 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[], 
		     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", 
   /** getLengthOfRows() == moving up 1 column */
   /* 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() ]);

	    fprintf(stdout, "%d %d %d\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 layer a layer of data to print
 */
void printLayerOfFloats(float layer[]) {

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

    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");
}
/**
 * @param c  the cube to print out
 */
void printCube(struct cube *c) {

    fprintf(stdout, "%f %f %f %f %f %f %f %f\n", 
	    c->p0, c->p1, c->p2, c->p3, c->p4, c->p5, c->p6, c->p7);
    fprintf(stdout, "%d %d %d \n", 
	    c->xIndex, c->yIndex, c->zIndex);
}