#ifndef _BRDF_GLSL_
#define _BRDF_GLSL_

// http://xlgames-inc.github.io/posts/improvedibl/
// http://blog.selfshadow.com/publications/s2013-shading-course/
vec3 f_schlick(const vec3 f0, const float vh) {
	return f0 + (1.0 - f0) * exp2((-5.55473 * vh - 6.98316) * vh);
}

float v_smithschlick(const float nl, const float nv, const float a) {
	return 1.0 / ((nl * (1.0 - a) + a) * (nv * (1.0 - a) + a));
}

//Uncorrelated masking/shadowing (info below) function
//Because it is uncorrelated, G1(NdotL, a) gives us shadowing, and G1(NdotV, a) gives us masking function.
//Approximation from: https://ubm-twvideo01.s3.amazonaws.com/o1/vault/gdc2017/Presentations/Hammon_Earl_PBR_Diffuse_Lighting.pdf
float g1_approx(const float NdotX, const float alpha)
{
	return (2.0 * NdotX) * (1.0 / (NdotX * (2.0 - alpha) + alpha));
}

//Uncorrelated masking-shadowing function
//Approximation from: https://ubm-twvideo01.s3.amazonaws.com/o1/vault/gdc2017/Presentations/Hammon_Earl_PBR_Diffuse_Lighting.pdf
float g2_approx(const float NdotL, const float NdotV, const float alpha)
{
	vec2 helper = (2.0 * vec2(NdotL, NdotV)) * (1.0 / (vec2(NdotL, NdotV) * (2.0 - alpha) + alpha));
	return max(helper.x * helper.y, 0.0); //This can go negative, let's fix that
}

float d_ggx(const float nh, const float a) {
	float a2 = a * a;
	float denom = nh * nh * (a2 - 1.0) + 1.0;
	denom = max(denom * denom, 0.00006103515625 /* 2^-14 = smallest possible half float value, prevent div by zero */);
	return a2 * (1.0 / 3.1415926535) / denom;
}

vec3 specularBRDF(const vec3 f0, const float roughness, const float nl, const float nh, const float nv, const float vh) {
	float a = roughness * roughness;
	vec3 result = d_ggx(nh, a) * g2_approx(nl, nv, a) * f_schlick(f0, vh) / max(4.0 * nv, 1e-5); //NdotL cancels out later
	return result;
}

// John Hable - Optimizing GGX Shaders
// http://filmicworlds.com/blog/optimizing-ggx-shaders-with-dotlh/
vec3 specularBRDFb(const vec3 f0, const float roughness, const float dotNL, const float dotNH, const float dotLH) {
	// D
	const float pi = 3.1415926535;
	float alpha = roughness * roughness;
	float alphaSqr = alpha * alpha;
	float denom = dotNH * dotNH * (alphaSqr - 1.0) + 1.0;
	float D = alphaSqr / (pi * denom * denom);
	// F
	const float F_a = 1.0;
	float F_b = pow(1.0 - dotLH, 5.0);
	// V
	float vis;
	float k = alpha / 2.0;
	float k2 = k * k;
	float invK2 = 1.0 - k2;
	vis = 1.0 / (dotLH * dotLH * invK2 + k2);
	vec2 FV_helper = vec2((F_a - F_b) * vis, F_b * vis);

	vec3 FV = f0 * FV_helper.x + FV_helper.y;
	vec3 specular = clamp(dotNL, 0.0, 1.0) * D * FV;
	return specular / 4.0; // TODO: get rid of / 4.0
}

vec3 orenNayarDiffuseBRDF(const vec3 albedo, const float roughness, const float nv, const float nl, const float vh) {
	float a = roughness * roughness;
	float s = a;
	float s2 = s * s;
	float vl = 2.0 * vh * vh - 1.0; // Double angle identity
	float Cosri = vl - nv * nl;
	float C1 = 1.0 - 0.5 * s2 / (s2 + 0.33);
	float test = 1.0;
	if (Cosri >= 0.0) test = (1.0 / (max(nl, nv)));
	float C2 = 0.45 * s2 / (s2 + 0.09) * Cosri * test;
	return albedo * max(0.0, nl) * (C1 + C2) * (1.0 + roughness * 0.5);
}

vec3 lambertDiffuseBRDF(const vec3 albedo, const float nl) {
	return albedo * (1.0 / 3.1415926535) * nl;
}

vec3 surfaceAlbedo(const vec3 baseColor, const float metalness) {
	return mix(baseColor, vec3(0.0), metalness);
}

vec3 surfaceF0(const vec3 baseColor, const float metalness) {
	return mix(vec3(0.04), baseColor, metalness);
}

float getMipFromRoughness(const float roughness, const float numMipmaps) {
	// First mipmap level = roughness 0, last = roughness = 1
	return roughness * numMipmaps;
}

