#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include "vcbutils.h"
#include "vcbmcube.h"
#include <mpi.h>

int main(int argc, char ** argv)
{
	vcbdatatype datatype;
	int ndims, orig[3], dsize[3], nattribs;
	int cutsize[3], lbounds[3], ubounds[3];
	int ncuts, k, interval, headlen, i;
	void * outdata;
	FILE * fp;

	/* application specific data starts here */
	int   nverts, nfacets, * fdata;
	float isoval, * vdata, * fnormals;
	/* end of application specific data */

	/* See if user passed an argument. */
	if (argc != 2) {
		fprintf(stderr, "Usage: %s datafile\n", argv[0]);
		exit(1);
	}

	MPI_Init(&argc,&argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &k);
	MPI_Comm_size(MPI_COMM_WORLD, &ncuts);

	vcbReadHeaderBinm(argv[1], &datatype, &ndims, orig, dsize, &nattribs);

	/* set up partition boundaries */
	interval = (int)ceil(((double)dsize[0])/ncuts);
	lbounds[0] = interval * k;
	ubounds[0] = interval * (k+1); 
	if (ubounds[0] >= dsize[0]) 
	ubounds[0] = dsize[0]-1;
	cutsize[0] = ubounds[0] - lbounds[0] + 1;

	/* set up other bounds */
	lbounds[1] = 0;
	ubounds[1] = dsize[1] - 1;
	lbounds[2] = 0;
	ubounds[2] = dsize[2] - 1;
	cutsize[1] = dsize[1];
	cutsize[2] = dsize[2];

	outdata = malloc(cutsize[0]*cutsize[1]*cutsize[2]*nattribs*vcbSizeOf(datatype));

	fp = fopen(argv[1],"rb");
	headlen = (3+2*ndims)*sizeof(int)+sizeof(datatype);

	vcbFileGrablk(fp, headlen, outdata, 
				nattribs,
				vcbSizeOf(datatype),
				ndims,
				dsize,
				lbounds, 
				ubounds, 
				NULL);

	for (i = 0; i < ndims; orig[i] += lbounds[i], i++);

	isoval = 45.f;
	if ((nfacets = vcbMarchingCubeBlk(datatype, outdata, isoval, cutsize, 
					VCB_WITHNORMAL, &nverts, &vdata, &fdata))<=0) {
		fprintf(stderr,
				"p%d: vcbMarchingCubeBlk: extract zero triangles at isovalue %f\n",
				k, isoval);
	}

	/* saving vertex array and triangle array into endian safe formats */
	for (i = 0; i < nverts; i ++) {
		vdata[i*6+0] += orig[0];
		vdata[i*6+1] += orig[1];
		vdata[i*6+2] += orig[2];
	}

	/* output the surface mesh */
	vcbGenBinm("vertexarray.bin", VCB_FLOAT, 1, NULL, &nverts, 6, vdata);
	vcbGenBinm("triangarray.bin", VCB_INT, 1, NULL, &nfacets, 3, fdata);

	free(vdata);
	free(fdata);
	free(fnormals);
	free(outdata);
	fclose(fp);

	MPI_Finalize();

	return 0;
}