Oren-Nayar Quantitative model - From the original paper

oren_nayar_diffuse.glsl

#define EPS    1e-4
#define PI     3.14159265
#define INV_PI 0.31830988

float maxEps(float x)  { return max(EPS, x);        }
float clamp01(float x) { return clamp(x, 0.0, 1.0); }

float pow2(float x) { return x*x; }
vec3 pow2(vec3 x)   { return x*x; }

// OREN-NAYAR MODEL - QUANTITATIVE
// http://www1.cs.columbia.edu/CAVE/publications/pdfs/Oren_CVPR93.pdf
vec3 orenNayarQuantitative(vec3 N, vec3 V, vec3 L, Material material) {
    float alpha = pow2(material.roughness);

    float NdotV = maxEps(dot(N, V));
    float NdotL = maxEps(dot(N, L));

    vec2 angles = acos(vec2(NdotL, NdotV));
    if(angles.x < angles.y) angles = angles.yx;
    float cosA = clamp01(dot(normalize(V - NdotV * N), normalize(L - NdotL * N)));

    float C1   = 1.0 + (alpha / (alpha + 0.4)) * (2.0 / PI - 1.0);

    float tmp = cosA >= 0.0 ? 0.0 : pow((2.0 * angles.y) * INV_PI, 2.0 + (1.0 / 2.0 * alpha));
    float C2 = 0.4 * (alpha / (alpha + 0.2)) * (pow(sin(angles.x), 0.56 + (0.2 / alpha)) - tmp);

    float C3 = 0.11 * (alpha / (alpha + 0.2)) * pow(4.0 * angles.x * angles.y / pow2(PI), 1.55 + (1.0 / 1.55 * alpha));

    vec3 Lr1 = (material.albedo * INV_PI) * NdotL * (C1 + (cosA * C2 * tan(angles.y)) + ((1.0 - abs(cosA)) * C3 * tan((angles.x + angles.y) * 0.5)));
    vec3 Lr2 = pow2(material.albedo) * NdotL * (0.045 * (alpha / (alpha + 0.2))) * ((1.0 - cosA) * pow((2.0 * angles.y) * INV_PI, 1.3 + (0.31 / alpha)));
    return clamp01(Lr1 + Lr2);
}