float wardSpecular(vec3 N, vec3 H, float dotNL, float dotNV, float dotNH, vec3 fiberDirection, float shinyParallel, float shinyPerpendicular) {
	if(dotNL < 0.0 || dotNV < 0.0) {
		return 0.0;
	}
	// fiberDirection - parse from rotation
	// shinyParallel - roughness
	// shinyPerpendicular - anisotropy

	vec3 fiberParallel = normalize(fiberDirection);
	vec3 fiberPerpendicular = normalize(cross(N, fiberDirection));
	float dotXH = dot(fiberParallel, H);
	float dotYH = dot(fiberPerpendicular, H);
	const float PI = 3.1415926535;
	float coeff = sqrt(dotNL/dotNV) / (4.0 * PI * shinyParallel * shinyPerpendicular);
	float theta = (pow(dotXH/shinyParallel, 2.0) + pow(dotYH/shinyPerpendicular, 2.0)) / (1.0 + dotNH);
	return clamp(coeff * exp(-2.0 * theta), 0.0, 1.0);
}

// https://www.unrealengine.com/en-US/blog/physically-based-shading-on-mobile
// vec3 EnvBRDFApprox(vec3 SpecularColor, float Roughness, float NoV) {
// 	const vec4 c0 = { -1, -0.0275, -0.572, 0.022 };
// 	const vec4 c1 = { 1, 0.0425, 1.04, -0.04 };
// 	vec4 r = Roughness * c0 + c1;
// 	float a004 = min( r.x * r.x, exp2( -9.28 * NoV ) ) * r.x + r.y;
// 	vec2 AB = vec2( -1.04, 1.04 ) * a004 + r.zw;
// 	return SpecularColor * AB.x + AB.y;
// }
// float EnvBRDFApproxNonmetal(float Roughness, float NoV) {
// 	// Same as EnvBRDFApprox( 0.04, Roughness, NoV )
// 	const vec2 c0 = { -1, -0.0275 };
// 	const vec2 c1 = { 1, 0.0425 };
// 	vec2 r = Roughness * c0 + c1;
// 	return min( r.x * r.x, exp2( -9.28 * NoV ) ) * r.x + r.y;
// }
float D_Approx(const float Roughness, const float RoL) {
	float a = Roughness * Roughness;
	float a2 = a * a;
	float rcp_a2 = 1.0 / a2;//rcp(a2);
	// 0.5 / ln(2), 0.275 / ln(2)
	float c = 0.72134752 * rcp_a2 + 0.39674113;
	return rcp_a2 * exp2( c * RoL - c );
}

#ifdef _ClearCoat
vec3 clearcoatBRDF(const float clearcoat, const float clearcoat_rough,
	const float coat_ior,
	const float dotNL, const float dotNH, const float dotNV, const float dotVH) {
	if (clearcoat <= 0.0) return vec3(0.0);
	float a = clearcoat_rough * clearcoat_rough;
	// F0 from Fresnel equation for dielectric
	float ccF0 = (coat_ior - 1.0) / (coat_ior + 1.0);
	ccF0 = ccF0 * ccF0;
	float F = ccF0 + (1.0 - ccF0) * pow(1.0 - dotVH, 5.0);
	float D = d_ggx(dotNH, a);
	float G = g2_approx(dotNL, dotNV, a);
	return vec3(clearcoat * D * G * F / max(4.0 * dotNV, 1e-5));
}

float coatAttenuation(const float clearcoat,
	const float coat_ior, const float dotNV) {
	if (clearcoat <= 0.0) return 1.0;
	float ccF0 = (coat_ior - 1.0) / (coat_ior + 1.0);
	ccF0 = ccF0 * ccF0;
	// Schlick with pow(.,5) - cheaper than exp2
	float F = ccF0 + (1.0 - ccF0) * pow(1.0 - dotNV, 5.0);
	return max(1.0 - F * clearcoat, 0.0);
}

vec3 coatTintAttenuation(const float clearcoat, const vec3 coat_tint, const float dotNV) {
	if (clearcoat <= 0.0) return vec3(1.0);
	float absorption = 1.0 / max(dotNV, 0.3);
	return mix(vec3(1.0), coat_tint, clamp(absorption * 0.2, 0.0, 1.0));
}
#endif

#ifdef _Sheen
// based on Blender sheen model/Frostbite PBR
vec3 sheenBRDF(const float sheen, const float sheen_rough,
	const vec3 sheen_tint, const float dotNL, const float dotNH, const float dotNV) {
	if (sheen <= 0.0) return vec3(0.0);
	float rough = clamp(sheen_rough, 1e-3, 1.0);
	float a = rough * rough;
	float sinNH2 = 1.0 - dotNH * dotNH;
	float a2 = a * a;
	float D = (2.0 + a2) * sinNH2 / (2.0 * 3.1415926535 * pow(1.0 + a2 * sinNH2, 2.0));
	float V = 1.0 / (4.0 * dotNL * dotNV + 1e-5);
	float sheenAlbedo = (1.0 - 0.5 * rough) * mix(1.0, dotNV, 0.5);
	return sheen_tint * sheen * D * V * dotNL * sheenAlbedo;
}

