/*
 * Copyright (C) 2015 The Android Open Source Project
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */


#include <ibl/CubemapIBL.h>

#include <ibl/Cubemap.h>
#include <ibl/CubemapUtils.h>
#include <ibl/utilities.h>

#include "CubemapUtilsImpl.h"

#include <utils/JobSystem.h>

#include <math/mat3.h>
#include <math/scalar.h>

#include <random>
#include <vector>

using namespace filament::math;
using namespace utils;

namespace filament {
namespace ibl {

static float pow5(float x) {
    const float x2 = x * x;
    return x2 * x2 * x;
}

static float pow6(float x) {
    const float x2 = x * x;
    return x2 * x2 * x2;
}

static float3 hemisphereImportanceSampleDggx(float2 u, float a) { // pdf = D(a) * cosTheta
    const float phi = 2.0f * (float) F_PI * u.x;
    // NOTE: (aa-1) == (a-1)(a+1) produces better fp accuracy
    const float cosTheta2 = (1 - u.y) / (1 + (a + 1) * ((a - 1) * u.y));
    const float cosTheta = std::sqrt(cosTheta2);
    const float sinTheta = std::sqrt(1 - cosTheta2);
    return { sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta };
}

static float3 UTILS_UNUSED hemisphereCosSample(float2 u) {  // pdf = cosTheta / F_PI;
    const float phi = 2.0f * (float) F_PI * u.x;
    const float cosTheta2 = 1 - u.y;
    const float cosTheta = std::sqrt(cosTheta2);
    const float sinTheta = std::sqrt(1 - cosTheta2);
    return { sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta };
}

static float3 UTILS_UNUSED hemisphereUniformSample(float2 u) { // pdf = 1.0 / (2.0 * F_PI);
    const float phi = 2.0f * (float) F_PI * u.x;
    const float cosTheta = 1 - u.y;
    const float sinTheta = std::sqrt(1 - cosTheta * cosTheta);
    return { sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta };
}

/*
 *
 * Importance sampling Charlie
 * ---------------------------
 *
 * In order to pick the most significative samples and increase the convergence rate, we chose to
 * rely on Charlie's distribution function for the pdf as we do in hemisphereImportanceSampleDggx.
 *
 * To determine the direction we then need to resolve the cdf associated to the chosen pdf for random inputs.
 *
 * Knowing pdf() = DCharlie(h) <n•h>
 *
 * We need to find the cdf:
 *
 * / 2pi     / pi/2
 * |         |  (2 + (1 / a)) * sin(theta) ^ (1 / a) * cos(theta) * sin(theta)
 * / phi=0   / theta=0
 *
 * We sample theta and phi independently.
 *
 * 1. as in all the other isotropic cases phi = 2 * pi * epsilon
 *    (https://www.tobias-franke.eu/log/2014/03/30/notes_on_importance_sampling.html)
 *
 * 2. we need to solve the integral on theta:
 *
 *             / sTheta
 * P(sTheta) = |  (2 + (1 / a)) * sin(theta) ^ (1 / a + 1) * cos(theta) * dtheta
 *             / theta=0
 *
 * By subsitution of u = sin(theta) and du = cos(theta) * dtheta
 *
 * /
 * |  (2 + (1 / a)) * u ^ (1 / a + 1) * du
 * /
 *
 * = (2 + (1 / a)) * u ^ (1 / a + 2) / (1 / a + 2)
 *
 * = u ^ (1 / a + 2)
 *
 * = sin(theta) ^ (1 / a + 2)
 *
 *             +-                          -+ sTheta
 * P(sTheta) = |  sin(theta) ^ (1 / a + 2)  |
 *             +-                          -+ 0
 *
 * P(sTheta) = sin(sTheta) ^ (1 / a + 2)
 *
 * We now need to resolve the cdf for an epsilon value:
 *
 * epsilon = sin(theta) ^ (a / ( 2 * a + 1))
 *
 *  +--------------------------------------------+
 *  |                                            |
 *  |  sin(theta) = epsilon ^ (a / ( 2 * a + 1)) |
 *  |                                            |
 *  +--------------------------------------------+
 */
static float3 UTILS_UNUSED hemisphereImportanceSampleDCharlie(float2 u, float a) { // pdf = DistributionCharlie() * cosTheta
    const float phi = 2.0f * (float) F_PI * u.x;

    const float sinTheta = std::pow(u.y, a / (2 * a + 1));
    const float cosTheta = std::sqrt(1 - sinTheta * sinTheta);

    return { sinTheta * std::cos(phi), sinTheta * std::sin(phi), cosTheta };
}

static float DistributionGGX(float NoH, float linearRoughness) {
    // NOTE: (aa-1) == (a-1)(a+1) produces better fp accuracy
    float a = linearRoughness;
    float f = (a - 1) * ((a + 1) * (NoH * NoH)) + 1;
    return (a * a) / ((float) F_PI * f * f);
}

static float UTILS_UNUSED DistributionAshikhmin(float NoH, float linearRoughness) {
    float a = linearRoughness;
    float a2 = a * a;
    float cos2h = NoH * NoH;
    float sin2h = 1 - cos2h;
    float sin4h = sin2h * sin2h;
    return 1.0f / ((float) F_PI * (1 + 4 * a2)) * (sin4h + 4 * std::exp(-cos2h / (a2 * sin2h)));
}

static float UTILS_UNUSED DistributionCharlie(float NoH, float linearRoughness) {
    // Estevez and Kulla 2017, "Production Friendly Microfacet Sheen BRDF"
    float a = linearRoughness;
    float invAlpha = 1 / a;
    float cos2h = NoH * NoH;
    float sin2h = 1 - cos2h;
    return (2.0f + invAlpha) * std::pow(sin2h, invAlpha * 0.5f) / (2.0f * (float) F_PI);
}

static float Fresnel(float f0, float f90, float LoH) {
    const float Fc = pow5(1 - LoH);
    return f0 * (1 - Fc) + f90 * Fc;
}

static float Visibility(float NoV, float NoL, float a) {
    // Heitz 2014, "Understanding the Masking-Shadowing Function in Microfacet-Based BRDFs"
    // Height-correlated GGX
    const float a2 = a * a;
    const float GGXL = NoV * std::sqrt((NoL - NoL * a2) * NoL + a2);
    const float GGXV = NoL * std::sqrt((NoV - NoV * a2) * NoV + a2);
    return 0.5f / (GGXV + GGXL);
}

static float UTILS_UNUSED VisibilityAshikhmin(float NoV, float NoL, float /*a*/) {
    // Neubelt and Pettineo 2013, "Crafting a Next-gen Material Pipeline for The Order: 1886"
    return 1 / (4 * (NoL + NoV - NoL * NoV));
}

/*
 *
 * Importance sampling GGX - Trowbridge-Reitz
 * ------------------------------------------
 *
 * Important samples are chosen to integrate Dggx() * cos(theta) over the hemisphere.
 *
 * All calculations are made in tangent space, with n = [0 0 1]
 *
 *             l        h (important sample)
 *             .\      /.
 *             . \    / .
 *             .  \  /  .
 *             .   \/   .
 *         ----+---o----+-------> n [0 0 1]
 *     cos(2*theta)     cos(theta)
 *        = n•l            = n•h
 *
 *  v = n
 *  f0 = f90 = 1
 *  V = 1
 *
 *  h is micro facet's normal
 *
 *  l is the reflection of v (i.e.: n) around h  ==>  n•h = l•h = v•h
 *
 *  h = important_sample_ggx()
 *
 *  n•h = [0 0 1]•h = h.z
 *
 *  l = reflect(-n, h)
 *    = 2 * (n•h) * h - n;
 *
 *  n•l = cos(2 * theta)
 *      = cos(theta)^2 - sin(theta)^2
 *      = (n•h)^2 - (1 - (n•h)^2)
 *      = 2(n•h)^2 - 1
 *
 *
 *  pdf() = D(h) <n•h> |J(h)|
 *
 *               1
 *  |J(h)| = ----------
 *            4 <v•h>
 *
 *
 * Pre-filtered importance sampling
 * --------------------------------
 *
 *  see: "Real-time Shading with Filtered Importance Sampling", Jaroslav Krivanek
 *  see: "GPU-Based Importance Sampling, GPU Gems 3", Mark Colbert
 *
 *
 *                   Ωs
 *     lod = log4(K ----)
 *                   Ωp
 *
 *     log4(K) = 1, works well for box filters
 *     K = 4
 *
 *             1
 *     Ωs = ---------, solid-angle of an important sample
 *           N * pdf
 *
 *              4 PI
 *     Ωp ~ --------------, solid-angle of a sample in the base cubemap
 *           texel_count
 *
 *
 * Evaluating the integral
 * -----------------------
 *
 *                    K     fr(h)
 *            Er() = --- ∑ ------- L(h) <n•l>
 *                    N  h   pdf
 *
 * with:
 *
 *            fr() = D(h)
 *
 *                       N
 *            K = -----------------
 *                    fr(h)
 *                 ∑ ------- <n•l>
 *                 h   pdf
 *
 *
 *  It results that:
 *
 *            K           4 <v•h>
 *    Er() = --- ∑ D(h) ------------ L(h) <n•l>
 *            N  h        D(h) <n•h>
 *
 *
 *              K
 *    Er() = 4 --- ∑ L(h) <n•l>
 *              N  h
 *
 *                  N       4
 *    Er() = ------------- --- ∑ L(v) <n•l>
 *             4 ∑ <n•l>    N
 *
 *
 *  +------------------------------+
 *  |          ∑ <n•l> L(h)        |
 *  |  Er() = --------------       |
 *  |            ∑ <n•l>           |
 *  +------------------------------+
 *
 */

UTILS_ALWAYS_INLINE
void CubemapIBL::roughnessFilter(
        utils::JobSystem& js, Cubemap& dst, const std::vector<Cubemap>& levels,
        float linearRoughness, size_t maxNumSamples, math::float3 mirror, bool prefilter,
        Progress updater, void* userdata) {
    roughnessFilter(js, dst, { levels.data(), uint32_t(levels.size()) },
            linearRoughness, maxNumSamples, mirror, prefilter, updater, userdata);
}

void CubemapIBL::roughnessFilter(
        utils::JobSystem& js, Cubemap& dst, const utils::Slice<Cubemap>& levels,
        float linearRoughness, size_t maxNumSamples, math::float3 mirror, bool prefilter,
        Progress updater, void* userdata)
{
    const float numSamples = maxNumSamples;
    const float inumSamples = 1.0f / numSamples;
    const size_t maxLevel = levels.size()-1;
    const float maxLevelf = maxLevel;
    const Cubemap& base(levels[0]);
    const size_t dim0 = base.getDimensions();
    const float omegaP = (4.0f * (float) F_PI) / float(6 * dim0 * dim0);
    std::atomic_uint progress = {0};

    if (linearRoughness == 0) {
        auto scanline = [&]
                (CubemapUtils::EmptyState&, size_t y, Cubemap::Face f, Cubemap::Texel* data, size_t dim) {
                    if (UTILS_UNLIKELY(updater)) {
                        size_t p = progress.fetch_add(1, std::memory_order_relaxed) + 1;
                        updater(0, (float)p / ((float) dim * 6.0f), userdata);
                    }
                    const Cubemap& cm = levels[0];
                    for (size_t x = 0; x < dim; ++x, ++data) {
                        const float2 p(Cubemap::center(x, y));
                        const float3 N(dst.getDirectionFor(f, p.x, p.y) * mirror);
                        // FIXME: we should pick the proper LOD here and do trilinear filtering
                        Cubemap::writeAt(data, cm.sampleAt(N));
                    }
        };
        // at least 256 pixel cubemap before we use multithreading -- the overhead of launching
        // jobs is too large compared to the work above.
        if (dst.getDimensions() <= 256) {
            CubemapUtils::processSingleThreaded<CubemapUtils::EmptyState>(
                    dst, js, std::ref(scanline));
        } else {
            CubemapUtils::process<CubemapUtils::EmptyState>(dst, js, std::ref(scanline));
        }
        return;
    }

    // be careful w/ the size of this structure, the smaller the better
    struct CacheEntry {
        float3 L;
        float brdf_NoL;
        float lerp;
        uint8_t l0;
        uint8_t l1;
    };

    std::vector<CacheEntry> cache;
    cache.reserve(maxNumSamples);

    // precompute everything that only depends on the sample #
    float weight = 0;
    // index of the sample to use
    // our goal is to use maxNumSamples for which NoL is > 0
    // to achieve this, we might have to try more samples than
    // maxNumSamples
    for (size_t sampleIndex = 0 ; sampleIndex < maxNumSamples; sampleIndex++) {

        // get Hammersley distribution for the half-sphere
        const float2 u = hammersley(uint32_t(sampleIndex), inumSamples);

        // Importance sampling GGX - Trowbridge-Reitz
        const float3 H = hemisphereImportanceSampleDggx(u, linearRoughness);

#if 0
        // This produces the same result that the code below using the the non-simplified
        // equation. This let's us see that N == V and that L = -reflect(V, H)
        // Keep this for reference.
        const float3 N = {0, 0, 1};
        const float3 V = N;
        const float3 L = 2 * dot(H, V) * H - V;
        const float NoL = dot(N, L);
        const float NoH = dot(N, H);
        const float NoH2 = NoH * NoH;
        const float NoV = dot(N, V);
#else
        const float NoH = H.z;
        const float NoH2 = H.z * H.z;
        const float NoL = 2 * NoH2 - 1;
        const float3 L(2 * NoH * H.x, 2 * NoH * H.y, NoL);
#endif

        if (NoL > 0) {
            const float pdf = DistributionGGX(NoH, linearRoughness) / 4;

            // K is a LOD bias that allows a bit of overlapping between samples
            constexpr float K = 4;
            const float omegaS = 1 / (numSamples * pdf);
            const float l = float(log4(omegaS) - log4(omegaP) + log4(K));
            const float mipLevel = prefilter ? clamp(float(l), 0.0f, maxLevelf) : 0.0f;

            const float brdf_NoL = float(NoL);

            weight += brdf_NoL;

            uint8_t l0 = uint8_t(mipLevel);
            uint8_t l1 = uint8_t(std::min(maxLevel, size_t(l0 + 1)));
            float lerp = mipLevel - (float) l0;

            cache.push_back({ L, brdf_NoL, lerp, l0, l1 });
        }
    }

    for (auto& entry : cache) {
        entry.brdf_NoL *= 1.0f / weight;
    }

    // we can sample the cubemap in any order, sort by the weight, it could improve fp precision
    std::sort(cache.begin(), cache.end(), [](CacheEntry const& lhs, CacheEntry const& rhs) {
        return lhs.brdf_NoL < rhs.brdf_NoL;
    });


    struct State {
        // maybe blue-noise instead would look even better
        std::default_random_engine gen;
        std::uniform_real_distribution<float> distribution{ -F_PI, F_PI };
    };

    auto scanline = [&](State& state, size_t y,
            Cubemap::Face f, Cubemap::Texel* data, size_t dim) {
        if (UTILS_UNLIKELY(updater)) {
            size_t p = progress.fetch_add(1, std::memory_order_relaxed) + 1;
            updater(0, (float) p / ((float) dim * 6.0f), userdata);
        }
        mat3 R;
        const size_t numSamples = cache.size();
        for (size_t x = 0; x < dim; ++x, ++data) {
            const float2 p(Cubemap::center(x, y));
            const float3 N(dst.getDirectionFor(f, p.x, p.y) * mirror);

            // center the cone around the normal (handle case of normal close to up)
            const float3 up = std::abs(N.z) < 0.999 ? float3(0, 0, 1) : float3(1, 0, 0);
            R[0] = normalize(cross(up, N));
            R[1] = cross(N, R[0]);
            R[2] = N;

            R *= mat3f::rotation(state.distribution(state.gen), float3{0,0,1});

            float3 Li = 0;
            for (size_t sample = 0; sample < numSamples; sample++) {
                const CacheEntry& e = cache[sample];
                const float3 L(R * e.L);
                const Cubemap& cmBase = levels[e.l0];
                const Cubemap& next = levels[e.l1];
                const float3 c0 = Cubemap::trilinearFilterAt(cmBase, next, e.lerp, L);
                Li += c0 * e.brdf_NoL;
            }
            Cubemap::writeAt(data, Cubemap::Texel(Li));
        }
    };

    // don't use the jobsystem unless we have enough work per scanline -- or the overhead of
    // launching jobs will prevail.
    if (dst.getDimensions() * maxNumSamples <= 256) {
        CubemapUtils::processSingleThreaded<State>(dst, js, std::ref(scanline));
    } else {
        CubemapUtils::process<State>(dst, js, std::ref(scanline));
    }
}

/*
 *
 * Importance sampling
 * -------------------
 *
 * Important samples are chosen to integrate cos(theta) over the hemisphere.
 *
 * All calculations are made in tangent space, with n = [0 0 1]
 *
 *                      l (important sample)
 *                     /.
 *                    / .
 *                   /  .
 *                  /   .
 *         --------o----+-------> n (direction)
 *                   cos(theta)
 *                    = n•l
 *
 *
 *  'direction' is given as an input parameter, and serves as tge z direction of the tangent space.
 *
 *  l = important_sample_cos()
 *
 *  n•l = [0 0 1] • l = l.z
 *
 *           n•l
 *  pdf() = -----
 *           PI
 *
 *
 * Pre-filtered importance sampling
 * --------------------------------
 *
 *  see: "Real-time Shading with Filtered Importance Sampling", Jaroslav Krivanek
 *  see: "GPU-Based Importance Sampling, GPU Gems 3", Mark Colbert
 *
 *
 *                   Ωs
 *     lod = log4(K ----)
 *                   Ωp
 *
 *     log4(K) = 1, works well for box filters
 *     K = 4
 *
 *             1
 *     Ωs = ---------, solid-angle of an important sample
 *           N * pdf
 *
 *              4 PI
 *     Ωp ~ --------------, solid-angle of a sample in the base cubemap
 *           texel_count
 *
 *
 * Evaluating the integral
 * -----------------------
 *
 * We are trying to evaluate the following integral:
 *
 *                     /
 *             Ed() =  | L(s) <n•l> ds
 *                     /
 *                     Ω
 *
 * For this, we're using importance sampling:
 *
 *                    1     L(l)
 *            Ed() = --- ∑ ------- <n•l>
 *                    N  l   pdf
 *
 *
 *  It results that:
 *
 *             1           PI
 *    Ed() = ---- ∑ L(l) ------  <n•l>
 *            N   l        n•l
 *
 *
 *  To avoid multiplying by 1/PI in the shader, we do it here, which simplifies to:
 *
 *  +----------------------+
 *  |          1           |
 *  |  Ed() = ---- ∑ L(l)  |
 *  |          N   l       |
 *  +----------------------+
 *
 */

void CubemapIBL::diffuseIrradiance(JobSystem& js, Cubemap& dst, const std::vector<Cubemap>& levels,
        size_t maxNumSamples, CubemapIBL::Progress updater, void* userdata)
{
    const float numSamples = maxNumSamples;
    const float inumSamples = 1.0f / numSamples;
    const size_t maxLevel = levels.size()-1;
    const float maxLevelf = maxLevel;
    const Cubemap& base(levels[0]);
    const size_t dim0 = base.getDimensions();
    const float omegaP = (4.0f * (float) F_PI) / float(6 * dim0 * dim0);

    std::atomic_uint progress = {0};

    struct CacheEntry {
        float3 L;
        float lerp;
        uint8_t l0;
        uint8_t l1;
    };

    std::vector<CacheEntry> cache;
    cache.reserve(maxNumSamples);

    // precompute everything that only depends on the sample #
    for (size_t sampleIndex = 0; sampleIndex < maxNumSamples; sampleIndex++) {
        // get Hammersley distribution for the half-sphere
        const float2 u = hammersley(uint32_t(sampleIndex), inumSamples);
        const float3 L = hemisphereCosSample(u);
        const float3 N = { 0, 0, 1 };
        const float NoL = dot(N, L);

        if (NoL > 0) {
            float pdf = NoL * (float) F_1_PI;

            constexpr float K = 4;
            const float omegaS = 1.0f / (numSamples * pdf);
            const float l = float(log4(omegaS) - log4(omegaP) + log4(K));
            const float mipLevel = clamp(float(l), 0.0f, maxLevelf);

            uint8_t l0 = uint8_t(mipLevel);
            uint8_t l1 = uint8_t(std::min(maxLevel, size_t(l0 + 1)));
            float lerp = mipLevel - (float) l0;

            cache.push_back({ L, lerp, l0, l1 });
        }
    }

    CubemapUtils::process<CubemapUtils::EmptyState>(dst, js,
            [&](CubemapUtils::EmptyState&, size_t y,
                    Cubemap::Face f, Cubemap::Texel* data, size_t dim) {

        if (updater) {
            size_t p = progress.fetch_add(1, std::memory_order_relaxed) + 1;
            updater(0, (float)p / ((float) dim * 6.0f), userdata);
        }

        mat3 R;
        const size_t numSamples = cache.size();
        for (size_t x = 0; x < dim; ++x, ++data) {
            const float2 p(Cubemap::center(x, y));
            const float3 N(dst.getDirectionFor(f, p.x, p.y));

            // center the cone around the normal (handle case of normal close to up)
            const float3 up = std::abs(N.z) < 0.999 ? float3(0, 0, 1) : float3(1, 0, 0);
            R[0] = normalize(cross(up, N));
            R[1] = cross(N, R[0]);
            R[2] = N;

            float3 Li = 0;
            for (size_t sample = 0; sample < numSamples; sample++) {
                const CacheEntry& e = cache[sample];
                const float3 L(R * e.L);
                const Cubemap& cmBase = levels[e.l0];
                const Cubemap& next = levels[e.l1];
                const float3 c0 = Cubemap::trilinearFilterAt(cmBase, next, e.lerp, L);
                Li += c0;
            }
            Cubemap::writeAt(data, Cubemap::Texel(Li * inumSamples));
        }
    });
}

// Not importance-sampled
static float2 UTILS_UNUSED DFV_NoIS(float NoV, float roughness, size_t numSamples) {
    float2 r = 0;
    const float linearRoughness = roughness * roughness;
    const float3 V(std::sqrt(1 - NoV * NoV), 0, NoV);
    for (size_t i = 0; i < numSamples; i++) {
        const float2 u = hammersley(uint32_t(i), 1.0f / numSamples);
        const float3 H = hemisphereCosSample(u);
        const float3 L = 2 * dot(V, H) * H - V;
        const float VoH = saturate(dot(V, H));
        const float NoL = saturate(L.z);
        const float NoH = saturate(H.z);
        if (NoL > 0) {
            // Note: remember VoH == LoH  (H is half vector)
            const float J = 1.0f / (4.0f * VoH);
            const float pdf = NoH / (float) F_PI;
            const float d = DistributionGGX(NoH, linearRoughness) * NoL / (pdf * J);
            const float Fc = pow5(1 - VoH);
            const float v = Visibility(NoV, NoL, linearRoughness);
            r.x += d * v * (1.0f - Fc);
            r.y += d * v * Fc;
        }
    }
    return r / numSamples;
}

/*
 *
 * Importance sampling GGX - Trowbridge-Reitz
 * ------------------------------------------
 *
 * Important samples are chosen to integrate Dggx() * cos(theta) over the hemisphere.
 *
 * All calculations are made in tangent space, with n = [0 0 1]
 *
 *                      h (important sample)
 *                     /.
 *                    / .
 *                   /  .
 *                  /   .
 *         --------o----+-------> n
 *                   cos(theta)
 *                    = n•h
 *
 *  h is micro facet's normal
 *  l is the reflection of v around h, l = reflect(-v, h)  ==>  v•h = l•h
 *
 *  n•v is given as an input parameter at runtime
 *
 *  Since n = [0 0 1], we also have v.z = n•v
 *
 *  Since we need to compute v•h, we chose v as below. This choice only affects the
 *  computation of v•h (and therefore the fresnel term too), but doesn't affect
 *  n•l, which only relies on l.z (which itself only relies on v.z, i.e.: n•v)
 *
 *      | sqrt(1 - (n•v)^2)     (sin)
 *  v = | 0
 *      | n•v                   (cos)
 *
 *
 *  h = important_sample_ggx()
 *
 *  l = reflect(-v, h) = 2 * v•h * h - v;
 *
 *  n•l = [0 0 1] • l = l.z
 *
 *  n•h = [0 0 1] • l = h.z
 *
 *
 *  pdf() = D(h) <n•h> |J(h)|
 *
 *               1
 *  |J(h)| = ----------
 *            4 <v•h>
 *
 *
 * Evaluating the integral
 * -----------------------
 *
 * We are trying to evaluate the following integral:
 *
 *                    /
 *             Er() = | fr(s) <n•l> ds
 *                    /
 *                    Ω
 *
 * For this, we're using importance sampling:
 *
 *                    1     fr(h)
 *            Er() = --- ∑ ------- <n•l>
 *                    N  h   pdf
 *
 * with:
 *
 *            fr() = D(h) F(h) V(v, l)
 *
 *
 *  It results that:
 *
 *            1                        4 <v•h>
 *    Er() = --- ∑ D(h) F(h) V(v, l) ------------ <n•l>
 *            N  h                     D(h) <n•h>
 *
 *
 *  +-------------------------------------------+
 *  |          4                  <v•h>         |
 *  |  Er() = --- ∑ F(h) V(v, l) ------- <n•l>  |
 *  |          N  h               <n•h>         |
 *  +-------------------------------------------+
 *
 */

static float2 DFV(float NoV, float linearRoughness, size_t numSamples) {
    float2 r = 0;
    const float3 V(std::sqrt(1 - NoV * NoV), 0, NoV);
    for (size_t i = 0; i < numSamples; i++) {
        const float2 u = hammersley(uint32_t(i), 1.0f / numSamples);
        const float3 H = hemisphereImportanceSampleDggx(u, linearRoughness);
        const float3 L = 2 * dot(V, H) * H - V;
        const float VoH = saturate(dot(V, H));
        const float NoL = saturate(L.z);
        const float NoH = saturate(H.z);
        if (NoL > 0) {
            /*
             * Fc = (1 - V•H)^5
             * F(h) = f0*(1 - Fc) + f90*Fc
             *
             * f0 and f90 are known at runtime, but thankfully can be factored out, allowing us
             * to split the integral in two terms and store both terms separately in a LUT.
             *
             * At runtime, we can reconstruct Er() exactly as below:
             *
             *            4                      <v•h>
             *   DFV.x = --- ∑ (1 - Fc) V(v, l) ------- <n•l>
             *            N  h                   <n•h>
             *
             *
             *            4                      <v•h>
             *   DFV.y = --- ∑ (    Fc) V(v, l) ------- <n•l>
             *            N  h                   <n•h>
             *
             *
             *   Er() = f0 * DFV.x + f90 * DFV.y
             *
             */
            const float v = Visibility(NoV, NoL, linearRoughness) * NoL * (VoH / NoH);
            const float Fc = pow5(1 - VoH);
            r.x += v * (1.0f - Fc);
            r.y += v * Fc;
        }
    }
    return r * (4.0f / numSamples);
}

static float2 DFV_Multiscatter(float NoV, float linearRoughness, size_t numSamples) {
    float2 r = 0;
    const float3 V(std::sqrt(1 - NoV * NoV), 0, NoV);
    for (size_t i = 0; i < numSamples; i++) {
        const float2 u = hammersley(uint32_t(i), 1.0f / numSamples);
        const float3 H = hemisphereImportanceSampleDggx(u, linearRoughness);
        const float3 L = 2 * dot(V, H) * H - V;
        const float VoH = saturate(dot(V, H));
        const float NoL = saturate(L.z);
        const float NoH = saturate(H.z);
        if (NoL > 0) {
            const float v = Visibility(NoV, NoL, linearRoughness) * NoL * (VoH / NoH);
            const float Fc = pow5(1 - VoH);
            /*
             * Assuming f90 = 1
             *   Fc = (1 - V•H)^5
             *   F(h) = f0*(1 - Fc) + Fc
             *
             * f0 and f90 are known at runtime, but thankfully can be factored out, allowing us
             * to split the integral in two terms and store both terms separately in a LUT.
             *
             * At runtime, we can reconstruct Er() exactly as below:
             *
             *            4                <v•h>
             *   DFV.x = --- ∑ Fc V(v, l) ------- <n•l>
             *            N  h             <n•h>
             *
             *
             *            4                <v•h>
             *   DFV.y = --- ∑    V(v, l) ------- <n•l>
             *            N  h             <n•h>
             *
             *
             *   Er() = (1 - f0) * DFV.x + f0 * DFV.y
             *
             *        = mix(DFV.xxx, DFV.yyy, f0)
             *
             */
            r.x += v * Fc;
            r.y += v;
        }
    }
    return r * (4.0f / numSamples);
}

static float UTILS_UNUSED DFV_LazanyiTerm(float NoV, float linearRoughness, size_t numSamples) {
    float r = 0;
    const float cosThetaMax = (float) std::cos(81.7 * F_PI / 180.0);
    const float q = 1.0f / (cosThetaMax * pow6(1.0f - cosThetaMax));
    const float3 V(std::sqrt(1 - NoV * NoV), 0, NoV);
    for (size_t i = 0; i < numSamples; i++) {
        const float2 u = hammersley(uint32_t(i), 1.0f / numSamples);
        const float3 H = hemisphereImportanceSampleDggx(u, linearRoughness);
        const float3 L = 2 * dot(V, H) * H - V;
        const float VoH = saturate(dot(V, H));
        const float NoL = saturate(L.z);
        const float NoH = saturate(H.z);
        if (NoL > 0) {
            const float v = Visibility(NoV, NoL, linearRoughness) * NoL * (VoH / NoH);
            const float Fc = pow6(1 - VoH);
            r += v * Fc * VoH * q;
        }
    }
    return r * (4.0f / numSamples);
}

static float DFV_Charlie_Uniform(float NoV, float linearRoughness, size_t numSamples) {
    float r = 0.0;
    const float3 V(std::sqrt(1 - NoV * NoV), 0, NoV);
    for (size_t i = 0; i < numSamples; i++) {
        const float2 u = hammersley(uint32_t(i), 1.0f / numSamples);
        const float3 H = hemisphereUniformSample(u);
        const float3 L = 2 * dot(V, H) * H - V;
        const float VoH = saturate(dot(V, H));
        const float NoL = saturate(L.z);
        const float NoH = saturate(H.z);
        if (NoL > 0) {
            const float v = VisibilityAshikhmin(NoV, NoL, linearRoughness);
            const float d = DistributionCharlie(NoH, linearRoughness);
            r += v * d * NoL * VoH; // VoH comes from the Jacobian, 1/(4*VoH)
        }
    }
    // uniform sampling, the PDF is 1/2pi, 4 comes from the Jacobian
    return r * (4.0f * 2.0f * (float) F_PI / numSamples);
}

/*
 *
 * Importance sampling Charlie
 * ---------------------------
 *
 * Important samples are chosen to integrate DCharlie() * cos(theta) over the hemisphere.
 *
 * All calculations are made in tangent space, with n = [0 0 1]
 *
 *                      h (important sample)
 *                     /.
 *                    / .
 *                   /  .
 *                  /   .
 *         --------o----+-------> n
 *                   cos(theta)
 *                    = n•h
 *
 *  h is micro facet's normal
 *  l is the reflection of v around h, l = reflect(-v, h)  ==>  v•h = l•h
 *
 *  n•v is given as an input parameter at runtime
 *
 *  Since n = [0 0 1], we also have v.z = n•v
 *
 *  Since we need to compute v•h, we chose v as below. This choice only affects the
 *  computation of v•h (and therefore the fresnel term too), but doesn't affect
 *  n•l, which only relies on l.z (which itself only relies on v.z, i.e.: n•v)
 *
 *      | sqrt(1 - (n•v)^2)     (sin)
 *  v = | 0
 *      | n•v                   (cos)
 *
 *
 *  h = hemisphereImportanceSampleDCharlie()
 *
 *  l = reflect(-v, h) = 2 * v•h * h - v;
 *
 *  n•l = [0 0 1] • l = l.z
 *
 *  n•h = [0 0 1] • l = h.z
 *
 *
 *  pdf() = DCharlie(h) <n•h> |J(h)|
 *
 *               1
 *  |J(h)| = ----------
 *            4 <v•h>
 *
 *
 * Evaluating the integral
 * -----------------------
 *
 * We are trying to evaluate the following integral:
 *
 *                    /
 *             Er() = | fr(s) <n•l> ds
 *                    /
 *                    Ω
 *
 * For this, we're using importance sampling:
 *
 *                    1     fr(h)
 *            Er() = --- ∑ ------- <n•l>
 *                    N  h   pdf
 *
 * with:
 *
 *            fr() = DCharlie(h) V(v, l)
 *
 *
 *  It results that:
 *
 *            1                          4 <v•h>
 *    Er() = --- ∑ DCharlie(h) V(v, l) ------------ <n•l>
 *            N  h                     DCharlie(h) <n•h>
 *
 *
 *  +---------------------------------------+
 *  |          4             <v•h>          |
 *  |  Er() = --- ∑ V(v, l) ------- <n•l>   |
 *  |          N  h          <n•h>          |
 *  +---------------------------------------+
 *
 */
static float UTILS_UNUSED DFV_Charlie_IS(float NoV, float linearRoughness, size_t numSamples) {
    float r = 0.0;
    const float3 V(std::sqrt(1 - NoV * NoV), 0, NoV);
    for (size_t i = 0; i < numSamples; i++) {
        const float2 u = hammersley(uint32_t(i), 1.0f / numSamples);
        const float3 H = hemisphereImportanceSampleDCharlie(u, linearRoughness);
        const float3 L = 2 * dot(V, H) * H - V;
        const float VoH = saturate(dot(V, H));
        const float NoL = saturate(L.z);
        const float NoH = saturate(H.z);
        if (NoL > 0) {
            const float J = 1.0f / (4.0f * VoH);
            const float pdf = NoH; // D has been removed as it cancels out in the previous equation
            const float v = VisibilityAshikhmin(NoV, NoL, linearRoughness);

            r += v * NoL / (pdf * J);
        }
    }
    return r / numSamples;
}

void CubemapIBL::brdf(utils::JobSystem& js, Cubemap& dst, float linearRoughness) {
    CubemapUtils::process<CubemapUtils::EmptyState>(dst, js,
            [ & ](CubemapUtils::EmptyState&, size_t y,
                    Cubemap::Face f, Cubemap::Texel* data, size_t dim) {
                for (size_t x=0 ; x<dim ; ++x, ++data) {
                    const float2 p(Cubemap::center(x, y));
                    const float3 H(dst.getDirectionFor(f, p.x, p.y));
                    const float3 N = { 0, 0, 1 };
                    const float3 V = N;
                    const float3 L = 2 * dot(H, V) * H - V;
                    const float NoL = dot(N, L);
                    const float NoH = dot(N, H);
                    const float NoV = dot(N, V);
                    const float LoH = dot(L, H);
                    float brdf_NoL = 0;
                    if (NoL > 0 && LoH > 0) {
                        const float D = DistributionGGX(NoH, linearRoughness);
                        const float F = Fresnel(0.04f, 1.0f, LoH);
                        const float V = Visibility(NoV, NoL, linearRoughness);
                        brdf_NoL = float(D * F * V * NoL);
                    }
                    Cubemap::writeAt(data, Cubemap::Texel{ brdf_NoL });
                }
            });
}

void CubemapIBL::DFG(JobSystem& js, Image& dst, bool multiscatter, bool cloth) {
    auto dfvFunction = multiscatter ? DFV_Multiscatter : DFV;
    auto job = jobs::parallel_for<char>(js, nullptr, nullptr, uint32_t(dst.getHeight()),
            [&dst, dfvFunction, cloth](char const* d, size_t c) {
                const size_t width = dst.getWidth();
                const size_t height = dst.getHeight();
                size_t y0 = size_t(d);
                for (size_t y = y0; y < y0 + c; y++) {
                    Cubemap::Texel* UTILS_RESTRICT data =
                            static_cast<Cubemap::Texel*>(dst.getPixelRef(0, y));

                    const float h = (float) height;
                    const float coord = saturate((h - y + 0.5f) / h);
                    // map the coordinate in the texture to a linear_roughness,
                    // here we're using ^2, but other mappings are possible.
                    // ==> coord = sqrt(linear_roughness)
                    const float linear_roughness = coord * coord;
                    for (size_t x = 0; x < height; x++, data++) {
                        // const float NoV = float(x) / (width-1);
                        const float NoV = saturate((x + 0.5f) / width);
                        float3 r = { dfvFunction(NoV, linear_roughness, 1024), 0 };
                        if (cloth) {
                            r.b = float(DFV_Charlie_Uniform(NoV, linear_roughness, 4096));
                        }
                        *data = r;
                    }
                }
            }, jobs::CountSplitter<1, 8>());
    js.runAndWait(job);
}

} // namespace ibl
} // namespace filament