/*
	File: MC33_util_grd.c
	Programed by: David Vega: dvega@uc.edu.ve
	March 2012
*/
#ifndef MC33_util_grd_c
#define MC33_util_grd_c
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#ifndef GRD_data_type
#define GRD_data_type float
#endif

#ifndef LAGRANGE3D4GRD_H
//*
//struct that store a tridimensional grid
typedef struct
{
	GRD_data_type ***F;
	int N[3];
	float L[3], Ang[3];//not necesary
	double r0[3], d[3];
#ifndef _ORTHO_GRD
	int nonortho; //necessary if _ORTHO_GRD is not defined
	double _A[3][3], A_[3][3]; //necessary if _ORTHO_GRD is not defined
#endif
	int periodic;//not necesary
	char title[128];//not necesary
} _GRD;

char stringaux[128];
#endif
//*/
//******************************************************************
//c = Ab, A is a 3x3 upper triangular matrix. If t != 0 a is transposed.
#ifndef _ORTHO_GRD
#ifndef ULTIMAT_H
void _multTSA_bf(double A[3][3], float *b, float *c, int t = 0)
{
	if(t)
	{
		c[2] = A[0][2]*b[0] + A[1][2]*b[1] + A[2][2]*b[2];
		c[1] = A[0][1]*b[0] + A[1][1]*b[1];
		c[0] = A[0][0]*b[0];
	}
	else
	{
		c[0] = A[0][0]*b[0] + A[0][1]*b[1] + A[0][2]*b[2];
		c[1] = A[1][1]*b[1] + A[1][2]*b[2];
		c[2] = A[2][2]*b[2];
	}
}

void (*mult_Abf)(double (*)[3], float *, float *, int);
mult_Abf = _multTSA_bf;
#endif
#endif

//******************************************************************
/*
lagrang3D calculates the function interpolated value (f), evaluated at the point
r[3], using a 4x4x4 points three-dimensional subgrid of the grid *Z.
*/
inline int _L3D_interpolationf(_GRD *Z, float *r, float &f)
{
	int iR[3], x, y, z, i, j, k;
	double t, w[3][4];
	f = 0;
	for(i=0;i<3;i++)
	{
		if((Z->periodic>>i)&0x01)	
		{
			t = (r[i] - Z->L[i]*floor(r[i]/Z->L[i]))/Z->d[i];
			iR[i] = (int)t;
			t -= iR[i];
			iR[i] = (iR[i]? iR[i] - 1: Z->N[i] - 1);
		}
		else
		{
			t = r[i]/Z->d[i];
			if(t < 0 || t > Z->N[i])
				return -1;
			iR[i] = (int)t;
			t -= iR[i];
			if(!iR[i])
				t--;
			else if(iR[i] > Z->N[i] - 2)
			{
				z = iR[i] - Z->N[i] + 2;
				iR[i]-= z + 1;
				t+= z;
			}
			else
				iR[i]--;
		}
		w[i][0] = (-2.0 + (3.0 - t)*t)*t*(1.0/6.0);
		w[i][1] = (2.0 + (-1.0 + (-2.0 + t)*t)*t)*0.5;
		w[i][2] = (2.0 +(1.0 - t)*t)*t*0.5;
		w[i][3] = (t*t - 1.0)*t*(1.0/6.0);
	}
	z = iR[2];
	for(k=0;k<4;k++)
	{
		y = iR[1];
		for(j=0;j<4;j++)
		{
			x = iR[0];
			for(i=0;i<4;i++)
			{
				f += w[0][i]*w[1][j]*w[2][k]*Z->F[z][y][x];
				if ((++x) > Z->N[0]) x = 1;
			}
			if ((++y) > Z->N[1]) y = 1;
		}
		if ((++z) > Z->N[2]) z = 1;
	}
	return 0;
}

int lagrang3Df(_GRD *Z, float *r, float &f)
{
	float r0[3];
	int i;
	for(i=0;i<3;i++)
		r0[i] = r[i] - Z->r0[i];
#ifndef _ORTHO_GRD
	if(Z->nonortho)
		mult_Abf(Z->A_,r0,r0,0);
#endif
	return _L3D_interpolationf(Z,r0,f);
}


