struct T_PBRMaterial
{
	vec3 cAlbedo, cSpecular;
	float roughness , anisotropy , subsurface , metallic;
};

struct T_PBRPrecomputed
{
	float nDotC;
	float ffndc;			// Fresnel from n.c
	vec3 tangent, bitangent;
	float specAlpha , viewSpecular;
	float aAspectX , aAspectY;	// Anisotropy
};


float PBR_FresnelFrom( in float dotProduct ) {
	float d = clamp( 1.0 - dotProduct , 0.0 , 1.0 );
	return d * d * d * d * d;
}

float PBR_GetSpecular( in float cosine , in float alpha )
{
	float cs = cosine * cosine;
	float as = alpha * alpha;
	return 1. / ( cosine + sqrt( cs + as - cs * as ) );
}

// Precompute some of the material's properties. This is independant of the
// light source.
T_PBRPrecomputed PBR_Precompute(
		in T_PBRMaterial material ,
		in vec3 rayDir ,
		in vec3 normal )
{
	T_PBRPrecomputed rv;

	rv.nDotC = dot( normal , rayDir );
	rv.ffndc = PBR_FresnelFrom( rv.nDotC );

	rv.tangent = cross( vec3( 0. , 1. , 0. ) , normal );
	if ( length( rv.tangent ) == 0.0 ) {
		rv.tangent = cross( vec3( 1. , 0. , 0. ) , normal );
	}
	rv.tangent = normalize( rv.tangent );
	rv.bitangent = normalize( cross( normal , rv.tangent ) );

	rv.specAlpha = pow( material.roughness * .5 + .5 , 2. );
	rv.viewSpecular = PBR_GetSpecular( rv.nDotC , rv.specAlpha );

	const float sRoughness = material.roughness * material.roughness;
	const float aspect = sqrt( 1.0 - material.anisotropy * .9 );
	rv.aAspectX = max( .001, sRoughness / aspect );
	rv.aAspectY = max( .001, sRoughness * aspect );

	return rv;
}

// Actually compute a light source's contribution
vec3 PBR_Shade(
		in T_PBRMaterial material ,
		in T_PBRPrecomputed precomputed ,
		in vec3 rayDir ,
		in vec3 normal ,
		in vec3 lightDir )
{
	const float nDotL = dot( normal , lightDir );
	if ( nDotL <= 0. ) {
		return vec3( 0. );
	}

	const vec3 halfVec = normalize( lightDir + rayDir );
	const float nDotH = dot( normal , halfVec );
	const float lDotH = dot( lightDir , halfVec );

	const float ffndl = PBR_FresnelFrom( nDotL );
	float grazingDiffuse = lDotH * lDotH * material.roughness;
	float dSubsurface = mix( 1.0 , grazingDiffuse , ffndl )
		* mix( 1.0 , grazingDiffuse , precomputed.ffndc );
	dSubsurface = 1.25 * ( dSubsurface * ( 1.0 / ( nDotL + precomputed.nDotC ) - .5 ) + .5 );
	grazingDiffuse = .5 + 2. * grazingDiffuse;
	const float dFresnel = mix( 1.0 , grazingDiffuse , ffndl )
		* mix( 1.0 , grazingDiffuse , precomputed.ffndc );

	float specular = PBR_GetSpecular( nDotL , precomputed.specAlpha )
		* precomputed.viewSpecular;
	specular = mix( specular , 1.0 , PBR_FresnelFrom( lDotH ) );

	const vec3 d = vec3(
		dot( halfVec , precomputed.tangent ) / precomputed.aAspectX ,
		dot( halfVec , precomputed.bitangent ) / precomputed.aAspectY ,
		nDotH );
	const float ds = dot( d , d );
	const float anisotropic = precomputed.aAspectX * precomputed.aAspectY
		* ds * ds * 3.14159265;

	return nDotL * ( material.cAlbedo
			* mix( dFresnel , dSubsurface , material.subsurface )
			* pow( 1.0 - material.metallic , 3 )
		+ specular * material.cSpecular / anisotropic
	);
}