/*============================================================================
	Title: HDR.cpp
	Module: Pi/Engine
	Author: Ignacio Castaño
	Date: 04/10/2003
============================================================================*/

/*----------------------------------------------------------------------------
	Doc
----------------------------------------------------------------------------*/

/** @file HDR.cpp
 * @brief Hight Dynamic Range environment maps.
**/


/*----------------------------------------------------------------------------
	Includes
----------------------------------------------------------------------------*/

// Core
#include "FileSystem.h"
#include "Parser.h"

// MathLib
#include "SphericHarm.h"
#include "Color.h"

// Engine
#include "HDR.h"


/*----------------------------------------------------------------------------
	Statics
----------------------------------------------------------------------------*/

static int offset_array[6][2] = {
	{2, 1}, {0, 1},	// X faces
	{1, 0}, {1, 2},	// Y faces
	{1, 1}, {1, 3},	// Z faces
};


/*----------------------------------------------------------------------------
	Functions
----------------------------------------------------------------------------*/

/** This formula is from Manne Öhrström's thesis.
 * Take two coordiantes in the range [-1, 1] that define a portion of a
 * cube face and return the area of the projection of that portion on the
 * surface of the sphere.
**/
static float AreaElement( float x, float y ) {
	return atan2(x * y, sqrt(SQ(x) + SQ(y) + 1));
}


/** Get direction from cube texel. @a x and @a y are in the [-1, 1] range. */
static void GetCubeDir( Vec3 * dir, int face, float x, float y ) {
	piDebugCheck( face >= 0 );
	piDebugCheck( face < 6 );

	// Set direction according to face.
	switch( face ) {
		case 0:
			dir->Set( 1, -y, -x );
		break;
		case 1:
			dir->Set( -1, -y, x );
		break;
		case 2:
			dir->Set( x, 1, y );
		break;
		case 3:
			dir->Set( x, -1, -y );
		break;
		case 4:
			dir->Set( x, -y, 1 );
		break;
		case 5:
			dir->Set( -x, -y, -1 );
		break;
	}

	dir->Normalize();
}


/** Sample a cubemap from a spherical function. */
PiTexture SampleCubeMap( const ISphericFunction & f, int size, float exposure ) {
	piCheck( size > 0 );
	piCheck( NextPowerOfTwo(size) == uint(size) );
	piCheck( exposure > 0 );


	PiTexture tex( "sampled_map", PI_TF_CUBEMAP|PI_TF_CLAMP|PI_TF_UPLOAD|PI_TF_HAS_COLOR );
	tex.Allocate( size, size );

	uint * map = tex->mem;

	for( int i = 0; i < 6; i++ ) {
		for( int y = 0; y < size; y++ ) {
			for( int x = 0; x < size; x++ ) {

				float fx = 2 * (float(x)+0.5) / size - 1;
				float fy = 2 * (float(y)+0.5) / size - 1;

				Vec3 dir;
				GetCubeDir( &dir, i, fx, fy );

				FColor color;
				f.Sample( &color, dir );

				color *= exposure;
				//color.Normalize();
				color.Clamp();

				*map++ = color;
			}
		}
	}

	return tex;
}


/*----------------------------------------------------------------------------
	Methods
----------------------------------------------------------------------------*/


/** Open an HDR map. */
bool PiHDRMap::Open( const char * fname ) {

	PiAutoPtr<PiVirtualFile> vfile = GFileSystem.OpenFile( fname );
	if( vfile == NULL ) {
		piDebug( "*** Cannot open hdr file '%s'.\n", fname );
		return false;
	}

	PiParser parser;
	parser.StartParseBuffer( vfile->GetMem(), vfile->GetSize() );

	while( true ) {
		parser.GetToken( true );

		if( !strcmp(parser.token, "FORMAT") ) {
			parser.GetToken( false );
			parser.GetToken( false );
		}
		if( !strcmp(parser.token, "EXPOSURE") ) {
			parser.GetToken( false );
			parser.GetToken( false );
		}
		if( !strcmp(parser.token, "-Y") ) {
			parser.GetToken( false );
			height = atoi(parser.token);
			piCheck( height != 0 );
			parser.GetToken( false );
			piCheck( !strcmp(parser.token, "+X") );
			parser.GetToken( false );
			width = atoi(parser.token);
			piCheck( width != 0 );
			break;
		}
	}

	if( !parser.NextLine() ) {
		piDebug( "*** Unexpected end of file in '%s'.\n", fname );
		return false;
	}

	map = new FColor[width * height];

	// Seek binary data.
	vfile->Seek( parser.GetPos(), SEEK_SET );

	PiRGBE * scanline = new PiRGBE[width];

	for( uint y = 0; y < height; y++ ) {

		if( !ReadColors( scanline, width, vfile ) ) {
			piDebug( "*** Unexpected end of file in '%s'.\n", fname );
			break;
		}

		// Convert scanline to floating point format.
		for( uint x = 0; x < width; x++ ) {
			scanline[x].ToVec( map[y * width + x] );
		}
	}

	delete [] scanline;

	return true;
}

