//BIG NOTE on SUBTRACING:
//
//Subtracing is when we do everything all round.  This is achieved
//by effectively sliding all geometry up by .5 on all axes.  This
//effectively encapsulates all geometry in at least 1/2 of a block
//of virtual geometry.
//
//We then trace into this virtual geometry using the normal voxel
//method, and then once we might have a hit drop out.
//
//That's when subtrace takes over.  It starts by setting the corners
//of the current virtual voxel to be the blocks that the corners
//are in the center of (Remember, everything is offset by .5 in
//each direction).  Once This is done, it can step through the volume
//in small steps and compute the density (in accordance with mixfn.)
//Once it thinks it may have found the surface, it drops out and
//does a binary search refinement to find the actual surface.
//If no surface is found, it will drop back to the largescale ray
//tracer. If a surface is found, it will continue executing normally.




//This file assumes 8-character tabs.

//Direction of ray for this pixel (shooting off from Inital Camera)
//This is not necessairly normalized.
varying vec3 RayDirection;

//Normalized ray direction
vec3 dir;

//Actual [REAL] value for the position of the camera.  This is
//before any of the tricks described below are applied.
varying vec3 InitialCamera;


//So, because of floating point error, we must always keep the
//camera close to the origin.  We can fix this by virtually
//moving the camera around. This means we have the camera between
//0..1 and we shift the entire map by whole voxels to give the
//illusion the camera is moving.

//Actual camera position [0..1]
vec3 CurrentCamera;

//Amount to slide texture by [0..1]
vec3 CameraOffset;

//Virtual offset of camera (in whole units)
vec3 FloorCameraPos;


//Geometry Texture.  This contains the following information:
//Red Channel:   Block Type < not 
//Green Channel: Jumpmap XY
//Blue Channel:  Jumpmap Z ***THIS IS BROKEN***
//Alpha Channel: Density of Block - largescale trace only.
//			(is there a block here or not)
//This map, and AddTex contains one "pixel" for every voxel.
uniform sampler3D GeoTex;

//Additional Information Texture.  Like above, one pixel per voxel
//Red Channel:   Actual Cell to Draw
//Green Channel: Metadata [0..16/256]  (However, we may be increasing this)
//Blue Channel:  ?????? (unused) (probably will be lighting information)
//Alpha Channel: "Density" of block according to addtex (use this to render)
uniform sampler3D AddTex;

//Cell Attribute Map
//This contains detailed information about all cells.  See it's use for
//more information.  TileAttributes.txt contains what is put into this
//texture.
uniform sampler2D AttribMap;



//Use this if on ATI (it's a bug)
//ATI Cards for some reasons address all of the textures by a miniscule amount
//sideways... This corrects that bug.
//#define ATI
#ifdef ATI
const vec2 offset=-vec2( 1./8192.);
const vec3 lshw = -vec3( 1./65536.);
#else
const vec2 offset=vec2( 0. );
const vec3 lshw = vec3( 0. );
#endif

//3D position of expected location along ray.
vec3 ptr;

//Last read voxel from GeoTex.  It contains 'This' cell's information.
vec4 lastvox;

//Cellhit (for offset to cell with actual block data in it)
//This is only really used in subtrace mode.
// [0..1],[0..1],[0..1] for where in the cell contact was made.
vec3 cellhit = vec3( 0.5 ); 

//Similar - used in subtrace, but contains the offset for the selected point.
vec3 CellPoint = vec3( 0. );


//Uncomment this to override subtrace (may speed it up on some systems)
//By enabling the override, the GLSL compiler can drop some code on the floor
//thus decreasing the shader's footprint.
//In practice, this doesn't affect too much.
//#define SUBTRACE_OVERRIDE
#ifdef SUBTRACE_OVERRIDE
const float do_subtrace = false;
#else
uniform float do_subtrace;
#endif

//Size of voxel texture in pixels.
uniform float msX;
uniform float msY;
uniform float msZ;

