/**
 * Copyright (C) 2013 Stellar Jockeys.
 * All Rights Reserved.
 *
 * NOTICE: All information contained herein is, and remains the property of
 * Stellar Jockeys and its suppliers, if any.  The intellectual and technical
 * concepts contained herein are proprietary to Stellar Jockeys and its suppliers.
 * Dissemination of this information or reproduction of this material is
 * strictly forbidden unless prior written permission is obtained from
 * Stellar Jockeys.
 *
 * @file fraglightpoint.glsl
 * @date 07/25/13
 */
#ifndef FRAGLIGHTPOINT_GLSL
#define FRAGLIGHTPOINT_GLSL

#ifndef FRAGNORMAL_GLSL
#error fragnormal.glsl should be included
#endif

#ifndef FRAGMATERIAL_GLSL
#error fragmaterial.glsl should be included
#endif

#ifdef FRAGSHADOW_GLSL
float Shadowing(in vec2 texcoord, in float depth, in vec3 lightdir)
{
    const int SAMPLES = 4;
    vec3 WEIGHTS[SAMPLES] = vec3[](vec3(0.25 , 1.0, 1.0),
                                   vec3(0.5  , 1.0, 1.0),
                                   vec3(0.625, 1.0, 0.75),
                                   vec3(0.75 , 1.0, 0.5));

    vec3 offsetDir = lightdir * uProjectionScale * 0.5;
    vec3 fragCoord = vec3 (texcoord, depth);
    float thickness = -uProjectionScale.z;
    float shadowing = 1.0;

    for (int index = 0; index < SAMPLES; ++index)
    {
        float check = ShadowCheck(offsetDir, fragCoord, thickness, WEIGHTS[index], uDepth, uNormals);
        shadowing = min(shadowing, check);
    }

    return shadowing;
}
#endif // FRAGSHADOW_GLSL

