Diffuse BRDF

gistfile1.txt

// BRDFs

#define PI 3.14159265359f
#define PI2 6.28318530718f
#define RECIPROCAL_PI 0.31830988618f
#define RECIPROCAL_PI2 0.15915494f
#define LOG2 1.442695f
#define EPSILON 1e-6
#define pow2(x) ((x)*(x))
#define pow3(x) pow2(x)*(x)
#define pow4(x) pow2(x)*pow2(x)
#define pow5(x) pow2(x)*pow2(x)*(x)

inline float saturate(float a)
{
	if (a > 1.0f) return 1.0f;
	if (a < 0.0f) return 0.0f;
	return a;
}

inline Vector3 saturate(const Vector3& v)
{
	return minPerElem(maxPerElem(v, Vector3(0.0f)), Vector3(1.0f));
}

// GLSL compatible
inline float step(float edge, float x)
{
	return (x < edge) ? 0.0f : 1.0f;
}

inline float mix(float x, float y, float a)
{
	return x*(1.0f - a) + y*a;
}

inline float smoothstep(float edge0, float edge1, float x)
{
	if (edge0 >= edge1) return 0.0f;
	float t = saturate((x - edge0) / (edge1 - edge0));
	return t * t * (3.0f - 2.0f * t);
}

float F_SchlickScalar(float f0, float f90, float u)
{
	return f0 + (f90 - f0) * pow(1.0f - u, 5.0f);
}

Vector3 LambertDiffuse(const Vector3& diffuseColor)
{
	return diffuseColor / PI;
}

Vector3 BurleyDiffuse(const Vector3& diffuseColor, float a, const Vector3& N, const Vector3& V, const Vector3& L)
{
	float dotNL = saturate(dot(N, L));
	float dotNV = saturate(dot(N, V));
	Vector3 H = normalize(L + V);
	float dotLH = saturate(dot(L, H));
	float fd90 = 0.5f + 2.0f * dotLH * dotLH * a;
	float nl = F_SchlickScalar(1.0f, fd90, dotNL);
	float nv = F_SchlickScalar(1.0f, fd90, dotNV);
	return diffuseColor * ((nl*nv) / PI);
}

// Frostbite
Vector3 RenormalizedBurleyDiffuse(const Vector3& diffuseColor, float roughness, const Vector3& N, const Vector3& V, const Vector3& L)
{
	float dotNL = saturate(dot(N, L));
	float dotNV = saturate(dot(N, V));
	Vector3 H = normalize(L + V);
	float dotLH = saturate(dot(L, H));
	float energyBias = mix(0.0f, 0.5f, roughness);
	float energyFactor = mix(1.0f, 1.0f / 1.51f, roughness);
	float fd90 = energyBias + 2.0f * dotLH * dotLH * roughness;
	float nl = F_SchlickScalar(1.0f, fd90, dotNL);
	float nv = F_SchlickScalar(1.0f, fd90, dotNV);
	return diffuseColor * ((nl*nv*energyFactor) / PI);
}

// http://ruh.li/GraphicsOrenNayar.html
Vector3 OrenNayarDiffuse(const Vector3& diffuseColor, float a, const Vector3& N, const Vector3& V, const Vector3& L)
{
	float dotNL = saturate(dot(N, L));
	float dotNV = saturate(dot(N, V));
	Vector3 H = normalize(L + V);
	float dotLH = saturate(dot(L, H));

	Vector3 v1 = V - N * dotNV;
	Vector3 v2 = L - N * dotNL;
	float theta_r = acosf(dotNV);
	float sigma2 = a;
	float cos_phi_diff = dot(normalize(V - N * dotNV), normalize(L - N * dotNL));
	if (isnan(cos_phi_diff)) cos_phi_diff = 1.0f;
	float theta_i = acosf(dotNL);
	float alpha = max(theta_i, theta_r);
	float beta = min(theta_i, theta_r);
	if (alpha > PI / 2) return Vector3(0.0f);

	float C1 = 1.0f - 0.5f * sigma2 / (sigma2 + 0.33f);
	float C2 = 0.45f * sigma2 / (sigma2 + 0.09f);
	if (cos_phi_diff >= 0.0f) C2 *= sinf(alpha);
	else C2 *= (sinf(alpha) - pow(2 * beta / PI, 3.0f));
	float C3 = 0.125f * sigma2 / (sigma2 + 0.09f) * pow((4 * alpha * beta) / (PI*PI), 2.0f);
	float L1 = (C1 + cos_phi_diff * C2 * tanf(beta) + (1.0f - abs(cos_phi_diff)) * C3 * tanf((alpha + beta) / 2));
	float L2 = 0.17f * (sigma2 / (sigma2 + 0.13f)) * (1.0f - cos_phi_diff * (4.0f * beta * beta) / (PI*PI));
	Vector3 rho_div_pi = diffuseColor / PI;
	Vector3 rho_rho_div_pi = mulPerElem(rho_div_pi, diffuseColor);
	return rho_div_pi * L1 + rho_rho_div_pi * L2;
}