//Multiplier to convert from world space coordinates into voxel map coordinates
vec3 msize = vec3( 1./msX, 1./msY, 1./msZ );

//total elapsed time.  Don't be shocked if this resets to zero. Considering
//this because over time, it could accumulate floating point error.
uniform float time;

//Normal of the surface we're hitting
vec3 normal = vec3( 0., 0., 1. );

//Total number of steps we've taken so far.
int step;

//Maximum steps we're permitted
const int maxsteps = 120;

//Maximum distance from camera we can go before it's treated as a "forever"
const float maxdist = 256.;



//For subtracing (these are where the local block pairs information goes).
//think of the values in these variables as the four points of a cube.
vec4 vecbot = vec4( 0., 0., 0., 0. );
vec4 vectop = vec4( 0., 0., 0., 0. );

//Local Distances - used for large-scale voxel tracing.
//This contains the farthest distance we can go on a given
//axis before we must consider checking to see if we hit something
vec3 dists;



//Subtrace tunable parameters:
const float mixval = .9;
const float densitylimit = .2;
const float linearstep = .1;
const int binaryRefinements = 8;
//This is the shape mixing function for subtrace.
vec3 mixfn( vec3 ins, float mixa )
{
//Wonky function...
//	return mix( cos( ins*3.14159 + 3.14159 )/2.0+0.5, stepss, mixa );

//Sharp edges...
//	return ins;

//Nice, smooth function we decided to use.
	vec3 coss = cos( ins*3.14159 + 3.14159 );
	vec3 sins = sign( coss );
	coss = abs( coss );
	coss = pow( coss, vec3( mixa ) );
	coss *= sins;

	return coss / 2.0 + 0.5;
}

//Large-scale traversal function
void TraverseIn( )
{
	//come up with vector to neutralize the sign on all computations
	//for traversal.  This makes it possible to always treat it like
	//we're tracing where all three direction components are positive.
	vec3 dircomps = -sign( dir ); //+1 if negative, -1 if positive

	//Floor behaves: -0.5 -> -1 / 0.5 -> 0
	//Frac behaves:  -0.1 -> .9, -0.5 -> .5 / -1 becomes 0

	//Load the firsxt voxel in.
	lastvox = texture3D( GeoTex, ( floor(ptr) )*msize  + CameraOffset );

	for( ; step < maxsteps && length( ptr )<maxdist; step++ )
	{
//		lastvox.g = 1./255.; //Disable X/Y jumpmap
		//Since the XY jumpmap works, leave this line commented.

		//Do not uncomment below line, because Z jumpmaps are in fact
		//broken!  If you uncomment this you will see visual artefacts.
		lastvox.b = 1./255.; //Disable Z jumpmap

		//Find the distance to the edges of our local cube.  These
		//are always positive values from 0 to 1.
		vec3 nextsteps = fract( ptr * dircomps  ) ;

		//Tack on the jumpmaps.  If we know it's clear in the XY direction
		//for 30 units, then enable the jumping in the XY direction by
		//that many units.
		nextsteps += vec3( lastvox.g, lastvox.g, lastvox.b ) *255. - 1.;

		//Find out how many units the intersection point between us and
		//the next intersection is in ray space.
		dists = nextsteps / abs(dir);

		//Find the closest axis.  We do this so we don't overshoot a hit.
		float mindist = 0.;
		mindist = min(dists.x,dists.y);
		mindist = min(mindist,dists.z);

		//Go there, plus a /tiny/ amount to prevent ourselves from hitting
		//an infinite loop.
		ptr += (mindist+0.0001) * dir;

		//Load the new voxel
		lastvox = texture3D( GeoTex, floor( ptr )*msize  + CameraOffset );

		//If it's a hit, we're good!  Return 
		if( lastvox.a > .1 )
		{
			return;
		}
	}

	//We ran out of runs, so we must pass a sentinal value into ptr.
	//-5000 means no geometry intersection.
	ptr = vec3( -5000. );
//	gl_FragColor = vec4( 0., 0., 1.,1. );
}

