#version 330

// Linear interpolated input values per fragment

in vec3 vertexEye;
in vec3 normalEye;

// Uniform values for shader selection

// 0: None
// 1: Lambert
// 2: OrenNayar
// 3: Schlick
// 4: Sandford
uniform int uniDiffuseShader;

// 0: None
// 1: BlinnPhong
// 2: BlinnPhongNormalized
// 3: CookTorrance
// 4: Ward
// 5: Sandford
uniform int uniSpecularShader;

// Uniform values for material properties

uniform vec4  uniMaterialEmission;
uniform vec4  uniMaterialAmbient;
uniform vec4  uniMaterialDiffuse;
uniform vec4  uniMaterialSpecular;
uniform float uniMaterialShininess;
uniform float uniMaterialRoughness;

// Uniform values for lightsource properties

uniform vec4 uniLightSourcePosition;
uniform vec4 uniLightSourceAmbient;
uniform vec4 uniLightSourceDiffuse;
uniform vec4 uniLightSourceSpecular;

// Output value is fragment color
layout(location = 0, index = 0) out vec4 fragColor;

float PI = 3.1415926535;

// Utility function for Sandford model
float g( float theta, float b )
{
    return 1.0 / ( 1.0 + pow( b, 2.0 ) * pow( tan( theta ), 2.0 ) );
}

vec4 Emission()
{
    return uniMaterialEmission;
}

vec4 Ambient()
{
    return uniMaterialAmbient * uniLightSourceAmbient;
}

vec4 Diffuse( vec3 N, vec3 L, float NdotL )
{
    vec4 diffuse = vec4( 0.0, 0.0, 0.0, 1.0 );

    if( uniDiffuseShader == 1 ) {
        // Lambert
        if( NdotL > 0 ) {
            diffuse = NdotL * uniMaterialDiffuse * uniLightSourceDiffuse;
        }
    } else if( uniDiffuseShader == 2 ) {
        // OrenNayar
        if( NdotL > 0 ) {
            vec3 E = -normalize( vertexEye );
            float NdotE   = max( dot( N, E ), 0.0 );
            float angleNE = acos( NdotE );
            float angleNL = acos( NdotL );
            float alpha = max( angleNE, angleNL );
            float beta = min( angleNE, angleNL );
            float gamma = max( cos( angleNE - angleNL ), 0.0 );
            float roughness = uniMaterialRoughness * uniMaterialRoughness;
            float A = 1.0 - 0.5 * ( roughness / ( roughness + 0.33 ) );
            float B = 0.45 * ( roughness / ( roughness + 0.09 ) );
            float OrenNayar = NdotL * ( A + ( B * gamma * sin( alpha ) * tan( beta ) ) );
            diffuse = OrenNayar * uniMaterialDiffuse * uniLightSourceDiffuse;
        }
    } else if( uniDiffuseShader == 3 ) {
        // Schlick
        if( NdotL > 0 ) {
            float spec = ( uniMaterialSpecular.r + uniMaterialSpecular.g + uniMaterialSpecular.b ) / 3.0;
            float OneMinusNdotL = 1.0 - NdotL;
            float Schlick = spec + ( 1.0 - spec ) * pow( OneMinusNdotL, 5.0 );
            diffuse = NdotL * ( 1.0 - Schlick ) * uniMaterialDiffuse * uniLightSourceDiffuse;
        }
    } else if( uniDiffuseShader == 4 ) {
        // Sandford
        // I_d = g(theta_r) * rho_d(lamda) * g(theta_i) * max(cos(theta_i), 0.0) * I_0/ (Pi^2 * G^2(b))
        vec3 E = -normalize( vertexEye );
        float NdotE = max( dot( N, E ), 0.0 );

        if( NdotL > 0 && NdotE > 0 ) {
            float b = uniMaterialRoughness;
            float e = uniMaterialShininess;
            vec4 Nominator = g( acos( NdotE ), b) * g( acos( NdotL ), b ) *
                             NdotL * uniMaterialDiffuse * uniLightSourceDiffuse; 

            float b2 = pow( b, 2.0 );
            float Denom = 1.0 / ( 1.0 - b2 ) * (1.0 - b2 / ( 1.0 - b2 ) * log( 1.0 / b2 ) );

            diffuse = Nominator / Denom;
        }
    }
    return diffuse;
}

