/* kernel_readin.c */
#include "mex.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

// #define batchfilename "beamletbatch0.bin"

#define BATCHNAME       (prhs[0])
#define BEAMLETS        (plhs[0])   // output structure

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	FILE *fid;

	mxArray *field, *struc;

	int j,k;
	int celldims[2];  // dimensions in the output cell array for the beamlets
	int dims[3];      // dimensions of sparse grids
	int dims2[2];     // dimensions of arrays that will hold indices and data
	int Nind;         // number of non-zero indices
	int Nbeamlets;    // number of beamlets in file
	int num;          // beamlet number
	int *intPtr;
	float *floatPtr;

	char errstr[200];

	const int one = 1;

	char *batchfilename;  // beamlet batch file name

	// information for sparse output data:
	int nfields = 6;
	const char *field_names[6] = {"num","x_count","y_count","z_count",
							"non_zero_indices","non_zero_values"};

	batchfilename = mxArrayToString(BATCHNAME);

	// open the beamlet batch file
	if ((fid = fopen(batchfilename,"rb")) == NULL)
	{
		sprintf(errstr,"Failed to open beamlet batch file %s\n",batchfilename);
		mexErrMsgTxt(errstr);
	}

	// read in the number of beamlets that are in the batch file
	if (fread(&Nbeamlets,sizeof(int),1,fid) != 1)
	{
		sprintf(errstr,"Error reading the number of beamlets from %s\n",batchfilename);
		mexErrMsgTxt(errstr);
	}

	// declare the BEAMLETS structure
	celldims[0] = 1;
	celldims[1] = Nbeamlets;

	if ((BEAMLETS = mxCreateCellArray(2,celldims)) == NULL)
		mexErrMsgTxt("Error in the declaration of the beamlets cell array.\n");

	// read in each of the beamlets from the file
	for (k=0;k<Nbeamlets;k++)
	{
		// create the structure that will hold the sparse beamlet
		if ((struc = mxCreateStructArray(1,&one,nfields,field_names)) == NULL)
		{
			sprintf(errstr,"Failed to allocate memory for the structure that will hold sparse beamlet %d\n",k);
			mexErrMsgTxt(errstr);
		}

		// read in the beamlet number
		if (fread(&num,sizeof(int),1,fid) != 1)
		{
			sprintf(errstr,"Failed to read in the beamlet number from beamlet %d.\n",k);
			mexErrMsgTxt(errstr);
		}

		// read in the beamlet dimensions
		if (fread(dims,sizeof(int),3,fid) != 3)
		{
			sprintf(errstr,"Failed to read in the dimensions of beamlet %d.\n",k);
			mexErrMsgTxt(errstr);
		}

		// assign fields to the structure for the beamlet number and beamlet dimensions
		dims2[0] = 1;
		dims2[1] = 1;

		if ((field = mxCreateNumericArray(2,dims2,mxINT32_CLASS,mxREAL)) == NULL)
		{
			sprintf(errstr,"Failed to allocate memory for the dummy field for beamlet %d\n",k);
			mexErrMsgTxt(errstr);
		}

		intPtr = (int *)mxGetData(field);

		intPtr[0] = num;
		mxSetField(struc,0,"num",mxDuplicateArray(field));

		intPtr[0] = dims[0];
		mxSetField(struc,0,"x_count",mxDuplicateArray(field));

		intPtr[0] = dims[1];
		mxSetField(struc,0,"y_count",mxDuplicateArray(field));

		intPtr[0] = dims[2];
		mxSetField(struc,0,"z_count",mxDuplicateArray(field));

		// read in the number of non-zero indices
		if (fread(&Nind,sizeof(int),1,fid) != 1)
		{
			sprintf(errstr,"Failed to read in the number of non-zero indices for beamlet k.\n",k);
			mexErrMsgTxt(errstr);
		}

		// assign fields to the structure for the beamlet dimensions
		dims2[0] = 1;
		dims2[1] = Nind;

		// allocate memory for non-zero indices and read-in the data
		if ((field = mxCreateNumericArray(2,dims2,mxINT32_CLASS,mxREAL)) == NULL)
		{
			sprintf(errstr,"Failed to allocate memory for the index array of beamlet %d\n",k);
			mexErrMsgTxt(errstr);
		}

		intPtr = (int *)mxGetData(field);
		if (fread(intPtr,sizeof(int),Nind,fid) != Nind)
		{
			sprintf(errstr,"Unable to read-in non-zero indices for beamlet %d.",k);
			mexErrMsgTxt(errstr);
		}

		// add unity to all of the indices to convert them to Matlab indexing notation
		for (j=0;j<Nind;j++) intPtr[j] += 1;

		mxSetField(struc,0,"non_zero_indices",mxDuplicateArray(field));

		// read in the non-zero values
		if ((field = mxCreateNumericArray(2,dims2,mxSINGLE_CLASS,mxREAL)) == NULL)
		{
			sprintf(errstr,"Failed to allocate memory for the data array of beamlet %d\n",k);
			mexErrMsgTxt(errstr);
		}

		floatPtr = (float *)mxGetData(field);
		if (fread(floatPtr,sizeof(float),Nind,fid) != Nind)
		{
			sprintf(errstr,"Unable to read-in non-zero data for beamlet %d.",k);
			mexErrMsgTxt(errstr);
		}
		mxSetField(struc,0,"non_zero_values",mxDuplicateArray(field));
		
		// copy the sparse beamlet structure to the beamlet batch cell array
		mxSetCell(BEAMLETS,k,mxDuplicateArray(struc));
	}

	fclose(fid);

	// I have not managed to free the memory allocated for the field and struc arrays
	// due to segmentation errors.  I'm not sure how to get around this problem either.
}