//******************************************************************
#ifndef LAGRANGE3D4GRD_H
void free_memory_grd(_GRD* &Z)
{
	int k, j;
	if(Z)
	{
		if(Z->F)
		{
			for(k=0;k<=Z->N[2];k++)
			{
				if(Z->F[k])
					for(j=0;j<=Z->N[1];j++)
						free(Z->F[k][j]);
				free(Z->F[k]);
			}
			free(Z->F);
		}
		free(Z);
		Z = 0;
	}
}

//******************************************************************
/*
read_grd read a filename file (the file must be a output *.grd file from the
DMol program), it returns a pointer to struct _GRD that contains all the grid
data.
*/
_GRD* read_grd(char *filename)
{
	_GRD *Z;
	FILE *in;
	int  i, j, k;
	double ca, cb, sg, cg, aux1, aux2, Arad[3];
	int grd_xi[3], grd_xf[3], grd_ordenij;
	float Ang[3];

	Z = (_GRD*)malloc(sizeof(_GRD));
	if(!Z) return 0;

#ifndef _ORTHO_GRD
	Z->nonortho = 0;
#endif
	Z->periodic = 0;

	in = fopen(filename,"r");
	if (!in) return 0;
	fgets(stringaux,100,in);
	fgets(stringaux,60,in);
	fgets(stringaux,60,in);
	sscanf(stringaux,"%f %f %f %f %f %f",&(Z->L[0]),&(Z->L[1]),&(Z->L[2]),&Ang[0],&Ang[1],&Ang[2]);
	fgets(stringaux,60,in);
	sscanf(stringaux,"%d %d %d",&(Z->N[0]),&(Z->N[1]),&(Z->N[2]));
	fgets(stringaux,60,in);
	sscanf(stringaux,"%d %d %d %d %d %d %d",&grd_ordenij,&grd_xi[0],&grd_xf[0],&grd_xi[1],
						&grd_xf[1],&grd_xi[2],&grd_xf[2]);
	if (Z->N[0] < 2 || Z->N[1] < 2 || Z->N[2] < 2) return 0;
	if (grd_ordenij != 1 && grd_ordenij != 3) return 0;
	for (i=0;i<3;i++)
	{
		Z->d[i] = Z->L[i]/Z->N[i];
		Z->r0[i] = grd_xi[i]*Z->d[i];
	}

	Z->periodic = (grd_xi[0] == 0)|((grd_xi[1] == 0)<<1)|((grd_xi[2] == 0)<<2);

#ifndef _ORTHO_GRD
	if(Ang[0] != 90 || Ang[1] != 90 || Ang[2] != 90)
	{
		Z->nonortho = 1;
		for (i=0;i<3;i++)
			Arad[i] = Ang[i]*(M_PI/180.0);
		ca = cos(Arad[0]);
		cb = cos(Arad[1]);
		sg = sin(Arad[2]);
		cg = cos(Arad[2]);
		aux1 = ca - cb*cg;
		aux2 = sqrt(sg*sg + 2*ca*cb*cg - ca*ca - cb*cb);
		Z->_A[0][0] = Z->A_[0][0] = 1.0;
		Z->_A[0][1] = cg;
		Z->_A[0][2] = cb;
		Z->_A[1][1] = sg;
		Z->A_[1][1] = cb = 1.0/sg;
		Z->A_[0][1] = -cg*cb;
		Z->_A[1][2] = aux1*cb;
		Z->_A[2][2] = aux2*cb;
		aux2 = 1.0/aux2;
		Z->A_[0][2] = (cg*aux1 - ca*sg*sg)*cb*aux2;
		Z->A_[1][2] = -aux1*cb*aux2;
		Z->A_[2][2] = sg*aux2;
	}
	else
	{
		Z->nonortho = 0;
		for (i=0;i<3;i++)
			Z->_A[i][i] = Z->A_[i][i] = 1.0;
		Z->_A[0][1] = Z->_A[0][2] = Z->_A[1][2] = 0;
		Z->A_[0][1] = Z->A_[0][2] = Z->A_[1][2] = 0;
	}
	Z->_A[1][0] = Z->_A[2][0] = Z->_A[2][1] = 0;
	Z->A_[1][0] = Z->A_[2][0] = Z->A_[2][1] = 0;
#endif

	Z->F = (GRD_data_type***)malloc((Z->N[2]+1)*sizeof(void*));
	if(!Z->F)
	{
		free(Z);
		return 0;
	}
	for(k=0;k<=Z->N[2];k++)
	{
		Z->F[k] = (GRD_data_type**)malloc((Z->N[1]+1)*sizeof(void*));
		if(!Z->F[k])
		{
			free_memory_grd(Z);
			return 0;
		}
		for(j=0;j<=Z->N[1];j++)
			Z->F[k][j] = (GRD_data_type*)malloc((Z->N[0]+1)*sizeof(GRD_data_type));
	}
	if(!Z->F[Z->N[2]][Z->N[1]])
	{
		free_memory_grd(Z);
		return 0;
	}
	for(k=0;k<=Z->N[2];k++)
		if (grd_ordenij == 1)
		{
			for(j=0;j<=Z->N[1];j++)
				for(i=0;i<=Z->N[0];i++)
				  fscanf(in,"%f",&Z->F[k][j][i]);
		}
		else
		{
			for(i=0;i<=Z->N[0];i++)
				for(j=0;j<=Z->N[1];j++)
					fscanf(in,"%f",&Z->F[k][j][i]);
		}
	fclose(in);
	return Z;
}