// http://ruh.li/GraphicsOrenNayar.html
Vector3 QualitativeOrenNayarDiffuse(const Vector3& diffuseColor, float a, const Vector3& N, const Vector3& V, const Vector3& L)
{
	// calculate intermediary values
	float dotNL = saturate(dot(N, L));
	float dotNV = saturate(dot(N, V));
	float dotLV = saturate(dot(L, V));
	Vector3 H = normalize(L + V);
	float dotLH = saturate(dot(L, H));

	float angleVN = acosf(dotNV);
	float angleLN = acosf(dotNL);

	float alpha = max(angleVN, angleLN);
	float beta = min(angleVN, angleLN);
	float gamma = dot(normalize(V - N * dotNV), normalize(L - N * dotNL));
	if (isnan(gamma)) gamma = 1.0f;

	float roughnessSquared = a;

	// calculate A and B
	float A = 1.0f - 0.5f * (roughnessSquared / (roughnessSquared + 0.57f));
	float B = 0.45f * (roughnessSquared / (roughnessSquared + 0.09f));
	float C = sinf(alpha) * tanf(beta);

	// put it all together
	//float L1 = max(0.0f, dotNL) * (A + B * max(0.0f, gamma) * C);
	float L1 = (A + B * max(0.0f, gamma) * C);

	Vector3 rho_div_pi = diffuseColor / PI;
	return rho_div_pi * L1;
}

// http://mimosa-pudica.net/improved-oren-nayar.html
Vector3 ImprovedOrenNayarDiffuse(const Vector3& diffuseColor, float a, const Vector3& N, const Vector3& V, const Vector3& L)
{
	// calculate intermediary values
	float dotNL = saturate(dot(N, L));
	float dotNV = saturate(dot(N, V));
	float dotLV = saturate(dot(L, V));
	Vector3 H = normalize(L + V);
	float dotLH = saturate(dot(L, H));

	float s = dotLV - dotNL * dotNV;
	float t = mix(1.0f, max(dotNL, dotNV), step(0.0f, s));
	float st = s * (1.0f / (t + EPSILON));

	float sigma2 = a;
	Vector3 A = diffuseColor * (0.17f * sigma2 / (sigma2 + 0.13f)) + Vector3(1.0f - 0.5f * sigma2 / (sigma2 + 0.33f));
	float B = 0.45f * sigma2 / (sigma2 + 0.09f);
	//return mulPerElem(diffuseColor * max(0.0f, dotNL), (A + Vector3(B * s / t) / PI));
	return mulPerElem(diffuseColor, (A * RECIPROCAL_PI + Vector3(B * st * RECIPROCAL_PI)));
}

// http://mimosa-pudica.net/improved-oren-nayar.html
Vector3 ImprovedFastOrenNayarDiffuse(const Vector3& diffuseColor, float a, const Vector3& N, const Vector3& V, const Vector3& L)
{
	// calculate intermediary values
	float dotNL = saturate(dot(N, L));
	float dotNV = saturate(dot(N, V));
	float dotLV = saturate(dot(L, V));
	Vector3 H = normalize(L + V);
	float dotLH = saturate(dot(L, H));

	float s = dotLV - dotNL * dotNV;
	float t = mix(1.0f, max(dotNL, dotNV), step(0.0f, s));
	float st = s * (1.0f / (t + EPSILON));

	float A = 1.0f / ((PI * 0.5f - 2.0f / 3.0f) * a + PI);
	float B = a * A;
	return diffuseColor * (A + B * st);
}