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

#define RADDIM		16
#define LATDIM		19
#define LNGDIM		37

#define MINRAD		  0.
#define MAXRAD		  1.5

#define MINLAT		(-90.)
#define MAXLAT		( 90.)

#define MINLNG		  0.
#define MAXLNG		360.


#define DRAD		( ( MAXRAD - MINRAD ) / (float)( RADDIM - 1 ) )
#define DLAT		( ( MAXLAT - MINLAT ) / (float)( LATDIM - 1 ) )
#define DLNG		( ( MAXLNG - MINLNG ) / (float)( LNGDIM - 1 ) )


#define DEFAULTXYZRES		128


#define DX		( ( 2. * MAXRAD ) / ( float)( XyzRes - 1 ) )
#define DY		( ( 2. * MAXRAD ) / ( float)( XyzRes - 1 ) )
#define DZ		( ( 2. * MAXRAD ) / ( float)( XyzRes - 1 ) )


#define D2R		( M_PI / 180. )
#define R2D		( 180. / M_PI )


#ifndef FALSE
#define FALSE		0
#endif
#ifndef TRUE
#define TRUE		( ! FALSE )
#endif


int	Debug;
float	Density[RADDIM][LATDIM][LNGDIM];
float	Smin, Smax;
int	XyzRes;



void	CartesianToSpherical( float, float, float, float *, float *, float * );
void	SphericalToCartesian( float, float, float, float *, float *, float * );
float	VolumeC( float, float, float, float, float, float );
float	VolumeS( float, float, float, float, float, float );