vec4 Specular( vec3 N, vec3 L, float NdotL )
{
    vec4 specular = vec4( 0.0, 0.0, 0.0, 1.0 );

    if( uniSpecularShader == 1 ) {
        // BlinnPhong
        if( NdotL > 0 ) {
            vec3 E = -normalize( vertexEye );
            vec3 H = normalize( E + L );
            float spec = pow( max( dot( H, N ), 0 ), uniMaterialShininess );
            specular = spec * uniMaterialSpecular * uniLightSourceSpecular;
        }
    } else if( uniSpecularShader == 2 ) {
        // BlinnPhongNormalized
        if( NdotL > 0 ) {
            vec3 E = -normalize( vertexEye );
            vec3 H = normalize( E + L );
            float norm = NdotL * ( uniMaterialShininess + 8 ) / 160;
            float spec = pow( max( dot( H, N ), 0 ), uniMaterialShininess );
            specular = spec * norm * uniMaterialSpecular * uniLightSourceSpecular;
        }
    } else if( uniSpecularShader == 3 ) {
        // CookTorrance
        if( NdotL > 0 ) {
            vec3 E = -normalize( vertexEye );
            vec3 H = normalize( E + L );

            float LdotH = max(dot(L, H),0.0);
            float NdotE = max(dot(N, E),0.0);
            float NdotH = max(dot(N, H),0.0);
            float NdotL = max(dot(N, L),0.0);
            if( LdotH > 0 && NdotH > 0 && NdotE > 0 ) {
                float spec = ( uniMaterialSpecular.r + uniMaterialSpecular.g + uniMaterialSpecular.b ) / 3.0;
		float OneMinusLdotH = 1 - LdotH;
		float F = spec + ( 1.0 - spec ) * pow( OneMinusLdotH, 5.0 );
		float g1 = 2.0 * NdotH * NdotE / LdotH;
		float g2 = 2.0 * NdotH * NdotL / LdotH;
		float geoAtt = min( 1.0, min( g1, g2 ) );
                float roughness = uniMaterialRoughness * uniMaterialRoughness;
		float r1 = 1.0 / ( roughness * pow( NdotH, 4.0 ) );
		float r2 = ( NdotH * NdotH - 1.0 ) / ( roughness * NdotH * NdotH );
		float NDF = r1 * exp( r2 );
		float CT_factor = F * geoAtt * NDF / NdotE;
		specular = CT_factor * uniMaterialSpecular * uniLightSourceSpecular / (4.0 * PI );
            }
        }
    } else if( uniSpecularShader == 4 ) {
        // Ward
        if( NdotL > 0 ) {
            vec3 E = -normalize( vertexEye );
            vec3 H = normalize( E + L );

            float NdotE = max( dot( N, E ), 0.0 );
            float NdotH = max( dot( N, H ), 0.0 );
            float roughness = uniMaterialRoughness * uniMaterialRoughness;
            float A = 1.0 / ( sqrt( NdotL * NdotE ) );
            float B = 1.0 / ( 4.0 * PI * roughness );
            float angle_delta = acos( NdotH );
            float tan_d_alpha2 = tan( angle_delta / roughness );
            float C = exp( -1.0 * tan_d_alpha2 * tan_d_alpha2 );
            float IsoWard = A * B * C;
            specular = NdotL * IsoWard * uniMaterialSpecular * uniLightSourceSpecular;
        }
    } else if( uniSpecularShader == 5 ) {
        // Sandford
        // I_s = h(alpha) * max(cos(theta_i), 0.0) * rho_s(theta_i) * I_0/ (4 * Pi * H(theta_i)
        vec3 E = -normalize( vertexEye );
        vec3 H = normalize( E + L );

        float NdotE = max( dot( N, E ), 0.0 );
        float NdotH = max( dot( N, H ), 0.0 );

        if( NdotL > 0 && NdotH > 0 && NdotE > 0 ) {
            float b = uniMaterialRoughness;
            float e = uniMaterialShininess;

            float e2    = pow( e, 2.0 );
            float e3    = pow( 1 - e2, 2.0 );
            float cos2  = pow( NdotL, 2.0 );
            float HNorm = 1.0 / (2.0 * e2) * 
                          ( ( 1.0 - e2 ) * NdotL + ( 2.0*e2 + e3*cos2 ) / sqrt( e3*cos2 + 4.0*e2 ) );

            float h = 1.0 / pow( pow( e * NdotH, 2.0 ) * ( 1.0 - pow( NdotH, 2.0 ) ), 2.0 );

            // Correction 1/(10000 * 2 * Pi), because area A must be converted from m2 to cm2.
            specular = h * NdotL * uniMaterialSpecular * uniLightSourceSpecular /
                       (10000.0 * 8.0 * PI * PI * HNorm );
            // Add limit factor as proposed in Sanford paper, but use b == 1.
            // Needed to avoid halos at the edges of objects.
            // * 1 / (1 + b^2 * tan^2(theta_r))^2
            float limitFactor = pow( g( acos( NdotE ), 1.0 ), 2.0 );
            specular *= limitFactor;
        }
    }
    return specular;
}

void main()
{
    vec3 N = normalize( normalEye );
    vec3 L = normalize( vec3( uniLightSourcePosition ) - vertexEye );
    float NdotL = max( dot( N, L ), 0.0 );

    fragColor = Emission() +
                Ambient()  +
                Diffuse(  N, L, NdotL ) +
                Specular( N, L, NdotL );
}