/** Close map. */
void PiHDRMap::Close() {
	width = 0;
	height = 0;
	delete [] map; map = NULL;
}

/** Sample the map. */
void PiHDRMap::Sample( FColor * out, int x, int y ) const {
	piCheck( x >= 0 && x < int(width) );
	piCheck( y >= 0 && y < int(height) );
	*out = map[y * width + x];
}


/** Read a scanline in mixed RLE format. */
bool PiHDRMap::ReadColors( PiRGBE * scanline, int len, PiVirtualFile * vfile ) {

	const int MINELEN = 8;
	const int MAXELEN = 0x7fff;

	// determine scanline type
	if( len < MINELEN || len > MAXELEN ) {
		return OldReadColors( scanline, len, vfile );
	}

	int i = vfile->Read8();
	if( vfile->EndOfFile() ) {
		return false;
	}

	if( i != 2 ) {
		vfile->Seek( -1, SEEK_CUR );
		return OldReadColors( scanline, len, vfile );
	}


	scanline->g = vfile->Read8();
	scanline->b = vfile->Read8();
	scanline->e = vfile->Read8();


	if( scanline->g != 2 && scanline->b & 128 ) {
		scanline->r = 2;
		return OldReadColors( scanline+1, len, vfile );
	}


	// read each component
	for( i = 0; i < 4; i++ ) {
	    for( int j = 0; j < len; ) {

			if( vfile->EndOfFile() ) {
				return false;
			}

			uint8 code = vfile->Read8();

			if( code > 128 ) {
			    code &= 127;
			    uint8 val = vfile->Read8();

				while( code-- ) {
					scanline[j++][i] = val;
				}
			}
			else  {
				while( code-- ) {
					if( vfile->EndOfFile() ) {
						return false;
					}

					scanline[j++][i] = vfile->Read8();
				}
			}
		}
    }

	if( vfile->OutOfBounds() ) {
		return false;
	}

	return true;
}


/** Read colors using the old scanline mode. */
bool PiHDRMap::OldReadColors( PiRGBE * scanline, int len, PiVirtualFile * vfile ) {

	int rshift = 0;

	while( len > 0 ) {

		if( vfile->EndOfFile() ) {
			return false;
		}

		scanline->r = vfile->Read8();
		scanline->g = vfile->Read8();
		scanline->b = vfile->Read8();
		scanline->e = vfile->Read8();

		if( vfile->OutOfBounds() ) {
			return false;
		}

		if( scanline->r == 1 && scanline->g == 1 && scanline->b == 1 ) {
			for( int i = int(scanline->e) << rshift; i > 0; i-- ) {
				scanline[0] = scanline[-1];
				scanline++;
				len--;
			}
			rshift += 8;
		}
		else {
			scanline++;
			len--;
			rshift = 0;
		}
	}

	return true;
}



/** Open an HDR light prove. */
bool PiLightProbe::Open( const char * fname, mode_t m ) {

	map.Close();

	const bool result = map.Open( fname );

	mode = m;
	if( mode == UNKNOWN ) {
		// try to guess it.
		const int w = map.GetWidth();
		const int h = map.GetHeight();

		if( w == h ) {
			mode = PROBE;
		}
		else if( 4 * w == 3 * h ) {
			mode = CROSS;
		}
	}

	return result;
}


/** Open an HDR light prove. */
void PiLightProbe::Close() {
	map.Close();
}