/*
internal format
*/
_GRD* read_grd_binary(char* filename)
{
	FILE* in;
	int i, j, k;
	char title[131];
	_GRD* Z;
	//char identifier[5] = "_GRD";
#define grb_identifier 0x4452475f

	in = fopen(filename,"rb");
	if(!in)
		return 0;
	fread(&i,sizeof(int),1,in);
//	if(i != *(int*)identifier)
	if(i != grb_identifier)
		return 0;

	fread(&i,sizeof(int),1,in);
	Z = (_GRD*)malloc(sizeof(_GRD));
	if(!Z) return 0;
	fread(title,i*sizeof(char),1,in);
//	title[i] = '\0';
	fread(Z->N,sizeof Z->N,1,in);
	fread(Z->L,sizeof Z->L,1,in);
	fread(Z->r0,sizeof Z->r0,1,in);
	fread(Z->d,sizeof Z->d,1,in);
#ifndef _ORTHO_GRD
	fread(&Z->nonortho,sizeof(int),1,in);
	if(Z->nonortho)
	{
		fread(title,3*sizeof(float),1,in);
		fread(Z->_A,sizeof Z->_A,1,in);
		fread(Z->A_,sizeof Z->A_,1,in);
	}
	else
	{
		memset(Z->_A,0,sizeof Z->_A);
		memset(Z->A_,0,sizeof Z->A_);
		for(j=0;j<3;j++)
			Z->_A[j][j] = Z->A_[j][j] = 1;
	}
#endif
	//Z->periodic = 0;
	if(Z->r0[0] == 0 && Z->r0[1] == 0 && Z->r0[2] == 0)
		Z->periodic = 1;
	Z->F = (GRD_data_type***)malloc((Z->N[2] + 1)*sizeof(void*));
	if(!Z->F)
	{
		free(Z);
		return 0;
	}
	for(k=0;k<=Z->N[2];k++)
	{
		Z->F[k] = (GRD_data_type**)malloc((Z->N[1] + 1)*sizeof(void*));
		if(!Z->F[k])
		{
			free_memory_grd(Z);
			return 0;
		}
		for(j=0;j<=Z->N[1];j++)
			Z->F[k][j] = (GRD_data_type*)malloc((Z->N[0] + 1)*sizeof(GRD_data_type));
	}
	if(!Z->F[Z->N[2]][Z->N[1]])
	{
		free_memory_grd(Z);
		return 0;
	}
	for(k=0;k<=Z->N[2];k++)
		for(j=0;j<=Z->N[1];j++)
			fread(Z->F[k][j],(Z->N[0] + 1)*sizeof(GRD_data_type),1,in);
	fclose(in);
	return Z;
}
#endif
//******************************************************************
/*
Read a set of files that contain a slab of res*res scan data points, the data
points are read as unsigned short int (if order is different from 0, the bytes
of the unsigned short are exchanged). The filename must end with a number, and
the fuction read all files with end number greater or equal to filename.
(Some datasets: http://www-graphics.stanford.edu/data/voldata/voldata.html)
*/
_GRD* read_scanfiles(char *filename, int res, int order = 0)
{
	_GRD *Z;
	FILE *in;
	char *nm, *nm2;
	int  i, j, k = 0, l;
	unsigned short int n;
	GRD_data_type ***pt, **p;

	Z = (_GRD*)malloc(sizeof(_GRD));
	if(!Z) return 0;

#ifndef _ORTHO_GRD
	Z->nonortho = 0;
#endif
	Z->periodic = 0;
	memset(Z->r0,0,sizeof Z->r0);
	for(i=0;i<2;i++)
	{
		Z->d[i] = 1;
		Z->L[i] = Z->N[i] = res - 1;
	}
	Z->d[2] = 1;
	Z->F = 0;
	l = strlen(filename) - 1;
	nm = (char*)malloc((l + 5)*sizeof(char));
	if(!nm)
		return 0;
	while(filename[l] >= '0' && filename[l] <= '9')
		l--;
	strncpy(nm,filename,l + 1);
	nm2 = nm + l + 1;
	l = atoi(filename + l + 1);
	while(1)
	{
		sprintf(nm2,"%-d",l++);
		in = fopen(nm,"rb");
		if(!in)
			break;
		if(!(k&63))
		{
			pt = (GRD_data_type***)realloc(Z->F,(k + 64)*sizeof(void*));
			if(!pt)
				break;
			Z->F = pt;
		}
		Z->F[k] = (GRD_data_type**)malloc(res*sizeof(void*));
		if(!Z->F[k])
			break;
		for(j=0;j<res;j++)
			Z->F[k][j] = (GRD_data_type*)malloc(res*sizeof(GRD_data_type));
		if(!Z->F[k][Z->N[1]])
			break;

		if(order)
			for(j=0;j<res;j++)
				for(i=0;i<res;i++)
				{
					fread(&n,sizeof(short int),1,in);
					Z->F[k][j][i] = (unsigned short int)((n>>8)|(n<<8));
				}
		else
			for(j=0;j<res;j++)
				for(i=0;i<res;i++)
				{
					fread(&n,sizeof(short int),1,in);
					Z->F[k][j][i] = n;
				}
		fclose(in);
		k++;
	}
	free(nm);
	Z->L[2] = Z->N[2] = --k;
	for(i=(k-1)/2;i>=0;i--)
	{
		p = Z->F[i];
		Z->F[i] = Z->F[k - i];
		Z->F[k - i] = p;
	}
#ifndef _ORTHO_GRD
	memset(Z->_A,0,sizeof Z->_A);
	memset(Z->A_,0,sizeof Z->A_);
	for(j=0;j<3;j++)
		Z->_A[j][j] = Z->A_[j][j] = 1;
#endif
	return Z;
}