void
main( int argc, char *argv[] )
{
	char *infile;
	char *outfile;
	FILE *fpin;
	FILE *fpout;
	float d, dmin, dmax;
	float x, y, z;
	float r, lng, lat;
	int i, j, k;
	int n;
	unsigned char ucd;
	int irad, ilat, ilng;
	float rad0, lat0, lng0;
	float rad1, lat1, lng1;
	float x0, y0, z0;
	float x1, y1, z1;
	float b000, b001, b010, b011;
	float b100, b101, b110, b111;
	float total;
	int setsmin, setsmax;


	/* parse the command line:					*/

	fpin = stdin;
	Debug = FALSE;
	setsmin = setsmax = FALSE;
	for( i = 1, n = 0; i < argc; i++ )
	{
		if( argv[i][0] == '-' )
		{
			switch( argv[i][1] )
			{
				case 'D':
					Debug = TRUE;
					break;

				case 'r':
				case 'R':
					XyzRes = atoi( argv[i+1] );
					fprintf( stderr, "Cartesian resolution = %d\n", XyzRes );
					i++;
					break;

				case 's':
				case 'S':
					switch( argv[i][4] )
					{
						case 'n':
						case 'N':
							setsmin = TRUE;
							Smin = atof( argv[i+1] );
							fprintf( stderr, "Smin = %f\n", Smin );
							i++;
							break;

						case 'x':
						case 'X':
							setsmax = TRUE;
							Smax = atof( argv[i+1] );
							fprintf( stderr, "Smax = %f\n", Smax );
							i++;
							break;

						default:
							fprintf( stderr, "Unknown command line argument: '%s'\n", argv[i] );
					}
					break;

				default:
					fprintf( stderr, "Unknown command line argument: '%s'\n", argv[i] );
					fprintf( stderr, "Usage: %s [-r res] [-smin s] [-smax s] infile outfile\n", argv[0] );
			}
		}
		else
		{
			if( n == 0 )
				infile = argv[i];
			else
				outfile = argv[i];
			n++;
		}
	}


	if( n == 0 )
	{
		fprintf( stderr, "No input or output file specified!\n" );
		exit( 1 );
	}

	if( n == 1 )
	{
		fprintf( stderr, "No output file specified!\n" );
		exit( 1 );
	}


	fpin = fopen( infile, "r" );
	if( fpin == NULL )
	{
		fprintf( stderr, "Cannot open input file '%s'\n", infile );
		exit( 1 );
	}


	fpout = fopen( outfile, "w" );
	if( fpout == NULL )
	{
		fprintf( stderr, "Cannot create output file '%s'\n", outfile );
		exit( 1 );
	}


	fprintf( fpout, "Vox1999a\n" );
	fprintf( fpout, "VolumeCount 1\n" );
	fprintf( fpout, "##%c\n", 0x0c );

	fprintf( fpout, "##\n" );
	fprintf( fpout, "VolumeSize %3d  %3d  %3d\n", XyzRes, XyzRes, XyzRes );
	fprintf( fpout, "VoxelSize %2d\n", 8 );
	fprintf( fpout, "Endian B\n" );
	fprintf( fpout, "Field 0 (\n" );
	fprintf( fpout, "\tPosition 0\n" );
	fprintf( fpout, "\tSize 8\n" );
	fprintf( fpout, "\tName Data\n" );
	fprintf( fpout, ")\n" );
	fprintf( fpout, "##%c\n", 0x0c );

	dmin =  1.0e+38;
	dmax = -1.0e+38;

	for( i = 0; i < RADDIM; i++ )
	{
		for( j = 0; j < LATDIM; j++ )
		{
			for( k = 0; k < LNGDIM; k++ )
			{
				n = fscanf( fpin, "%f", &d );
				if( n != 1 )
				{
					fprintf( stderr, "Premature EOF at (%d,%d,%d)\n", i, j, k );
					exit( 1 );
				}
				if( d < -4000. )
					d = 0.;

				Density[i][j][k] = d;

				if( d < dmin )
					dmin = d;
				if( d > dmax )
					dmax = d;
			}
		}
	}
	fclose( fpin );

	fprintf( stderr, "Raw data range = [ %e , %e ]\n", dmin, dmax );


	/* do we want to keep this range, or use one from the command line?	*/

	if( ! setsmin )
		Smin = dmin;

	if( ! setsmax )
		Smax = dmax;



	fprintf( stderr, "\n" );
	for( i = 0, x = -MAXRAD; i < XyzRes; i++, x += DX )
	{
		fprintf( stderr, "\b\b\b%3d", i );
		for( j = 0, y = -MAXRAD; j < XyzRes; j++, y += DY )
		{
			for( k = 0, z = -MAXRAD; k < XyzRes; k++, z += DZ )
			{
				CartesianToSpherical( x, y, z, &r, &lat, &lng );

				/* lat, lng are now in degrees so they can correctly	*/
				/* produce indices into the spherical data:		*/


				/* where in the Density array do these spherical values fall?	*/

				irad = (int)( (   r - MINRAD ) / DRAD );
				ilat = (int)( ( lat - MINLAT ) / DLAT );
				ilng = (int)( ( lng - MINLNG ) / DLNG );

				if( irad < 0  ||  irad >= RADDIM-1  ||
				    ilat < 0  ||  ilat >= LATDIM-1  ||
				    ilng < 0  ||  ilng >= LNGDIM-1 )
				{
					d = 0.;
					ucd = (unsigned char) d;
					fwrite( (char *) &ucd, 1, 1, fpout );
					continue;
				}



				/* values at the lower bound of this subvolume:		*/

				rad0 = MINRAD + (float)irad * DRAD;
				lat0 = MINLAT + (float)ilat * DLAT;
				lng0 = MINLNG + (float)ilng * DLNG;


				/* values at the upper bound of this subvolume:		*/

				rad1 = rad0 + DRAD;
				lat1 = lat0 + DLAT;
				lng1 = lng0 + DLNG;


				/* turn the boundaries into cartesian:			*/
				/* lat0, lng0, lat1, lng1 are all in degrees		*/

				SphericalToCartesian( rad0, lat0, lng0, &x0, &y0, &z0 );
				SphericalToCartesian( rad1, lat1, lng1, &x1, &y1, &z1 );


				/* basis coefficients before normalizing:		*/

#ifdef OLDWAY
				b000 = VolumeC( x, y, z, x1, y1, z1 );
				b001 = VolumeC( x, y, z, x1, y0, z0 );
				b010 = VolumeC( x, y, z, x1, y0, z1 );
				b011 = VolumeC( x, y, z, x1, y0, z0 );
				b100 = VolumeC( x, y, z, x0, y1, z1 );
				b101 = VolumeC( x, y, z, x0, y1, z0 );
				b110 = VolumeC( x, y, z, x0, y0, z1 );
				b111 = VolumeC( x, y, z, x0, y0, z0 );
#endif

				/* lat, lng, lat0, lng0, lat1, lng1 are all in degrees	*/

				b000 = VolumeS( r, lat, lng, rad1, lat1, lng1 );
				b001 = VolumeS( r, lat, lng, rad1, lat0, lng0 );
				b010 = VolumeS( r, lat, lng, rad1, lat0, lng1 );
				b011 = VolumeS( r, lat, lng, rad1, lat0, lng0 );
				b100 = VolumeS( r, lat, lng, rad0, lat1, lng1 );
				b101 = VolumeS( r, lat, lng, rad0, lat1, lng0 );
				b110 = VolumeS( r, lat, lng, rad0, lat0, lng1 );
				b111 = VolumeS( r, lat, lng, rad0, lat0, lng0 );


				/* normalize the basis coefficients:			*/

				total = b000 + b001 + b010 + b011 + b100 + b101 + b110 + b111;
				b000 /= total;
				b001 /= total;
				b010 /= total;
				b011 /= total;
				b100 /= total;
				b101 /= total;
				b110 /= total;
				b111 /= total;


				/* use the basis coefficients:				*/

				d  = b000 * Density[irad  ][ilat  ][ilng  ];
				d += b001 * Density[irad  ][ilat  ][ilng+1];
				d += b010 * Density[irad  ][ilat+1][ilng  ];
				d += b011 * Density[irad  ][ilat+1][ilng+1];
				d += b100 * Density[irad+1][ilat  ][ilng  ];
				d += b101 * Density[irad+1][ilat  ][ilng+1];
				d += b110 * Density[irad+1][ilat+1][ilng  ];
				d += b111 * Density[irad+1][ilat+1][ilng+1];

#ifdef BIGDEBUG
fprintf( stderr, "%3d %3d %3d (%3d,%3d,%3d): %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f : %4.1f\n",
	i, j, k, irad, ilat, ilng,
	Density[irad  ][ilat  ][ilng  ], Density[irad  ][ilat  ][ilng+1], Density[irad  ][ilat+1][ilng  ],
	Density[irad  ][ilat+1][ilng+1], Density[irad+1][ilat  ][ilng  ], Density[irad+1][ilat  ][ilng+1],
	Density[irad+1][ilat+1][ilng  ], Density[irad+1][ilat+1][ilng+1], d );
#endif


				/* scale into a volume value:				*/

				d = 255. * ( d - Smin ) / ( Smax - Smin );
				if( d < 0. )
					d = 0.;
				if( d > 255. )
					d = 255.;
				ucd = (unsigned char) d;

				fwrite( (char *) &ucd, 1, 1, fpout );
			}
		}
	}
	fprintf( stderr, "\n" );

	fclose( fpout );
}