//You can make tileattributes update densities immediately by using the alternate (commented) code here.

float Density( vec3 ltexptr )
{
	return texture3D( AddTex, ltexptr*msize + CameraOffset ).a * 2.;
//	vec2 ltm = texture3D( AddTex, ltexptr*msize +CameraOffset ).rg;
//	return texture2D( AttribMap, vec2( ltm.g + (7./128.), ltm.r ) ).r;
}


void main()
{
	vec3 npos;

	//Configure all variables as described in the beginning.
	CurrentCamera = mod( InitialCamera, 1.0 );
	FloorCameraPos = floor( InitialCamera );
	CameraOffset = mod( (FloorCameraPos * msize), 1.0 ) + lshw;
	ptr = CurrentCamera + vec3(do_subtrace*0.5);
	dir = normalize(RayDirection);
	step = 0;

	//For subtracing
	vec3 nlc;
	vec3 dists;
	float minq;
	float minm=10000.;
	float mixtot=0.;
	vec3 lc = vec3( 0. );
	vec3 texptr;
	bool bfound;
	float fmarch;
	npos = vec3( 0. );

	do
	{
		//Find the next cell
		TraverseIn( );

		//We ran off the end of the map.
		if( ptr.x < -4000. )
		{
			//May want to do something funny here... Right now we
			//just set the sky to blue.
			gl_FragData[0] = vec4( 0., 0., 0., step );
			gl_FragData[1] = vec4( 0. );
			return;
		}

		//If we're not subtracing, we need to just calculate the normal
		//and get out!  The rest of this function is used for subtracing.
		if( do_subtrace < .5 )
		{
			cellhit = fract( ptr.xyz );

			vec3 lch = abs( cellhit-0.50 );
			if( lch.x > lch.y && lch.x > lch.z )
			{
				normal = vec3( sign(- dir.x ), 0., 0. );
			}
			else if ( lch.y > lch.z )
			{
				normal = vec3( 0., sign(- dir.y ), 0. );
			}
			else
			{
				normal = vec3( 0., 0., sign(- dir.z ) );
			}

			break;
		}

		//Calculate parameters for /this/ cell
		lc = fract(ptr.xyz);
		texptr = floor(ptr);

		//Handle inside-cell interpolation -> Find the values for the four surrounding cells.
/*		float ptA = texture3D( AddTex, ( texptr - vec3(0.0,0.0,0.0) + 0.0 )*msize +CameraOffset ).r>0.?1.:0.;
		float ptB = texture3D( AddTex, ( texptr - vec3(1.0,0.0,0.0) + 0.0 )*msize +CameraOffset ).r>0.?1.:0.;
		float ptC = texture3D( AddTex, ( texptr - vec3(0.0,1.0,0.0) + 0.0 )*msize +CameraOffset ).r>0.?1.:0.;
		float ptD = texture3D( AddTex, ( texptr - vec3(1.0,1.0,0.0) + 0.0 )*msize +CameraOffset ).r>0.?1.:0.;
		float ptE = texture3D( AddTex, ( texptr - vec3(0.0,0.0,1.0) + 0.0 )*msize +CameraOffset ).r>0.?1.:0.;
		float ptF = texture3D( AddTex, ( texptr - vec3(1.0,0.0,1.0) + 0.0 )*msize +CameraOffset ).r>0.?1.:0.;
		float ptG = texture3D( AddTex, ( texptr - vec3(0.0,1.0,1.0) + 0.0 )*msize +CameraOffset ).r>0.?1.:0.;
		float ptH = texture3D( AddTex, ( texptr - vec3(1.0,1.0,1.0) + 0.0 )*msize +CameraOffset ).r>0.?1.:0.;
		vectop = min(vec4( ptA, ptB, ptC, ptD )*1000.,1.);
		vecbot = min(vec4( ptE, ptF, ptG, ptH )*1000.,1.);
*/

		float ptA = Density( texptr - vec3(0.0,0.0,0.0) ); 
		float ptB = Density( texptr - vec3(1.0,0.0,0.0) ); 
		float ptC = Density( texptr - vec3(0.0,1.0,0.0) ); 
		float ptD = Density( texptr - vec3(1.0,1.0,0.0) ); 
		float ptE = Density( texptr - vec3(0.0,0.0,1.0) ); 
		float ptF = Density( texptr - vec3(1.0,0.0,1.0) ); 
		float ptG = Density( texptr - vec3(0.0,1.0,1.0) ); 
		float ptH = Density( texptr - vec3(1.0,1.0,1.0) ); 

		vectop = vec4( ptA, ptB, ptC, ptD );
		vecbot = vec4( ptE, ptF, ptG, ptH );

		npos = vec3( 0. );
		bfound = false;
		//Find the distance from one side of the cube to the other...
		nlc = (sign(dir)+1.0)/2.0 - sign(dir)* lc;
		dists = nlc / abs(dir);
		minq = -1000.;
		minm = 0.0;
		mixtot = 0.0;

		float mindist = 0.;
		if( dists.x <= dists.y && dists.x <= dists.z )
			mindist = dists.x;
		else if( dists.y <= dists.x && dists.y <= dists.z )
			mindist = dists.y;
		else
			mindist = dists.z;



		//Linearly search through the block, trying to find the intersection.
		for( fmarch = 0.0; fmarch <= mindist+linearstep*.5; fmarch+=linearstep )
		{
			npos = lc + dir * (fmarch);

			//You may notice here - we don't actually shoot a ray through.
			//We warp the ray as a function of mixfn.
			vec3 tnpos = mixfn( npos, mixval );
			mixtot = mix( 	mix( mix(vecbot.a, vecbot.b, tnpos.x ),
					mix( vecbot.g, vecbot.r, tnpos.x ), tnpos.y ),
					mix( mix(vectop.a, vectop.b, tnpos.x ),
					mix( vectop.g, vectop.r, tnpos.x ), tnpos.y ),
					tnpos.z );
			if( mixtot > densitylimit ) { bfound = true; minm = fmarch; break; }
			if( mixtot > minq ) { minq = mixtot; minm = fmarch; }
		}

		if( bfound )
			break;

		//If not found, keep going on...
		step+=2;
	} while( step < maxsteps );
	if( step >= maxsteps )
	{
		gl_FragData[0] = vec4( 0., 0., 0., step );
		gl_FragData[1] = vec4( 0., 0., 0., -1. );
		return;
	}

	if( do_subtrace > 0.5 )
	{
		//Binary search the remaining space.
		//This means we're a hit, and just want to get our XYZ location of the hit very precice.
		minm -= linearstep*0.5;
		fmarch = minm;

		float fmultiplier = linearstep * 0.5;
		for( int i = 0; i < binaryRefinements; i++ )
		{
			npos = lc + dir * (fmarch);
			npos = 1.-npos;
			vec3 tnpos = mixfn( npos, mixval );
			mixtot = mix( 	mix( mix( vectop.r, vectop.g, tnpos.x ),
					mix(vectop.b, vectop.a, tnpos.x ), tnpos.y ),
					mix( mix( vecbot.r, vecbot.g, tnpos.x ),
					mix(vecbot.b, vecbot.a, tnpos.x ), tnpos.y ),
					tnpos.z );

			float fmuxsign = .5;
			if( mixtot > densitylimit ) { bfound = true; fmuxsign = -.5; }
			fmultiplier=abs(fmultiplier) * fmuxsign;
			fmarch += fmultiplier;
		}

		//Find normal...
		//We do this by inching a very small amount in each direction
		//to compute the gradiant of the density of the metasurface.
		//The normal is actually the normalized gradient.
		vec3 npx = npos + vec3( 0.01, 0., 0. );
		npx = mixfn( npx, mixval );

		float mixtotx = mix( 	mix( mix( vectop.r, vectop.g, npx.x ), mix(vectop.b, vectop.a, npx.x ), npx.y ),
					mix( mix( vecbot.r, vecbot.g, npx.x ), mix(vecbot.b, vecbot.a, npx.x ), npx.y ), npx.z );

		vec3 npy = npos + vec3( 0.0, 0.01, 0. );
		npy = mixfn( npy, mixval );

		float mixtoty = mix( 	mix( mix( vectop.r, vectop.g, npy.x ), mix(vectop.b, vectop.a, npy.x ), npy.y ),
					mix( mix( vecbot.r, vecbot.g, npy.x ), mix(vecbot.b, vecbot.a, npy.x ), npy.y ), npy.z );

		vec3 npz = npos + vec3( 0., 0., 0.01 );
		npz = mixfn( npz, mixval );

		float mixtotz = mix( 	mix( mix( vectop.r, vectop.g, npz.x ), mix(vectop.b, vectop.a, npz.x ), npz.y ),
					mix( mix( vecbot.r, vecbot.g, npz.x ), mix(vecbot.b, vecbot.a, npz.x ), npz.y ), npz.z );
		normal = normalize( vec3( mixtotx, mixtoty, mixtotz ) - mixtot );

		ptr += dir * fmarch;
	}

	gl_FragData[0] = vec4( ptr, step );
	gl_FragData[1] = vec4( normal.xyz, 1. );
	return;
/*
	if( do_subtrace >= 0.5 )
	{
//#define DISABLE_FUZZYEDGES
#ifndef DISABLE_FUZZYEDGES
		vec3 texptr = floor( ptr.xyz );
		float ptA = texture3D( GeoTex, ( texptr - vec3(0.0,0.0,0.0) + 0.0 )*msize + CameraOffset ).r>0.?1.:0.;
		float ptB = texture3D( GeoTex, ( texptr - vec3(1.0,0.0,0.0) + 0.0 )*msize + CameraOffset ).r>0.?1.:0.;
		float ptC = texture3D( GeoTex, ( texptr - vec3(0.0,1.0,0.0) + 0.0 )*msize + CameraOffset ).r>0.?1.:0.;
		float ptD = texture3D( GeoTex, ( texptr - vec3(1.0,1.0,0.0) + 0.0 )*msize + CameraOffset ).r>0.?1.:0.;
		float ptE = texture3D( GeoTex, ( texptr - vec3(0.0,0.0,1.0) + 0.0 )*msize + CameraOffset ).r>0.?1.:0.;
		float ptF = texture3D( GeoTex, ( texptr - vec3(1.0,0.0,1.0) + 0.0 )*msize + CameraOffset ).r>0.?1.:0.;
		float ptG = texture3D( GeoTex, ( texptr - vec3(0.0,1.0,1.0) + 0.0 )*msize + CameraOffset ).r>0.?1.:0.;
		float ptH = texture3D( GeoTex, ( texptr - vec3(1.0,1.0,1.0) + 0.0 )*msize + CameraOffset ).r>0.?1.:0.;

		vec4 PN1 = vec4( pNoise( ptr.xyz + 1.0 ), pNoise( ptr.xyz + 2.0 ), pNoise( ptr.xyz + 3.0 ), pNoise( ptr.xyz + 4.0 ) );
		vec4 PN2 = (1.-PN1);

		vectop = vec4( ptA, ptB, ptC, ptD ) * (PN1*.2+1.);
		vecbot = vec4( ptE, ptF, ptG, ptH ) * (PN2*.2+1.);
#endif
		vec4 vectopsA = vectop * vec4( 1.-ppmod.x,ppmod.x,1.-ppmod.x,ppmod.x ) * (1.-ppmod.z) * vec4( vec2( 1.-ppmod.y ), vec2(ppmod.y ) );
		vec4 vecbotsA = vecbot * vec4( 1.-ppmod.x,ppmod.x,1.-ppmod.x,ppmod.x ) * (ppmod.z) * vec4( vec2( 1.-ppmod.y ), vec2(ppmod.y ) );
		vec4 minsA = max( vectopsA, vecbotsA );
		vec2 minsB = max( minsA.xy, minsA.zw );
		float amin = max( minsB.x, minsB.y );
		if( vectopsA.x == amin ) CellPoint = vec3( 0., 0., 0. );
		else if( vectopsA.y == amin ) CellPoint = vec3( 1., 0., 0. );
		else if( vectopsA.z == amin ) CellPoint = vec3( 0., 1., 0. );
		else if( vectopsA.w == amin ) CellPoint = vec3( 1., 1., 0. );
		else if( vecbotsA.x == amin ) CellPoint = vec3( 0., 0., 1. );
		else if( vecbotsA.y == amin ) CellPoint = vec3( 1., 0., 1. );
		else if( vecbotsA.z == amin ) CellPoint = vec3( 0., 1., 1. );
		else if( vecbotsA.w == amin ) CellPoint = vec3( 1., 1., 1. );
		else
		{
//			gl_FragColor = vec4( 1., 0., 1., 0. );
			return;
		}
	}
	else
	{
		CellPoint = vec3(0.0);
	}
//	gl_FragColor = vec4( ppmod, 1.0 );
	return;

	vec3 GlobalPos = FloorCameraPos + ptr;
//	gl_FragColor = vec4( CellPoint, 1. );
//	return;
	ppmod = mod( ppmod -(do_subtrace*.5), 1.0 );

	vec4 ExtraData = texture3D( AddTex, floor( ptr - CellPoint)*msize + CameraOffset );
	ExtraData.g = 0.;
//gl_FragColor = vec4( mod( ExtraData.r * 16., 1. ), ExtraData.r, ExtraData.g*16., 0. );
//return;
	vec4 BaseColor = texture2D( AttribMap, vec2( ExtraData.g + (0./128.), ExtraData.r ) );
	vec4 NoiseColor = texture2D( AttribMap, vec2( ExtraData.g + (1./128.), ExtraData.r ) );
	vec4 NoiseSet = texture2D( AttribMap, vec2( ExtraData.g + (2./128.), ExtraData.r ) );
	vec4 NoiseMux = texture2D( AttribMap, vec2( ExtraData.g + (3./128.), ExtraData.r ) );

	vec4 CoreData = texture2D( AttribMap, vec2( ExtraData.g + (4./128.), ExtraData.r ) ); //For trees, etc.
	vec4 TimeSettings = texture2D( AttribMap, vec2( ExtraData.g + (5./128.), ExtraData.r ) );

	vec4 Speckles = texture2D( AttribMap, vec2( ExtraData.g + (6./128.), ExtraData.r ) );
//	vec4 ExtraSet = texture2D( AttribMap, vec2( ExtraData.g + (7./128.), ExtraData.r ) );

//	gl_FragColor = vec4( BaseColor );

	vec3 noiseplace = NoiseSet.xyz * vec3( GlobalPos ) + TimeSettings.xyz * time;
	float noise = pNoise( noiseplace ) * NoiseMux.r + pNoise( noiseplace * 2. ) * NoiseMux.g + pNoise( noiseplace * 4. ) * NoiseMux.b + pNoise( noiseplace * 8. ) * NoiseMux.a + NoiseColor.a;

	float core = length( (ppmod.xyz-.5) * vec3( 1., 1., CoreData.y ) ) * CoreData.x;
	core = mod( core+noise, CoreData.z ) * CoreData.a;

	noise = noise+core;

	//if( noise > NoiseSet.a-10.0 )
	if( true )
	{

		vec3 NoiseOut = noise * ((noise<Speckles.a)?Speckles.rgb:NoiseColor.rgb);
		noise = abs( noise );

		vec3 OutColor = BaseColor.rgb + NoiseOut;

		float normalcomponent = max( BaseColor.a, dot(normal, vec3( .58, .58, .58 )));

		OutColor *= normalcomponent;
//		gl_FragColor = vec4( OutColor, 1. );
		return;


	}

return;
//In case we have an issue dropping out...
//} while( step < maxsteps );


	return;
*/
}