//******************************************************************
/*
Read a file that contains only the data points. byte plus 1 is the number of
bytes of each data point (1 to 4) and the bit order (little or big endian, if
it is big endian, byte is negative) . The vector N[3] must contains the number
of points in each dimension, and the size of file must be byte*N[0]*N[1]*N[2].
*/
_GRD* read_raw_file(char *filename, int *N, int byte = 0)
{
	_GRD *Z;
	FILE *in;
	unsigned int ui = 0;
	int  i, j, k;
	if(byte < -3 || byte > 3)
		return 0;
	Z = (_GRD*)malloc(sizeof(_GRD));
	if(!Z)
		return 0;

#ifndef _ORTHO_GRD
	Z->nonortho = 0;
#endif
	Z->periodic = 0;
	memset(Z->r0,0,sizeof Z->r0);
	for(i=0;i<3;i++)
		Z->d[i] = 1;
	in = fopen(filename,"rb");
	if(!in)
		return 0;
	Z->L[0] = Z->N[0] = N[0] - 1;
	Z->L[1] = Z->N[1] = N[1] - 1;
	Z->L[2] = Z->N[2] = N[2] - 1;
	Z->F = (GRD_data_type***)malloc(N[2]*sizeof(void*));
	if(!Z->F)
	{
		free(Z);
		return 0;
	}
	for(k=0;k<N[2];k++)
	{
		Z->F[k] = (GRD_data_type**)malloc(N[1]*sizeof(void*));
		if(!Z->F[k])
		{
			free_memory_grd(Z);
			return 0;
		}
		for(j=0;j<N[1];j++)
			Z->F[k][j] = (GRD_data_type*)malloc(N[0]*sizeof(GRD_data_type));
	}
	if(!Z->F[Z->N[2]][Z->N[1]])
	{
		free_memory_grd(Z);
		return 0;
	}
	if (byte < 0)
	{
		byte = 1 - byte;
		for(k=N[2]-1;k>=0;k--)
			for(j=0;j<N[1];j++)
				for(i=0;i<N[0];i++)
				{
					fread(&ui,byte,1,in);
					if (byte == 2)
						Z->F[k][j][i] = (unsigned short int)((ui>>8)|(ui<<8));
					else if (byte == 4)
						Z->F[k][j][i] = (ui>>24)|((ui>>8)&0x00f0)|((ui<<8)&0x0f00)|(ui<<24);
					else
						Z->F[k][j][i] = (ui>>16)|(ui&0x00f0)|(ui<<16);
				}
	}
	else
	{
		byte++;
		for(k=N[2]-1;k>=0;k--)
			for(j=0;j<N[1];j++)
				for(i=0;i<N[0];i++)
				{
					fread(&ui,byte,1,in);
					Z->F[k][j][i] = ui;
				}
	}
	fclose(in);
#ifndef _ORTHO_GRD
	memset(Z->_A,0,sizeof Z->_A);
	memset(Z->A_,0,sizeof Z->A_);
	for(j=0;j<3;j++)
		Z->_A[j][j] = Z->A_[j][j] = 1;
#endif
	return Z;
}