/* Note: Latitude and Phi are the same ; Longitude and Theta are the same!	*/
/* angles come in and go out in degrees						*/


void
CartesianToSpherical( float x, float y, float z, float *r, float *phi, float *theta )
{
	*r = sqrt( x*x + y*y + z*z );
#ifdef ORIGINAL_WAY
	*theta = R2D * atan2( x, z );
#endif
	*theta = R2D * atan2( z, x );
	*phi =   R2D * atan2( y, sqrt( x*x + z*z ) );
	if( *theta < 0. )
		*theta += 360.;

	if( *phi < MINLAT  ||  *phi > MAXLAT )
		fprintf( stderr, "Warning: lat = %f\n", *phi );

	if( *theta < MINLNG  ||  *theta > MAXLNG )
		fprintf( stderr, "Warning: lng = %f\n", *theta );
}



/* Note: Latitude and Phi are the same ; Longitude and Theta are the same!	*/
/* angles come in and go out in degrees						*/

void
SphericalToCartesian( float r, float phi, float theta, float *x, float *y, float *z )
{
	float xz;

	phi   *= D2R;
	theta *= D2R;

	*y =  r * sin(phi);
	xz =  r * cos(phi);

#ifdef ORIGINAL_WAY
	*x = xz * sin(theta);
	*z = xz * cos(theta);
#endif
	*x = xz * cos(theta);
	*z = xz * sin(theta);
}


float
VolumeC( float x0, float y0, float z0, float x1, float y1, float z1 )
{
	return fabs( ( x1 - x0 ) * ( y1 - y0 ) * ( z1 - z0 ) );
}


/* angles come in and go out in degrees						*/

float
VolumeS( float r0, float phi0, float theta0, float r1, float phi1, float theta1 )
{
	float fr, ftheta, fphi;

	phi0   *= D2R;
	phi1   *= D2R;
	theta0 *= D2R;
	theta1 *= D2R;

	fr = r1*r1*r1 - r0*r0*r0;
	ftheta = theta1 - theta0;
	fphi = sin(phi1) - sin(phi0);

	return fabs( fr*ftheta*fphi/3. );
}