/** Sample the light prove. */
void PiLightProbe::Sample( FColor * out, const Vec3 & dir ) const {

	const int w = map.GetWidth();
	const int h = map.GetHeight();
	int x = 0;
	int y = 0;


	// Compute coordinates.
	if( mode == CROSS ) {
		// Get face size.
		int size = h / 4;
		piCheck( size == w / 3 );

		// Swap axis.
		Vec3 d( dir.x, dir.z, dir.y);

		// Determine face and project coordinates.
		int f;
		float fx, fy;
		if( fabs(d.x) > fabs(d.y) && fabs(d.x) > fabs(d.z) ) {
			if( d.x > 0 ) {
				f = 0;
				fx = -d.z / d.x;
				fy = -d.y / d.x;
			}
			else {
				f = 1;
				fx = d.z / -d.x;
				fy = -d.y / -d.x;
			}
		}
		else if( /*fabs(d.y) > fabs(d.x) &&*/ fabs(d.y) > fabs(d.z) ) {
			if( d.y > 0 ) {
				f = 2;
				fx = d.x / d.y;
				fy = d.z / d.y;
			}
			else {
				f = 3;
				fx = d.x / -d.y;
				fy = -d.z / -d.y;
			}
		}
		else /*if( fabs(d.z) > fabs(d.x) && fabs(d.z) > fabs(d.y) )*/ {
			if( d.z > 0 ) {
				f = 4;
				fx = d.x / d.z;
				fy = -d.y / d.z;
			}
			else {
				f = 5;
				fx = d.x / -d.z;	// in the vertical cross section this face is double-mirrored.
				fy = d.y / -d.z;
			}
		}

		// Offset coordinates.
		x = int( (fx * 0.5 + 0.5) * size ) + offset_array[f][0] * size;
		y = int( (fy * 0.5 + 0.5) * size ) + offset_array[f][1] * size;

		// Sample map.
		x = clamp(x, 0, w-1);
		y = clamp(y, 0, h-1);
		map.Sample( out, x, y );
	}
	else if( mode == PROBE ) {

		/*float m = 2 * sqrt(SQ(dir.z) + SQ(dir.y) + SQ(1 - dir.x));
		x = int((w-1) * (-dir.y / m + 0.5f));
		y = int((h-1) * (-dir.z / m + 0.5f));*/

		// The probes do not use the Lambert equal-area projection!
		float m = (1 / TWO_PI) * acos(dir.y) / sqrt(SQ(dir.x) + SQ(dir.z));
		x = int((w-1) * (dir.x * m + 0.5f));
		y = int((h-1) * (-dir.z * m + 0.5f));

		x = clamp(x, 0, w-1);
		y = clamp(y, 0, h-1);

		map.Sample( out, x, y );
	}
	else {
		*out = FColor::Black;
	}

}


/** Filter the map and output a spherical harmonic. */
bool PiLightProbe::PreFilter( SphericHarm * out, float exposure ) {

	const int w = map.GetWidth();
	const int h = map.GetHeight();

	if( mode == PROBE ) {
		piCheck( w == h );

		// This routine is from Ramamoorthi and Hanrahan.

		// The main integration routine.  Of course, there are better ways
		// to do quadrature but this suffices.  Calls updatecoeffs to
		// actually increment the integral. Width is the size of the
		// environment map

		int i, j;
		for( i = 0; i < w; i++ ) {
			for( j = 0; j < h; j++ ) {

				// We now find the cartesian components for the point (i,j)
				float u, v, r, theta, phi, domega;

				u = (j - w/2.0) / (w/2.0);	// u ranges from -1 to 1.
				v = (h/2.0 - i) / (h/2.0);  // v ranges from -1 to 1.

				// The radius.
				r = sqrt( SQ(u) + SQ(v) );

				if( r > 1.0 ) {
					// Consider only circle with r < 1.
					continue ;
				}

				// theta parameter of (i,j).
				theta = PI * r;

				// phi parameter.
				phi = atan2(v, u) + PI;

				// Cartesian components.
				Vec3 dir;
				dir.x = sin(theta) * sin(phi);
				dir.y = cos(theta);
				dir.z = sin(theta) * cos(phi);

				// Sample color.
				FColor color;
				map.Sample( &color, i, j );
				color *= exposure;
				color.Normalize();

				// Computation of the solid angle.  This follows from some
				// elementary calculus converting sin(theta) d theta d phi into
				// coordinates in terms of r.  This calculation should be redone
				// if the form of the input changes.
				domega = (2 * PI / w) * (2 * PI / h) * sinc(theta);

				// Update Integration.
				out->AddLightWeight( dir, color, domega );
			}
		}

		return true;
	}
	else if( mode == CROSS ) {

		int size = h / 4;
		piCheck( size == w/3 );

		for( int f = 0; f < 6; f++ ) {

			// compute offsets
			int offset_x = size * offset_array[f][0];
			int offset_y = size * offset_array[f][1];
			int x, y;

			for( y = 0; y < size; y++ ) {
				for( x = 0; x < size; x++ ) {

					float fx = 2 * (float(x)+0.5) / size - 1;
					float fy = 2 * (float(y)+0.5) / size - 1;

					if( f == 5 ) {
						fx = -fx;	// negative z is mirrored.
						fy = -fy;
					}

					// Get Cube dir.
					Vec3 dir;
					GetCubeDir( &dir, f, fx, fy );

					// Swap axis.
					Vec3 d(dir.x, dir.z, dir.y);


					// Sample color.
					FColor color;
					map.Sample( &color, x+offset_x, y+offset_y );
					color *= exposure;
					color.Normalize();

					// Compute the corresponding area element.
					float x0 = fx - (1.0f / size);
					float y0 = fy - (1.0f / size);
					float x1 = fx + (1.0f / size);
					float y1 = fy + (1.0f / size);
					float A = AreaElement(x0, y0) - AreaElement(x0, y1) - AreaElement(x1, y0) + AreaElement(x1, y1);

					// Update Integration.
					out->AddLightWeight( d, color, A );
				}
			}
		}
		return true;
	}

	return false;
}