void main()
{
    vec2 textureCoord = vScreenPosition * 0.5 + 0.5;

    // Un-normalize frag position
    float depth = texture2D(uDepth, textureCoord).r;
    vec3 fragPosition = vec3(vScreenPosition, depth * 2.0 - 1.0);
    fragPosition.z -= uProjectionBias;
    fragPosition /= uProjectionScale;

    vec3 lUnnormalized = vLightPosition - fragPosition;
    float distSq = dot(lUnnormalized, lUnnormalized);
    float lightRadiusSq = vLightRadius * vLightRadius;

    // Volmetric
    vec3 scattering = vec3(0.0);
#if 1
    {
        // Let r be a point on the projection line relative to the light.
        //
        // We are solving for r.z on the surface of the cone-sphere intersection,
        // but leave the fragment depth as it is used to further prune the results.
        vec3 r = -lUnnormalized.xyz;

        // Use the bounding sphere for the starting values.
        // Degenerates when out of bounds. Set front and back to the same (zero thickness).
        float rxy_rxy = dot(r.xy, r.xy);
        float rz_rz_sphere = lightRadiusSq - rxy_rxy;
        float rz_front = sqrt(max(rz_rz_sphere, 0.0));
        float rz_back = -rz_front;

        float maxThickness = vLightRadius * 2.0;

        // NOTE: This branch is the same for a given light.
        // NOTE: There are numerical issue in data loading/processing which cause
        // slightly negative results when it should actual be zero....
        if (vLightWidth > -0.001)
        {
            maxThickness = vLightRadius;

            // Let d be the unit vector of the direction of the cone.
            float dz_dz = vLightDirection.z * vLightDirection.z;

            // Let w be the cosine of the half-aperture of the cone.
            float w_w = vLightWidth * vLightWidth;

            // Solve for r.z where:
            // dot(d, r) = ||r|| * w
            // dot(d.xy, r.xy) + (d.z * r.z) = sqrt(dot(r, r)) * w
            //
            // Sqaure both sides. Note that this introduces solutions on the negative cone.
            // dot(d.xy, r.xy)^2 + 2(d.z * r.z)dot(d.xy, r.xy) + (d.z * r.z)^2 = (dot(r.xy, r.xy) + r.z^2) * w^2
            //
            // Let a = d.z^2 - w^2
            // Let b = 2 * d.z * dot(d.xy, r.xy)
            // Let c = dot(d.xy, r.xy)^2 - w^2 * dot(r.xy, r.xy)
            float a = dz_dz - w_w;
            float dxy_rxy = dot(vLightDirection.xy, r.xy);

            // NOTE: |d.z| = |w| regresses to linear.
            // An epsilon for ease of testing the case.
            float epsilon = 0.0001;

            // NOTE: This branch is the same for a given light.
            if (a > epsilon)
            {
                // Case 1: Double cone strictly contains the projection vector (Z axis)
                // There are always two real solutions. One on the positive cone and
                // one on the negative cone.
                //
                // Since a is positive, the positive cone solutions use the same sign as d.z.
                // A positive d.z is a back solution. A negative d.z is a front solution.

                //scattering *= vec3(1.0, 0.0, 0.0);

                float discriminant = (dxy_rxy * dxy_rxy * w_w) + (rxy_rxy * a * w_w);
                float sqrt_discriminant = sqrt(discriminant);

                // NOTE: This branch is the same for a given light.
                if (vLightDirection.z > 0.0)
                {
                    float rz = (-vLightDirection.z * dxy_rxy + sqrt_discriminant) / a;
                    rz_back = max(rz_back, rz);
                }
                else
                {
                    float rz = (-vLightDirection.z * dxy_rxy - sqrt_discriminant) / a;
                    rz_front = min(rz_front, rz);
                }
            }
            else if (a < -epsilon)
            {
                // Case 2: Double cone strictly excludes the projection vector (Z axis)
                // If discriminant is negative, the line does not intersect the double cone.
                // Otherwise, the two solutions are both on the positive cone or both on
                // the negative cone.
                //
                // Since a is negative, the negative solution is the front solution.

                //scattering *= vec3(0.0, 1.0, 0.0);

                // TODO(harry): This branch is fragment based, but spatially coherent...
                // Could avoid checking with better light bounding (at least an OBB).
                if (dxy_rxy > 0.0)
                {
                    float discriminant = (dxy_rxy * dxy_rxy * w_w) + (rxy_rxy * a * w_w);

                    // Have degenerate solutions regress to the tangent solution (zero thickness).
                    float sqrt_discriminant = sqrt(max(discriminant, 0.0));

                    rz_front = min(rz_front, (-vLightDirection.z * dxy_rxy - sqrt_discriminant) / a);
                    rz_back = max(rz_back, (-vLightDirection.z * dxy_rxy + sqrt_discriminant) / a);
                }
                else
                {
                    // Solutions would be on negative cone, imaginary or the apex.
                    discard;
                }
            }
            else
            {
                // Case 3: The projection vector (Z axis) lies on the surface of the cone.
                // If d.xy and r.xy are perpendicular there is no solution or it is the apex.
                // Otherwise, there is one solution on either the positive or negative cone.
                //
                // If d.z is positive, the solution is a back solution.
                // If d.z is negative, the solution is a front solution.
                // If d.z is zero, the cone degenerated into a hemisphere.

                //scattering *= vec3(0.0, 0.0, 1.0);

                // TODO(harry): This branch is fragment based, but spatially coherent...
                // Could avoid checking with better light bounding (at least an OBB).
                if (dxy_rxy > 0.0)
                {
                    // NOTE: This branch is the same for a given light.
                    if (vLightDirection.z > 0.0)
                    {
                        float minus_c = (rxy_rxy * w_w) - (dxy_rxy * dxy_rxy);
                        float b = 2.0 * vLightDirection.z * dxy_rxy;
                        float rz = minus_c / b;

                        rz_back = max(rz_back, rz);
                    }
                    else if (vLightDirection.z < 0.0)
                    {
                        float minus_c = (rxy_rxy * w_w) - (dxy_rxy * dxy_rxy);
                        float b = 2.0 * vLightDirection.z * dxy_rxy;
                        float rz = minus_c / b;

                        rz_front = min(rz_front, rz);
                    }
                    else
                    {
                        // There are no conical solutions. Leave it as is.
                    }
                }
                else
                {
                    // Solutions would be on negative cone, imaginary or the apex.
                    discard;
                }
            }
        } // end if conical

        rz_front = max(rz_front, r.z);
        rz_back = max(rz_back, r.z);

        if (rz_front <= rz_back)
        {
            discard;
        }

        float thickness = rz_front - rz_back;
        thickness /= maxThickness;

        vec3 middle = -vec3(r.xy, (rz_front + rz_back) * 0.5);
        float middle_dist_sq = dot(middle, middle);

        // Inverse square-ish falloff
        float density = vVolumetricRadialFalloff * (lightRadiusSq * vVolumetricRadialFalloff + lightRadiusSq);
        density /= lightRadiusSq * vVolumetricRadialFalloff + middle_dist_sq;
        density -= vVolumetricRadialFalloff;
        density *= vLightColor.a * uVolumetricColor.a * uVolumetricIntensity * vVolumetricIntensity;

        vec3 middle_norm = normalize(middle);

        // Conical falloff
        float conicalDensity = (dot(middle_norm, vLightDirection) + vLightWidth) * -vVolumetricConicalFalloff;
        conicalDensity = clamp(conicalDensity, 0.0, 1.0);

        thickness *= clamp(density * conicalDensity, 0.0, 1.0);
        scattering = thickness * vLightColor.rgb * uVolumetricColor.rgb;
    } // end volumetric
#endif

    // Inverse square-ish falloff
    float intensity = vRadialFalloff * (lightRadiusSq * vRadialFalloff + lightRadiusSq);
    intensity /= lightRadiusSq * vRadialFalloff + distSq;
    intensity -= vRadialFalloff;
    intensity = max(intensity * vLightColor.a * uLightIntensity, 0.0);

    const vec3 v = vec3(0.0, 0.0, 1.0);
    vec3 l = normalize(lUnnormalized);
    vec3 n = ColorToNormal(texture2D(uNormals, textureCoord).rgb);
    vec3 h = normalize(l + v);

#ifdef FRAGSHADOW_GLSL
    vec3 shadowL = mix(lUnnormalized, l, step(1.0, distSq));
    intensity *= Shadowing(textureCoord, depth, shadowL);
#endif

    // Conical falloff
    float conicalFalloff = (dot(l, vLightDirection) + vLightWidth) * -vConicalFalloff;
    conicalFalloff = clamp(conicalFalloff, 0.0, 1.0);

    vec4 material = texture2D(uMaterial, textureCoord);
    float fresnel = UnpackFresnel(material);
    float specExp = UnpackSpecularExponent(material);
    float shading = max(dot(l, n), 0.0) * conicalFalloff;
    float specular = Specularity(n, h, specExp, fresnel);

    // Diffuse albedo and specular
    vec4 color = texture2D(uDiffuse, textureCoord);
    color.rgb = color.rgb * uDiffuseAmount + specular;

    // Light properties and shading
    color.rgb *= vLightColor.rgb * (intensity * shading);
    color.rgb += scattering;

    gl_FragColor = vec4(color.rgb, dot(color.rgb, DESATURATE));
}

#endif // FRAGLIGHTPOINT_GLSL
