/**
 * @file        step1.c
 * @author      Robert S Laramee
 * @version     1.0
 * @see         step0.c
 * <P>
 * start date   2 July 1999
 * finish date  7 July 1999
 * <P>
 * Description  This is step 1 of 0-2 steps.  Step 1 is to take the output
 *              from step 0 and compute a (wavelet) transformation in
 * 2 dimensions, averaging blocks of 4 (as opposed to 8 in 3D) to reduce
 * x-256, y-256, z-128 slices to x-128, y-128, z-128 slices.
 * <P>
 * This is an adaptation and simplification of Phil and Mike's build_mr.c 
 * program.
 *
 * Some handy numbers:
 * 129 x 129 x 129 = 2,146,689
 * 128 x 128 x 128 = 2,097,152 (this is the one we want as output)
 *  65 x  65 x  65 =   274,625
 *  64 x  64 x  64 =   262,144
 *  32 x  32 x  32 =    32,768
 *  16 x  16 x  16 =     4,096
 */

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

#define XHEADROWS    256
#define YHEADCOLUMNS 256
#define ZHEADLAYERS  128

#define XKNEEROWS    256
#define YKNEECOLUMNS 256
#define ZKNEELAYERS  128

/* the lobster data does not need this procedure
 * #define XLOBSTERROWS    128
 * #define YLOBSTERCOLUMNS 128
 * #define ZLOBSTERLAYERS  64
 */

/**
 * Compile the program using:
 * % gcc -Wall step1.c -Iinclude -lm -o step1.exe 
 * Start the program using:
 * % step1 [input file] > [output file]
 * @param argv[1] -the input data file name
 */
int main(int argc , char *argv[]) {

    int    index   = 0;
    int    offset  = 0;
    int    x,y,z;

    float  average = 0;
    float *dataArray_ptr, *index_ptr;
    FILE  *inputFile;

    /*  check command line arguments  */
    if (argc != 2) {
	fprintf(stderr,"***Error, usage: step1 [inputfile] > [output file]\n");
	exit(1);
    }
    /* don't foget to close the file */
    if ( !( inputFile = fopen(argv[1], "rb"))){
	fprintf(stderr, "***Error, couldn't open input file: %s\n", argv[1]);
	exit(1);
    }
    fprintf(stderr, "# main(): reading from file: %s...\n", argv[1]);

    /**
     * Memory allocation
     * 1.  an array of floats for each number in the volume data set
     */
    if ( !(dataArray_ptr = (float *)malloc(XKNEEROWS    *
					   YKNEECOLUMNS *
					   ZKNEELAYERS  *
					   sizeof(float))
	)) {
	fprintf(stderr, "*** Error: couldn't allocate space for data. \n");
	exit(1);
    }

    /** 
     * This procedure reads in the data and stores it in the data array.
     */
    for (z = 0; z < ZKNEELAYERS; z++) {
	for (y = 0; y < YKNEECOLUMNS; y++) {
	    for (x = 0; x < XKNEEROWS; x++) {

		index_ptr = dataArray_ptr + index;
		fscanf(inputFile, "%f", index_ptr);
		index++;
	    }
	}
    }

    index = 0;
    /**
     * Here is where we will do our transformation:
     * FOR EACH LAYER
     *    We'll reduce the size of each layer from 256x256 to 128x128
     *    by doing an averaging in 2 dimensions.  Average 4 neighboring
     *    xy points to 1 xy point.
     */
    /* FOR EVERY ZLAYER */
    for (z = 0; z < ZKNEELAYERS; z++) {
	/*  fprintf(stdout, "begin layer z = %d\n", z);  */
	
	/* FOR EVERY OTHER YCOLUMN  */
	for (y = 0; y < YKNEECOLUMNS; y = y+2) {
	    /*   fprintf(stdout, "begin column y = %d\n", y);  */

	    /* FOR EVERY OTHER XROW  */
	    for (x = 0; x < XKNEEROWS; x = x+2)  {

		offset = (z * XKNEEROWS * YKNEECOLUMNS) + 
		         (y * YKNEECOLUMNS) + x;

		average = (dataArray_ptr[offset                      ] + 
			   dataArray_ptr[offset + 1                  ] +
			   dataArray_ptr[offset + YKNEECOLUMNS       ] +
			   dataArray_ptr[offset + YKNEECOLUMNS + 1   ])/4;

		fprintf(stdout, "%f\n", average);
		index++;
	    }  /*  end FOR EVERY OTHER XROW  */
	}  /*  end FOR EVERY OTHER YCOL  */
    }  /*  end FOR EVERY XLAYER  */

    free(dataArray_ptr);  /*  release memory back to heap  */
    fprintf(stderr, "# number of lines is %d\n", index);
    /*  fprintf(stdout, "# number of lines is %d\n", index);  */
    return 0;
}  /*  end main()  */