//******************************************************************
/*
read a dat file:
http://www.cg.tuwien.ac.at/research/vis/datasets/
*/
_GRD* read_dat_file(char *filename)
{
	_GRD *Z;
	FILE *in;
	int  i, j, k;
	unsigned short int n, nx, ny, nz;

	Z = (_GRD*)malloc(sizeof(_GRD));
	if(!Z)
		return 0;

#ifndef _ORTHO_GRD
	Z->nonortho = 0;
#endif
	Z->periodic = 0;
	memset(Z->r0,0,sizeof Z->r0);
	for(i=0;i<3;i++)
		Z->d[i] = 1;
	in = fopen(filename,"rb");
	if(!in)
		return 0;
	fread(&nx,sizeof(short int),1,in);
	fread(&ny,sizeof(short int),1,in);
	fread(&nz,sizeof(short int),1,in);
	Z->L[0] = Z->N[0] = nx - 1;
	Z->L[1] = Z->N[1] = ny - 1;
	Z->L[2] = Z->N[2] = nz - 1;

	Z->F = (GRD_data_type***)malloc(nz*sizeof(void*));
	if(!Z->F)
	{
		free(Z);
		return 0;
	}
	for(k=0;k<nz;k++)
	{
		Z->F[k] = (GRD_data_type**)malloc(ny*sizeof(void*));
		if(!Z->F[k])
		{
			free_memory_grd(Z);
			return 0;
		}
		for(j=0;j<ny;j++)
			Z->F[k][j] = (GRD_data_type*)malloc(nx*sizeof(GRD_data_type));
	}
	if(!Z->F[Z->N[2]][Z->N[1]])
	{
		free_memory_grd(Z);
		return 0;
	}
	for(k=nz-1;k>=0;k--)
		for(j=0;j<ny;j++)
			for(i=0;i<nx;i++)
			{
				fread(&n,sizeof(short int),1,in);
				Z->F[k][j][i] = n;
			}
	fclose(in);
#ifndef _ORTHO_GRD
	memset(Z->_A,0,sizeof Z->_A);
	memset(Z->A_,0,sizeof Z->A_);
	for(j=0;j<3;j++)
		Z->_A[j][j] = Z->A_[j][j] = 1;
#endif
	return Z;
}

//******************************************************************
#endif