float sheenAttenuation(const float sheen, const float sheen_rough,
	const vec3 sheen_tint, const float dotNV) {
	if (sheen <= 0.0) return 1.0;
	float rough = clamp(sheen_rough, 1e-3, 1.0);
	float sheenAlbedo = (1.0 - 0.5 * rough) * mix(1.0, dotNV, 0.5);
	float maxComp = sheen * max(max(sheen_tint.r, sheen_tint.g), sheen_tint.b) * sheenAlbedo;
	return max(1.0 - maxComp, 0.0);
}
#endif

#ifdef _Anisotropy
// anisotropic GGX Burley 2012
vec3 anisotropicBRDF(const vec3 f0, const float roughness, const float anisotropy,
	const float aniso_rot, const vec3 tangent, const vec3 bitangent,
	const vec3 n, const vec3 l, const vec3 v,
	const float dotNL, const float dotNV) {
	if (abs(anisotropy) <= 0.001) return vec3(0.0);
	float rot = aniso_rot * 3.1415926535 * 2.0;
	float cr = cos(rot);
	float sr = sin(rot);
	vec3 t = normalize(tangent * cr + bitangent * sr);
	vec3 b = normalize(bitangent * cr - tangent * sr);
	float aniso_abs = abs(anisotropy);
	float at = max(roughness * (1.0 + aniso_abs), 1e-5);
	float ab = max(roughness * (1.0 - aniso_abs), 1e-5);
	if (anisotropy < 0.0) { vec3 tmp = t; t = b; b = tmp; }
	float at2 = at * at;
	float ab2 = ab * ab;
	vec3 h = normalize(l + v);
	float dotTH = dot(t, h);
	float dotBH = dot(b, h);
	float dotTV = dot(t, v);
	float dotBV = dot(b, v);
	float dotTL = dot(t, l);
	float dotBL = dot(b, l);
	float denom = max(dotTH * dotTH / at2 + dotBH * dotBH / ab2, 1e-7);
	float D = 1.0 / (3.1415926535 * at * ab * denom * denom);
	float V = 1.0 / max(dotNL * (dotTL / at + dotBL / ab) * (dotTV / at + dotBV / ab), 1e-5);
	float dotVH = max(dot(v, h), 0.0);
	vec3 F = f_schlick(f0, dotVH);
	return D * V * F / max(4.0 * dotNV, 1e-5);
}
#endif

#ifdef _Subsurface
// Blenders bssrdf_burley implementation
vec3 subsurfaceBRDF(const vec3 albedo, const vec3 sss_color,
	const vec3 sss_radius, const float subsurface, const float sss_anisotropy,
	const float dotNL) {
	if (subsurface <= 0.0) return vec3(0.0);
	vec3 mfp = sss_radius * (0.25 / 3.1415926535);
	vec3 A = clamp(albedo, 0.0, 1.0);
	vec3 d = 1.9 - A + 3.5 * (A - 0.8) * (A - 0.8);
	d = mfp / max(d, 1e-5);
	float aniso = clamp(sss_anisotropy, 0.0, 0.9);
	float scatter = subsurface * (1.0 / 3.1415926535);
	vec3 sssDiffuse = sss_color * scatter * dotNL;
	float backScatter = max(0.0, 1.0 - dotNL) * (1.0 - aniso) * 0.5;
	vec3 sssBack = sss_color * subsurface * backScatter;
	float dist = max(0.0, 1.0 - dotNL);
	vec3 rcp_d = 1.0 / max(d, vec3(1e-5));
	vec3 x = vec3(dist) * rcp_d;
	vec3 extinction = 1.0 / (1.0 + x + 0.5 * x * x);
	return sssDiffuse * extinction + sssBack * (vec3(1.0) - extinction);
}
#endif

#ifdef _Transmission
// Blenders microfacet glass/refraction model
vec3 transmissionBRDF(const vec3 albedo, const float transmission,
	const float trans_rough, const float ior, const float thin_wall,
	const float dotNL, const float dotNV, const float dotVH) {
	if (transmission <= 0.0) return vec3(0.0);
	float F0 = (ior - 1.0) / (ior + 1.0);
	F0 = F0 * F0;
	float F = F0 + (1.0 - F0) * pow(1.0 - dotVH, 5.0);
	float transmittance = 1.0 - F;
	if (thin_wall > 0.5) {
		return albedo * transmission * transmittance * dotNL;
	}
	float a = trans_rough * trans_rough;
	float rough_atten = mix(1.0, 1.0 / max(dotNV, 0.1), a);
	return albedo * transmission * transmittance * rough_atten * dotNL;
}
#endif

#endif